Geostatistical method for inferring RMR ahead of tunnel face excavation using dynamically exposed geological information

Accepted Manuscript Geostatistical method for predicting RMR ahead of tunnel face excavation using dynamically exposed geological information

Jianqin Chen, Xiaojun Li, Hehua Zhu PII: DOI: Reference:

S0013-7952(16)30849-3 doi: 10.1016/j.enggeo.2017.08.004 ENGEO 4613

To appear in:

Engineering Geology

Received date: Revised date: Accepted date:

21 December 2016 13 July 2017 1 August 2017

Please cite this article as: Jianqin Chen, Xiaojun Li, Hehua Zhu , Geostatistical method for predicting RMR ahead of tunnel face excavation using dynamically exposed geological information. The address for the corresponding author was captured as affiliation for all authors. Please check if appropriate. Engeo(2017), doi: 10.1016/j.enggeo.2017.08.004

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Geostatistical method for predicting RMR ahead of tunnel face excavation using dynamically exposed geological information Jianqin Chena, Xiaojun Lia,b,c,* , Hehua Zhua,c a

Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, 1239

Siping Road, Shanghai 200092, China State Key Laboratory for Disaster Reduction in Civil Engineering, Tongji University, 1239 Siping

PT

b

Road, Shanghai 200092, China c

College of Engineering, Tibet University, Tibet 850000, China

SC

RI

* Corresponding author

Abstract:

NU

Rock Mass Rating (RMR) is a rock mass classification system that is often used

MA

to select appropriate excavation methods and rock support systems in tunnel projects. This paper presents a geostatistical method for predicting RMR values quantitatively

ED

ahead of excavation of the tunnel face. The study makes full use of geological

EP T

information exposed on excavated tunnel faces to capture the spatial correlation of rock mass quality and later predicts the RMR value using the kriging method.

AC C

Predictions are constantly updated during the tunnel construction process. The advantages of the proposed method are as follows: (1) The RMR prediction uncertainty is quantified by accounting for spatial variability and model uncertainty; therefore, the resulting prediction can consider the geological conditions of the worst scenario; (2) The spatial variability of the geological condition is represented as a variogram model that is updated by observation data on the new excavated faces, and as the tunnel advances, the RMR prediction accuracy improves; and (3) The periodicity of geological conditions can be considered in RMR prediction. The 1

ACCEPTED MANUSCRIPT

proposed method is applied to a rock tunneling project, the Mingtang tunnel in Anhui province, China. The method achieves approximately 80% prediction accuracy; therefore, it has high potential as a tool for predicting RMR information ahead of

PT

tunnel face excavation.

RI

Keywords: Rock Mass Rating Prediction; Spatial Variability; Model Uncertainty;

AC C

EP T

ED

MA

NU

SC

Dynamic Updating

2

ACCEPTED MANUSCRIPT List of Tables Table. 1 Measured RMR values provided by the excavated tunnel faces during construction. Table 2. Variogram models parameters of the components using two kriging

PT

estimation methods.

RI

Table 3. Cross-validation criteria for two estimation methods.

SC

Table 4. Period of the variogram models for two kriging estimation methods using

NU

nested hole-effect model.

AC C

EP T

ED

MA

Table 5. Cross-validation criteria for two estimation methods considering periodicity.

3

ACCEPTED MANUSCRIPT List of Figures Fig. 1. Location of the Mingtang Tunnel in China. Fig. 2. (a) Three-dimensional geological model of the Mingtang Tunnel (granite is

PT

represented by cyan, and gneiss is represented by yellow); and (b) Positions of the measured RMR values provided by excavated tunnel faces.

RI

Fig. 3. Frequency distribution histogram of the Rock Mass Rating.

SC

Fig. 4. RMR estimation interval.

MA

circles are detailed in Table 1.

NU

Fig. 5. RMR prediction with increasing number of data; the actual positions of the

Fig. 6. Variogram fitting at four positions (K21+170, K21+422, K21+612 and

ED

K22+381); for each prediction position, three experimental variogram models are calculated with the number of data equal to eight, fourteen and twenty; for

EP T

each experimental variogram model, the bottom line represented by purple is fitted using the lower bound of  2f , the top line represented by green is fitted

AC C

using the upper bound of  2f and the middle line represented by blue is fitted by automatic fitting; (a) K21+170; (b) K21+422; (c) K21+612; and (d) K22+381.

Fig. 7. Kriging prediction interval considering model uncertainty; the blue line and black line are the estimation results using the upper and lower bounds of the variogram models; the red line is the true value calculated by the in situ tunnel face mapping; (a) K21+170; (b) K21+422; (c) K21+612; and (d) K22+381. 4

ACCEPTED MANUSCRIPT

Fig. 8. Kriging variance at four prediction positions; the blue line and black line are Kriging variance using the upper and lower bounds of the variogram models; (a) K21+170; (b) K21+422; (c) K21+612; and (d) K22+381. Fig. 9. Prediction interval at four prediction positions with consideration of model

PT

uncertainty and spatial variability; the blue line is the upper bound of the

RI

estimation, the black line is the lower bound of estimation and the red line is the

SC

true value calculated by in situ tunnel face mapping; (a) K21+170; (b) K21+422;

NU

(c) K21+612; and (d) K22+381.

Fig. 10. RMR prediction intervals; (a) Model One: model uncertainty; and (b) Model

MA

Two: model uncertainty and spatial variability.

Fig. 11. Variogram model fitting; (a) RMRbasis; (b) joint condition (JC); (c) RQD; (d)

ED

joint and bedding spacing (JS); (e) uniaxial compressive strength (UCS); and (f)

EP T

groundwater condition (GW).

Fig. 12. Scatterplot of the measured values and predicted values for two prediction

AC C

methods.

Fig. 13. Variogram model fitting considering periodicity; (a) RMRbasis; (b) joint condition (JC); (c) RQD; (d) joint and bedding spacing (JS); (e) uniaxial compressive strength (UCS); and (f) groundwater condition (GW). Fig. 14. Spectral analysis on RMRbasis and RMR components; (a) RMRbasis; (b) joint condition (JC); (c) RQD; (d) joint and bedding spacing (JS); (e) uniaxial compressive strength (UCS); and (f) groundwater condition (GW). 5

ACCEPTED MANUSCRIPT 1. Introduction Geological bodies

are

characterized

by

complexity,

uncertainty

and

heterogeneity after complex geological processes. Drilling investigations, in situ tests and surface geophysical investigations in the exploratory phase of a tunneling project

PT

are often limited and suffer from rough topography and difficult logistic conditions in

RI

mountainous terrain (Petronio et al., 2007). This may lead to large uncertainty and

SC

major geological risk in the construction phase of tunneling. Therefore, the rock mass

engineering design and construction.

NU

classification systems known as empirical methods are still used in current

MA

One of the most common classification systems is the Rock Mass Rating (RMR) that was developed by Bieniawski (1973, 1989). The RMR system incorporates six

ED

parameters including uniaxial compressive strength (UCS) of rocks, Rock Quality

EP T

Designation (RQD), joint and bedding spacing (JS), joint condition (JC), groundwater condition (GW), and orientation of discontinuities with respect to the opening axis.

AC C

Each of these parameters is classified and given a value or rating to express its influence on tunnel stability (Palmström, 2009). The RMR system is effective for evaluating the quality of rock mass surrounding the tunnel. It can be used to assist engineers in proposing suitable excavation techniques and support systems in the process of tunneling. Regarding RMR prediction ahead of tunnel excavation, it is well known that actual rock mass class may not be consistent with the design value because of the 6

ACCEPTED MANUSCRIPT

variability of rock mass quality. Therefore, predicting the RMR value ahead of tunnel excavation is crucial for selection of an appropriate rock support system, as well as to plan a safe and economically efficient excavation method. Typically, RMR prediction values range from 10 to 50 meters ahead of the tunnel face, and this is significant

PT

because about one week is needed to adjust the excavation method to the geological

RI

conditions ahead. Currently, tunnel engineers usually rely on the seismic prediction

SC

method or ground penetration radar to predict rock mass quality and hazardous

NU

geological structures ahead of tunnel excavation. However, the inspection results of the seismic prediction method or ground penetration radar are qualitative, so they

MA

must be interpreted by experienced geophysicists. The dynamic drill and blast excavation method provides a great opportunity for engineers and geologists to

ED

observe and sample the rock mass on the tunnel faces. The RMR values can be

EP T

calculated according to the exposed tunnel face. Moreover, tunnel construction is a dynamic process. As the tunnel advances, new geological information is added

AC C

constantly. However, this information has not been put into good use. A geostatistical approach has been applied to estimate RMR in different kinds of site characterizations. For tunnel site characterization, most current research uses drillhole data (Kaewkongkaew et al., 2015) to interpolate RMR values directly for the whole domain in the exploratory phase of tunnel project. Ferrari et al. (2014) performed the ordinary kriging estimation to a valley area. These studies focus mainly on the exploratory phase of a project. This means that the point-wise measurements 7

ACCEPTED MANUSCRIPT from drill holes or geological field surveys are used to interpolate the whole domain. As the estimation process is conducted only once, the variogram model is constant. However, new geological information could be gathered constantly as the tunnel advances. Therefore, the variogram model can be updated by using observation data

PT

obtained from newly excavated tunnel faces. In addition, model uncertainty that is

RI

introduced in the process of experimental variogram fitting is not considered in the

SC

conventional analysis, which generates prediction uncertainty. The lower bound of the

NU

RMR estimation value represents the possible worst scenario in tunneling. Tunnel excavation plan based on this result would be considered extremely conservative. The

MA

upper bound of RMR estimation value represents the best possible scenario, so tunnel excavation plans based on this result would be considered a high-risk venture. In

ED

addition, geological phenomenon often occurs repetitively over geological time,

EP T

leading to repetitive or cyclic variations in the facies and petrophysical properties (Pyrcz and Deutsch, 2014). The variogram fitting of RMR and its components must

AC C

also consider this periodicity.

In summary, this research focuses on applying geostatistics to quantitatively predict RMR values ahead of the excavation face during tunnel construction. This research offers several contributions: (1) The RMR prediction is based on dynamic observations of tunnel faces during excavation; (2) The RMR prediction uncertainty is quantified by taking into account model uncertainty and spatial variability; and (3) The periodicity of geological conditions can be considered in the RMR prediction. 8

ACCEPTED MANUSCRIPT

This paper is organized by introducing the background of a rock tunneling project in Section 2, presenting a geostatistical method for predicting RMR values quantitatively ahead of the tunnel face in Section 3, discussing the application of the

PT

method in Section 4, and presenting several conclusions in Section 5.

RI

2. Site Description

SC

To demonstrate application of the geostatistical method for RMR prediction for

NU

tunnel engineering, the Mingtang Tunnel is selected as a case study. Mingtang Tunnel is a rock tunnel situated in Yuexi County, Anhui Province, China (Fig. 1). This is part

MA

of the Yue-wu highway that connects Wuhan city and Shanghai city. The tunnel was excavated using the drill and blast method. The total length of the tunnel is 7.548 km,

ED

and the maximum burial depth is approximately 562 m. As shown in Fig. 2(a), the

EP T

geology of the site area is mainly composed of granite and gneiss. Three major faults transect rock formations in the tunnel alignment. The study area is in the interval

AC C

between K17+599 and K24+708. During tunnel construction, geological information is sampled from the excavation faces in forty-three discrete positions, and their corresponding rock mass quality is evaluated. The positions and measured RMR values are shown in Table 1. Fig. 2(b) shows the positions of the tunnel face relative to where the in situ tunnel face mapping is obtained.

3. Methodology 9

ACCEPTED MANUSCRIPT

3.1 Photogrammetry to obtain RMR values Dynamic drill and blast excavation provides a useful opportunity for engineers and geologists to observe and sample rock mass. Most difficult, however, is obtaining the discontinuity information safely and efficiently. In our study, RQD value, joint and

PT

bedding spacing (JS) and joint condition (JC) (except infilling and weathering) of

RI

actual RMR values can be obtained using the photogrammetry-based tunnel face

SC

mapping technique (Haneberg, 2008; Sturzenegger et al., 2009; Chen et al., 2016; Li

NU

et al., 2016; Zhu et al., 2016). Uniaxial compressive strength (UCS) of rock mass is measured through the uniaxial compression test. Groundwater condition (GW),

MA

infilling and weathering are obtained by field observations. 3.2 Geostatistical method

ED

Adjacent areas show similar geological conditions, so the RMR data sampled

EP T

from locations closer together tend to be more similar than those far apart. This spatial correlation can be quantified through a variogram, which is a measure of “geological

AC C

variability”. Variogram increases as data become more dissimilar. Therefore, the parameters of the variogram best fitting the data are first conducted using R software (R. Core Team, 2015) with the GSTAT package (Pebesma, 2004). Then, Kriging is used to provide an unbiased estimation, with uncertain quantification. 3.2.1 Exploratory spatial data analysis The RMR values range from 44 to 74 with a mean of 61.5 and a standard deviation of 7.77. The frequency distribution is shown in Fig. 3. The normality of 10

ACCEPTED MANUSCRIPT RMR distribution is verified using the Shapiro-Wilk test (Shapiro and Wilk, 1965). The Gaussian distribution of RMR is confirmed with a significance level of 5%. Moreover, trend analysis results show that there is no trend in the RMR. This confirms the stationarity hypothesis of RMR in the studied domain.

PT

3.2.2 Prediction with model uncertainty and spatial variability

RI

3.2.2.1 Model uncertainty

SC

Construction of the semivariogram, which is a mathematical model that captures

NU

the spatial correlation among data, is a very important step in any geostatistical analysis. The semivariogram is a measure of variability that increases as samples

MA

become more dissimilar. The variogram is defined as the expected value of a squared difference (Isaaks and Srivastava, 1989):

ED

2  h   Var[Z ( x)  Z ( x  h)]  E{[ Z ( x)  Z ( x  h)]2}

(1)

EP T

where Z is a stationary random function with known mean m and variance  2 , which is independent of location, so m( x)  m and  2 ( x )   2 for all locations x

AC C

in the study area.

Spatial correlation is assessed by calculating the experimental variogram, then fitting with different parametric variogram models, such as the Spherical model, Exponential model, and Gaussian model. There are four main approaches for estimating the parameters of the variogram model: visual, least square regression (Trangmar et al., 1985), weighted least squares (Cressie, 1985), and likehood methods (Cressie, 1991). However, the process of variogram model fitting introduces some 11

ACCEPTED MANUSCRIPT

degree of smoothing and subjectivity. These elements can introduce parametric uncertainty (Rubin, 2003). The method of least squares is used to fit the rising limb of the parametric variogram to the experimental variogram, while holding the sill value equal to the

PT

sample variance, so that range can be calculated. The sample variance is treated as a

RI

random variable. Hence, the uncertainty in determining the range is also considered

SC

by least square regression. Estimation of the sample variance is subject to uncertainty, which is estimated by Rehfeldt et al. (1992):

is the sample variance and

MA

where ˆ 2f

2 4 ˆ f N*

NU

Var (ˆ 2f ) 

(2)

N * is the number of independent

measurements that can be obtained by dividing the tunnel into independent segments.

ED

When assuming the sample variance to be chi-square distributed, the statistical

EP T

distribution of the estimation error is normally distributed and an approximately 95% confidence region can be constructed about  2f , as (3)

AC C

ˆ 2f  2[Var(ˆ 2f )]1/2   2f  ˆ 2f  2[Var(ˆ 2f )]1/2

Two theoretical variogram models can be fitted by using the lower and upper limits of  2f , and determining the corresponding ranges automatically. Kriging estimation of these two theoretical variogram models generates the corresponding upper and lower bounds of RMR values ( RMRupper and RMRlower ) that accounts for the model uncertainty. 3.2.2.2 Spatial variability 12

ACCEPTED MANUSCRIPT

In geostatistical estimation theory, the associated standard error for each estimation is a measure of the possible variation in value, and hence the uncertainty of the estimated value. If the standard error is small, then the possible statistical variation is corresponding small; if it is large, then the possible variation is large.

PT

When only considering the spatial variability, the estimation interval of RMR can

RI

be expressed as following with 68% confidence:

(4)

SC

RMR   RMR  RMRPred  RMR   RMR

NU

where RMR is the mean and  RMR is standard deviation of the kriging estimation. This prediction result only considers the spatial variability.

MA

3.2.2.3 Model uncertainty and spatial variability

As is shown in Fig. 4, when considering both the model uncertainty and spatial

ED

variability, the estimation interval of RMR can be expressed as following: (5)

EP T

min{RMRlower   lower , RMRupper   upper }  RMRPred  max{RMRlower   lower , RMRupper   upper }

where RMRupper and RMRlower are the kriging estimated values using the bounds of

AC C

theoretical variogram models (defined in Eq. 3), and  upper and  lower are the corresponding standard errors. 3.2.3 RMR prediction The spatial variability can be updated using observation data on newly excavated tunnel faces. To evaluate the dynamic update effect of RMR prediction, different number of observation data were used to construct the variogram model. Then, the prediction results at the same position were compared. 13

ACCEPTED MANUSCRIPT

As illustrated in Fig. 5, red circles represent the tunnel positions where predictions are made and green circles represent the excavated tunnel positions with field measured RMR values. The actual positions of the circles are detailed in Table 1. The positions of the circles relative to the tunnel are shown in Fig. 2(b). There are

PT

three clouds with different colors next to each of the red circles, and each cloud

RI

contains different numbers of observation data ranging from eight, fourteen and

SC

twenty. Three variogram models are constructed based on the observation data in the

NU

clouds and they predict the RMR values at the position of the red circle. The goal is to show how the RMR prediction changes as the number of measurements increase. The

MA

prediction position advances to different red circles, such as tunnel excavation. For each prediction, the number of data surrounded by the clouds is the same.

ED

3.2.4 Two kriging estimation methods: Direct Method and Indirect Method

EP T

Two types of prediction strategies are conducted in our study. The first approach is called the Direct Method. In this approach, the base RMR value is fit with a

AC C

variogram model and then estimated directly. The base RMR is composed of five parameters:

RMRbasis  UCS  RQD  JS  JC  GW

(6)

The last parameter (orientation of discontinuities with respect to the opening axis) is excluded from the kriging estimation because it depends on the characteristics of rock discontinuities, as well as their relation to the tunnel structure, which is unknown (Pinheiro et al., 2016). 14

ACCEPTED MANUSCRIPT

The second approach is called the Indirect Method. In this approach, parameters composed of RMR are estimated separately as individual variables. The final estimated RMR is obtained by adding the ratings of each variable. The kriging 2 estimation Z ( RMR2 ) and kriging variance  RMR are represented as 2

PT

Z ( RMR2 )  Z (UCS )  Z ( RQD)  Z ( JS )  Z ( JC )  Z (GW ) 2 2 2 2 2  RMR   UCS   RQD   JS2   JC   GW 2

(8)

SC

RI

2cov(UCS , RQD )  2cov(UCS , JS )  2cov(UCS , JC )  2cov(UCS , GW ) 2cov( RQD, JS )  2cov( RQD, JC )  2cov( RQD, GW ) 2cov( JS , JC )  2cov( JS , GW )  2cov( JC , GW )

(7)

4. Results and Discussion

MA

NU

2 2 2 2 2 where  UCS ,  RQD ,  JS ,  JC and  GW are the kriging variances for each component.

ED

4.1 Model uncertainty and spatial variability

EP T

Fig. 6 shows the fitted results of variogram models at four positions, respectively. For each prediction position, three experimental variogram models are calculated with

AC C

the number of observation data equaling to eight, fourteen and twenty. The aim of these models is to investigate how prediction accuracy changes with the number of observation data.

In addition, the experimental variogram model is fitted with three theoretical variogram models. The bottom line (represented by purple) is generated using the lower bound of  2f , the top line (represented by green) is generated using the upper bound of  2f , and the middle line (represented by blue) is generated by automatic fitting. The bottom line and top line almost encapsulate the variation range of the sill 15

ACCEPTED MANUSCRIPT

of the variogram model. As illustrated in Fig. 6, the variogram models change with increasing measurement data. In Fig. 6(a), the sills of the variogram models decrease gradually. The sills in Fig. 6(b) and Fig. 6(c) increase first and decrease later. In contrast, the

PT

sills increase gradually in Fig. 6(d). Therefore, with newly added data, the variogram

RI

model is more consistent with the actual situation, meaning that the uncertainty of the

SC

variogram model decreases.

NU

Fig. 7 shows the Kriging estimation result with a different number of observation data. The blue line and black line are the estimation results using the upper and lower

MA

bounds of the variogram models. The red line is the true value calculated by the in situ tunnel face mapping. The prediction interval considers the variogram model

ED

uncertainty described in Section 3.2.2. The prediction interval covers the true value

EP T

gradually with an increase of the measured data, which indicates that prediction accuracy increases.

AC C

Fig. 8 shows the corresponding kriging variance using upper and lower bounds of variogram models. The upper and lower bounds of kriging variance are used to quantify spatial variability. The final RMR estimation interval allows for both model uncertainty and spatial variability, as shown in Fig. 9. When the spatial variability is incorporated into the prediction result, the prediction interval increases to cover the true value. The prediction result of the eight positions is summarized in Fig. 10. Model One 16

ACCEPTED MANUSCRIPT

refers to prediction that considers only model uncertainty. Model Two refers to prediction that considers both model uncertainty and spatial variability. Fig. 10(a) shows the RMR prediction intervals for Model One. Fig. 10(b) shows the prediction intervals for Model Two. In Fig. 10, the prediction intervals decrease with an increase

PT

of measurement data, especially at the position K21+170. For the position K21+170,

RI

when the number of measurement data is eight, the upper bound of  2f is large,

SC

which leads to a large RMR estimation interval with the lower bound of the estimated

NU

RMR value close to zero. When the number of measurement data increases to twenty, the upper bound of  2f decreases, which leads to a smaller RMR estimation interval

MA

with the lower bound of estimated RMR value approximating to the true value. Reduction of the prediction interval does not decrease the prediction accuracy.

ED

4.2 Comparison of the Direct Method and Indirect Method

EP T

For the Indirect Method, the components of RMR are estimated separately as individual variables and fitted with separate variogram models from Fig. 11(b) to Fig.

AC C

11(f). The parameters of these variogram models are listed in Table 2. For the variogram model of each component, the range of GW is the largest, the ranges of UCS are secondary, and the ranges of JC, RQD and JS are the smallest. Because ground water can flow in the fractures of rock mass, the influence area is larger when comparing the fracture parameters (JC, RQD and JS) and strength parameter (UCS). However, the fractures are generated by tectonic actions and thus show larger variability and heterogeneities in the domain. Moreover, the sill of the GW is largest 17

ACCEPTED MANUSCRIPT

because three faults traverse the Mingtang tunnel and dramatically change the groundwater conditions. There is abundant ground water neighboring the faults, while there is limited ground water in the regions far away from the tunnel. In the Direct Method, the basic RMR value (defined in equation 6) is fitted with

PT

a variogram model, as is shown in Fig. 11(a). The sill of RMRbasis (67.00) is affected

RI

by the largest sill of its component, which is GW with the value of 55.55. In addition,

SC

the range of RMRbasis (171.33) is affected by the smallest sill of its component,

NU

which is JC with a value of 175.45.

Fig. 12 shows the scatterplot between measured values and predicted values for

MA

two estimation methods. It reveals a near- linear positive relationship between the measured value and predicted value for the two models. The correlation coefficients

ED

for the Indirect Method and Direct Method are 0.73 and 0.7, respectively. Three

EP T

performance metrics, the mean error (ME), the root mean squared error (RMSE) and the mean standardized prediction error (MSPE), are used to evaluate the overall

AC C

predictive accuracy. The ME should be close to zero. The RMSE measures the difference between the predicted values and measured values, and this should be as small as possible. The MSPE is standardized by dividing ME by the kriging variance. An accurate model should have an MSPE close to zero. As shown in Table 3, the Indirect Method has the smaller absolute values of ME and RMSE, with the MSPE value being closer to zero. Therefore, the Indirect Method has better estimation accuracy for the case study. 18

ACCEPTED MANUSCRIPT

4.3 Periodicity in the variogram models As is shown in Fig. 13, the experimental variograms exhibit a sinusoidal wave that conveys the cyclicity of the underlying geological phenomenon. This phenomenon is sometimes referred to as the hole effect (Journel and Huijbregts, 1978;

PT

Webster and Oliver, 2007). Hole effect is representative of a “periodic” phenomenon

RI

and is characterized by undulations in the variogram. Therefore, nested variograms

SC

composed of the Spherical model, Periodic model and Hole-effect model are

NU

constructed to fit the experimental variogram models of the whole RMR and each component of RMR.

MA

Spectral analysis (Lomb, 1976) is used to detect cyclicities in the RMRbasis and RMR component. As shown in Fig. 14(a), the power spectral density of RMRbasis is

ED

dominated by three significant peaks at frequencies of 0.0003, 0.0006 and 0.0021 that

EP T

correspond to periods of 1421.8, 507.8 and 122.6 meters, respectively. In the RMR component, there are also strong cyclicities indicated by peaks in the power spectral

AC C

density. As shown in Table 4, the periods computed by the variogram model are similar to the periods computed by the spectral analysis. This cyclicity suggests a strong influence of tectonic actions of the geological body on the quality of rock mass. Cross validation is performed to test the fit of the nested models to the data. Their correlation coefficients both increase to 0.76. Three performance metrics are used to evaluate the overall predictive accuracy. As is shown in Table 5, the Indirect and Direct Methods have almost similar estimation accuracies for ME, RMSE and 19

ACCEPTED MANUSCRIPT

MSPE values after considering the periodicity in the variogram models. Furthermore, prediction results of the Indirect Method, with the exception of RMSE, is better after incorporating the periodicity model, in comparison to Table 3. This indicates that integrating periodicity into the variogram increases prediction accuracy. However,

PT

when the Indirect Method considers periodicity, the variability of the prediction

SC

RI

increases as each component is predicted separately.

NU

5. Conclusions

This paper presents a geostatistical method for predicting RMR values

MA

quantitatively ahead of tunnel face excavation by making full use of the geological information collected from the tunnel face and considering the dynamic construction

ED

process of tunneling. RMR prediction accuracy is evaluated by considering spatial

EP T

variability and model uncertainty, newly added geological data, and periodicity of the geological condition.

AC C

First, RMR prediction uncertainty is quantified by accounting for spatial variability and model uncertainty. Model uncertainty is derived from the variogram model fitting process, and spatial variability originates from variation in the estimated value. The resulting RMR prediction interval incorporates the best scenario case and worst scenario case. The tunnel excavation method is based on the limits of the interval representing different levels of risk. Furthermore, considering the model and spatial variability also increases the prediction accuracy. 20

ACCEPTED MANUSCRIPT

Second, RMR prediction accuracy increases with newly added geological data. Model uncertainty decreases with newly added geological data. Finally, RMR prediction is implemented in two approaches: the Direct Method and Indirect Method. The Direct Method estimates the RMR value directly, while the

PT

Indirect Method estimates the components constituting the RMR system separately. In

RI

our case study, the Indirect Method shows better estimation accuracy than the Direct

SC

Method. Moreover, the nested variograms composed of Spherical model, Periodic

NU

model and Hole-effect models perform relatively well in terms of prediction accuracy for both the Indirect Method and Direct Method. When the Indirect Method considers

Acknowledgments

ED

MA

periodicity, variability increases because each component is predicted separately.

EP T

The authors gratefully acknowledge support from the Natural Science Foundation of China (NSFC 41272289, 41130751), the Science and Technology Plan

AC C

Project of the Ministry of Transport of China (2013318J02120), and the Fundamental Research Funds for Central Universities. Special thanks are due to Prof. Yoram Rubin from University of California, Berkeley for his valuable instructions regarding the geostatistical method.

References Alimoradi, A., Moradzadeh, A., Naderi, R., Salehi, M.Z., Etemadi, A., 2008. 21

ACCEPTED MANUSCRIPT

Prediction of geological hazardous zones in front of a tunnel face using TSP-203 and artificial neural networks. Tunnelling and Underground Space Technology 23(6), 711-717. Amberg Measuring Technique, 2002. TSP 203 Processing. Amberg Co.,Switzerland.

PT

Bieniawski, Z.T., 1973, Engineering classification of jointed rock masses: Transaction

RI

of the South African Institution of Civil Engineers, v. 15, pp. 335-344.

SC

Bieniawski, Z.T., 1989, Engineering rock mass classifications: a complete manual for

NU

engineers and geologists in mining, civil, and petroleum engineering, New York, Wiley, xii, 251pp.

MA

Cressie, N., 1985. Fitting variogram models by weighted least squares.Journal of the International Association for Mathematical Geology, 17(5), 563-586.

ED

Cressie, N., 1991. Statistics for spatial data. John Wiley & Sons, NY, pp.900.

EP T

Chen, J., Zhu, H., Li, X., 2016. Automatic extraction of discontinuity orientation from rock mass surface 3D point cloud. Computers & Geosciences 95, 18-31.

AC C

Choi, J.Y., Lee, C.I., 2007. An estimation of rock mass rating using 3D- indicator kriging approach with uncertainty assessment of rock mass classification. In: Proceedings of the 11th congress of the international society for rock mechanics, Lisbon, vol. 2, pp. 1285-88. Choi, Y., Yoon, S.Y., Park, H.D., 2009. Tunneling analyst: a 3D GIS extension for rock mass classification and fault zone analysis in tunneling. Computers & Geosciences 35(6), 1322-1333. 22

ACCEPTED MANUSCRIPT

Exadaktylos, G., Stavropoulou, M., 2008. A specific upscaling theory of rock mass parameters exhibiting spatial variability: analytical relations and computational scheme. International Journal of Rock Mechanics and Mining Sciences 45(7), 1102-1125.

PT

Ferrari, F., Apuani, T., Giani, G.P., 2014. Rock Mass Rating spatial estimation by

RI

geostatistical analysis. International Journal of Rock Mechanics and Mining

SC

Sciences 70, 162-176.

NU

Gholami, R., Rasouli, V., Alimoradi, A., 2013. Improved RMR rock mass

engineering 46(5), 1199-1209.

MA

classification using artificial intelligence algorithms. Rock mechanics and rock

Haneberg, W.C., 2008. Using close range terrestrial digital photogrammetry for 3-D

ED

rock slope modeling and discontinuity mapping in the United States. Bulletin of

EP T

Engineering Geology and the Environment 67(4), 457-469. Isaaks, E.H., Srivastava, R.M., 2001. An introduction to applied geostatistics. 1989.

AC C

New York, USA: Oxford University Press. Jones DR, A taxonomy of global optimization methods based on response surfaces. Journal of Global Optimization, 23, 345-383.

Journel, A.G, Huijbregts, C.J., 1978. Mining Geostat. NewYork, NY: Academic Press, Inc. Kaewkongkaew, K., Phien-Wej, N., Kham-ai, D., 2015. Prediction of rock mass along tunnels by geostatistics. KSCE Journal of Civil Engineering 19(1), 81-90. 23

ACCEPTED MANUSCRIPT

Lomb, N.R., 1976. Least-squares frequency analysis of unequally spaced data. Astrophysics and space science, 39(2), 447-462. La Pointe, P.R., 1980, January. Analysis of the spatial variation in rock mass properties through geostatistics. In The 21st US Symposium on Rock Mechanics

PT

(USRMS). American Rock Mechanics Association.

RI

Palmström, A., 2009. Combining the RMR, Q, and RMi classification systems.

SC

Tunnelling and Underground Space Technology, 24(4), 491-492.

NU

Pebesma, E.J., 2004. Multivariable geostatistics in S: the gstat package. Computers & Geosciences, 30(7), 683-691.

MA

Pyrcz, M.J. and Deutsch, C.V., 2003. The whole story on the hole effect.Geostatistical Association of Australasia, Newsletter, 18, pp.3-5.

ED

Pyrcz, M.J. and Deutsch, C.V., 2014. Geostatistical reservoir modeling. Oxford

EP T

university press, New York, NY, pp.87. R Core Team: R: A Language and Environment for Statistical Computing, R for

Statistical

Computing,

Vienna,

Austria,

available

at:

AC C

Foundation

https://www.R-project.org/ (last access: 16 December 2015), 2015. Santos, V., da Silva, A.P.F., Brito, M.G., 2015. Prediction of RMR Ahead Excavation Front

in D&B

Tunnelling.

In Engineering

Geology

for Society and

Territory-Volume 6 (pp. 415-419). Springer International Publishing. Stavropoulou, M., Exadaktylos, G., Saratsis, G., 2007. A combined three-dimensional geological- geostatistical- numerical model of underground excavations in rock. 24

ACCEPTED MANUSCRIPT

Rock mechanics and rock engineering 40(3), 213-243. Trangmar, B.B., Yost, R.S. and Uehara, G., 1986. Application of geostatistics to spatial studies of soil properties. Advances in agronomy, 38, 45-94. You, K., Lee, J. S., 2006. Estimation of rock mass classes using the 3-dimensional

PT

multiple indicator kriging technique. Tunnelling and Underground Space

RI

Technology 21(3), 229-229.

deformation

analysis

(DDA)

with

binocular

NU

dimensional discontinuous

SC

Zhu, H., Wu, W., Chen, J., Ma, G., Liu, X., Zhuang, X., 2016. Integration of three

photogrammetry for stability analysis of tunnels in blocky rockmass. Tunnelling

MA

and Underground Space Technology 51, 30-40.

Li, X., Chen, J., Zhu, H., 2016. A new method for automated discontinuity trace

ED

mapping on rock mass 3D surface model. Computers & Geosciences 89, 118-131.

excavation

EP T

Petronio, L., Poletto, F., Schleifer, A., 2007. Interface prediction ahead of the front

by

the

tunnel-seismic-while-drilling

(TSWD)

method.

AC C

Geophysics 72(4), 39-44.

Egaña, M., Ortiz, J.M., 2013. Assessment of RMR and its uncertainty by using geostatistical simulation in a mining project. Journal of GeoEngineering 8(3), 83-90. Rehfeldt, K.R., Boggs, J.M., Gelhar, L.W., 1992. Field study of dispersion in a heterogeneous aquifer: 3. Geostatistical analysis of hydraulic conductivity. Water Resources Research 28(12), 3309-3324. 25

ACCEPTED MANUSCRIPT

Sturzenegger, M., Stead, D., 2009. Quantifying discontinuity orientation and persistence on high mountain rock slopes and large landslides using terrestrial remote sensing techniques. Natural Hazards and Earth System Science 9(2), 267-287.

PT

Pinheiro, M., Vallejos, J., Miranda, T., Emery, X., 2016. Geostatistical simulation to

SC

mass rating. Engineering Geology 205, 93-103.

RI

map the spatial heterogeneity of geomechanical parameters: A case study with rock

NU

Webster, R., Oliver, M.A., 2007. Geostatistics for environmental scientists. John Wiley & Sons.

MA

Zhang, L., Einstein, H.H., 1998. Estimating the mean trace length of rock discontinuities. Rock Mechanics and Rock Engineering 31(4), 217-235.

ED

Rubin, Y., 2003. Applied stochastic hydrogeology. Oxford University Press, New

Shaphiro,

S.S.,

EP T

York, NY, 56pp. Wilk,

M.B.,

1965.

An

AC C

normality. Biometrika, 52(3), 591-611.

26

analysis

of

variance

test

for

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC C

EP T

ED

MA

Fig. 1. Location of the Mingtang Tunnel in China.

Fig. 2. (a) Three-dimensional geological model of the Mingtang Tunnel (granite is represented by cyan, and gneiss is represented by yellow); and (b) Positions of the measured RMR values provided by excavated tunnel faces.

27

SC

RI

PT

ACCEPTED MANUSCRIPT

AC C

EP T

ED

MA

NU

Fig. 3. Frequency distribution histogram of the Rock Mass Rating.

Fig. 4. RMR estimation interval.

28

ACCEPTED MANUSCRIPT

EP T

ED

MA

NU

SC

RI

the circles are detailed in Table 1.

PT

Fig. 5. RMR prediction with increasing number of data; the actual positions of

AC C

Fig. 6. Variogram fitting at four positions; (a) K21+170; (b) K21+422; (c) K21+612; and (d) K22+381.

29

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Fig. 7. Kriging prediction interval considering model uncertainty; (a) K21+170;

AC C

EP T

ED

MA

(b) K21+422; (c) K21+612; and (d) K22+381.

Fig. 8. Kriging variance at four prediction positions; (a) K21+170; (b) K21+422; (c) K21+612; and (d) K22+381.

30

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Fig. 9. Prediction interval at four prediction positions with consideration of model uncertainty and spatial variability; (a) K21+170; (b) K21+422; (c) K21+612;

AC C

EP T

ED

MA

and (d) K22+381.

31

EP T

ED

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

AC C

Fig. 10. RMR prediction intervals; (a) Model One: model uncertainty; and (b) Model Two: model uncertainty and spatial variability.

32

RI

PT

ACCEPTED MANUSCRIPT

SC

Fig. 11. Variogram model fitting; (a) RMRbasis; (b) joint condition (JC); (c) RQD; (d) joint and bedding spacing (JS); (e) uniaxial compressive strength (UCS);

AC C

EP T

ED

MA

NU

and (f) groundwater condition (GW).

Fig. 12. Scatterplot of the measured values and predicted values for two prediction methods.

33

RI

PT

ACCEPTED MANUSCRIPT

SC

Fig. 13. Variogram model fitting considering periodicity; (a) RMRbasis; (b) joint condition (JC); (c) RQD; (d) joint and bedding spacing (JS); (e) uniaxial

AC C

EP T

ED

MA

NU

compressive strength (UCS); and (f) groundwater condition (GW).

Fig. 14. Spectral analysis on RMRbasis and RMR components; (a) RMRbasis; (b) joint condition (JC); (c) RQD; (d) joint and bedding spacing (JS); (e) uniaxial compressive strength (UCS); and (f) groundwater condition (GW).

34

ACCEPTED MANUSCRIPT Table. 1 Measured RMR values provided by the excavated tunnel faces during construction.

AC C

PT

RI

64 56 61 59 59 54 63 69 66 66 62 68 68 74 60 74 74 74 72 74 60 50 60 57 57 64 64 44 45 47 54 52 61 59 57 52 58 68

SC NU

K17+599 K17+620 K17+640 K17+681 K19+676 K20+124 K20+205 K20+241 K20+347 K20+375 K20+402 K20+432 K20+456 K20+475 K20+583 K20+650 K20+880 K20+905 K20+930 K20+960 K20+977 K21+037 K21+062 K21+170 K21+200 K21+405 K21+422 K21+503 K21+547 K21+612 K21+955 K21+987 K22+381 K22+533 K22+575 K22+881 K22+939 K23+790

Measured RMR value

ED

EP T

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

Position

MA

Number

35

ACCEPTED MANUSCRIPT 39 40 41 42 43

K24+090 K24+120 K24+149 K24+177 K24+708

65 65 65 65 59

SC

RI

PT

Table 2. Variogram models parameters of the components using two kriging estimation methods. Parameters in Direct Indirect Method Method estimation Methods JC RQD JS UCS GW Model parameters RMRbasis Sph*

Sph*

Sph*

Sph*

Sph*

Sph*

Sill

67.00

16.55

1.81

5.17

1.04

55.55

Range

171.33

175.45

260

163.68

350

1107.85

Nugget

0

0

0

0

0

0

MA

NU

Model

ED

*Sph stands for Spherical model.

EP T

Table 3. Cross-validation criteria for two estimation methods. Method Indirect Method Direct Method Criterion ME 0.0558 -0.2035

AC C

RMSE MSPE

1.1179

5.4462

0.0038

-0.0082

Table 4. Period of the variogram models for two kriging estimation methods using nested hole-effect model. Direct Parameters in Indirect Method Method estimation Methods RMRbasis JC RQD JS UCS GW Period (m) Variogram model

225

280

500

36

450

700

1500

ACCEPTED MANUSCRIPT Spectral analysis

123

284

711

296

592

1185

Table 5. Cross-validation criteria for two estimation methods considering periodicity. Method Indirect Method Direct Method Criterion ME 0.0125 -0.0865 5.0961

4.9934

MSPE

0.0011

-0.0051

AC C

EP T

ED

MA

NU

SC

RI

PT

RMSE

37

ACCEPTED MANUSCRIPT

Highlights 

>A geostatistical method for predicting RMR values quantitatively ahead of the tunnel face which makes full use of geological information exposed on tunnel face is presented.



>The spatial variability of geological condition is updated by observation data



PT

on the new excavated faces. The RMR prediction uncertainty is quantified by accounting for spatial variability and model uncertainty.

EP T

ED

MA

NU

SC

RI

The prediction of RMR values allows for periodicity of geological conditions.

AC C



38