Estimation of deformation modulus of rock masses based on Bayesian model selection and Bayesian updating approach

Engineering Geology 199 (2015) 19–27

Contents lists available at ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

Estimation of deformation modulus of rock masses based on Bayesian model selection and Bayesian updating approach Xianda Feng a,1, Rafael Jimenez b,⁎ a b

School of Civil Engineering and Architecture, University of Jinan, No. 336, West Road of Nan Xinzhuang Xilu Jinan 250022 Shandong Province, PR China Technical University of Madrid, ETSI Caminos, C. y P; C/Profesor Aranguren s/n; 28040 Madrid; Spain

a r t i c l e

i n f o

Article history: Received 14 April 2015 Received in revised form 20 August 2015 Accepted 8 October 2015 Available online 14 October 2015 Keywords: Deformation modulus Predictive uncertainty Bayesian updating approach Model selection Rock mass

a b s t r a c t The deformation modulus is one of the most important parameters to model the behavior of rock masses, but its direct measurement by in situ tests is costly, time-consuming and sometimes infeasible. For that reason, many models have been proposed to estimate the deformation moduli of rock masses based on geotechnical classification indices, such as the Rock Mass Rating (RMR), the Geological Strength Index (GSI), the Tunneling Quality Index (Q), or the Rock Quality Designation (RQD). We present an approach, based on model selection criteria —such as Akaike information criterion (AIC), Bayesian information criterion (BIC) and deviance information criterion (DIC)— to select the most appropriate model, among a set of four candidate models —linear, power, exponential and logistic—, to estimate the deformation modulus of a rock mass, given a set of observed data. Once the most appropriate model is selected, a Bayesian framework is employed to develop predictive distributions of the deformation moduli of rock masses, and to update them with new project-specific data that significantly reduce the associated predictive uncertainty. Such Bayesian updating approach can, therefore, affect our computed estimates of probability of failure, which is of significant interest to reliability-based rock engineering design. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The deformation modulus of a rock mass is one of the most important parameters influencing its mechanical behavior (Kayabasi et al., 2003). Therefore, it is one of the most commonly required input parameters for both analytical and numerical methods in geotechnical engineering, and it is the basis for many geotechnical analyses (Palmström and Singh, 2001). According to the Commission of Terminology of the International Society of Rock Mechanics (ISRM, 1975), the deformation modulus of a rock mass, Erm, is defined as “the ratio of stress to corresponding strain during loading of a rock mass including elastic and inelastic behavior”. There are many in situ tests to directly measure Erm such as plate loading or plate jacking tests, dilatometer tests, flat jack tests, pressure chamber tests, etc. (Bieniawski, 1978; Palmström and Singh, 2001; Hoek and Diederichs, 2006; Kang et al., 2012). However, such in situ tests are time consuming, expensive and sometimes even infeasible. For that reason, many models have been proposed to indirectly estimate the deformation moduli of rock masses based on geotechnical classification indices such as the Rock Mass Rating (RMR), the Geological Strength Index (GSI), the Tunneling Quality Index (Q), or the Rock Quality

⁎ Corresponding author. E-mail address: [email protected] (R. Jimenez). 1 Formerly at Technical University of Madrid, Spain.

http://dx.doi.org/10.1016/j.enggeo.2015.10.002 0013-7952/© 2015 Elsevier B.V. All rights reserved.

Designation (RQD) (Bieniawski, 1978; Serafim and Pereira, 1983; Nicholson and Bieniawski, 1990; Barton, 1996; Sonmez et al., 2004; Hoek and Diederichs, 2006; Zhang, 2010). The most commonly used models are summarized in Table 1. Since dozens of models (see Table 1) have been proposed during the past decades, engineering practitioners are often concerned about “how to choose the most appropriate model” to be used in a particular project for which there is a set of available data. To tackle this problem, a suitable model selection method should be adopted (see e.g., Burnham and Anderson (2002) for details). Since most of the models (i.e., empirical correlations) listed in Table 1 are obtained using regression methods, one may naturally think that comparing the R2 of different models is the most convenient approach for model selection. But this approach only measures the goodness of fit of the model; model complexity is ignored, hence always favoring “fuller” models with more parameters. In other words, “neglecting the principle of parsimony makes it a poor technique for model selection” (Johnson and Omland, 2004). In contrast, model selection criteria based on “information criteria” consider both model fit and complexity, hence being more suitable for model selection. While model selection criteria have been previously used in geotechnical engineering (Honjo et al., 1994; Li et al., 2013; Tang et al., 2013a, 2013b; Cao and Wang, 2013; Wang et al., 2013; Cao and Wang, 2014a, 2014b; Wang et al., 2014), we hope that the application presented herein can help, in conjunction with other recent contributions (Wang and Aladejare, 2015), to promote a wider use of model selection in rock engineering and, in

20

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27

Table 1 A summary of models to estimate the deformation moduli of rock masses based on geotechnical classification indices. Proposed by

Empirical correlations

Input parameters

Types

Limitations

Bieniawski (1978) Serafim and Pereira (1983) Kim (1993) Read et al. (1999) Jašarević and Kovačević (1996) Aydan et al. (1997) Verman et al. (1997) Diederichs and Kaiser (1999)

Erm = 2RMR − 100 Erm = 10((RMR − 10) / 40) Erm = 300 × 10−3exp(0.07RMR) Erm = 0.1(RMR / 10)3 Erm = exp.(4.407 + 0.081 ∗ RMR) Erm = 0.0097RMR3.54 Erm = 0.4Hα10((RMR-20)/38 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Erm = 7(±3) 10ðRMR44Þ=21 Erm = 0.0876RMR Erm = 0.0876RMR + 1.056(RMR − 50) + 0.015 (RMR − 50)2 Erm = 100exp[−((RMR − 100) / 37)2] Erm = 9 × 10−7RMR3.868 Erm = Ei(0.0028RMR2 + 0.9exp(RMR / 22.82)) Erm = Ei[0.5(1 − cos(πRMR / 100))] Erm = Ei[RMR / (RMR + β(100 − RMR))] Erm = Ei exp.((RMR − 100) / 36) Erm = Ei10[(RMR − 100)(100 − RMR) / 4000 ∗ exp.(−RMR / 100)] Erm = Ei exp.[−((RMR − 116) / 41)2] qffiffiffiffiffiffi σ ci 10((GSI − 10) / 40) Erm = 100 qffiffiffiffiffiffi σ ci 10((GSI − 10) / 40) Erm = (1 − D / 2) 100

RMR RMR RMR RMR RMR RMR RMR, H, α RMR

Linear Exponential Exponential Power Exponential Power Exponential Non-linear

RMR N 50

RMR RMR RMR RMR RMR, Ei RMR, Ei RMR, Ei, β RMR, Ei RMR, Ei RMR, Ei

Linear Polynomial Gaussian function Power Non-linear Trigonometric Fractional Exponential Exponential Gaussian function

RMR N 50 RMR ≤ 50

GSI, σci

Exponential

GSI, σci, D

Exponential

σci b 100

GSI, D

Exponential

σci N 100

GSI, Ei, D

Sigmoid equation

GSI, D

Sigmoid equation

GSI, Ei, D

Non-linear

Q Q Q Q Q, σci RQD, Ei RQD, Ei, WD RQD, Ei, σci, WD RQD, Ei RMi

Logarithmic Power Power Exponential Power Linear Non-linear Non-linear Exponential Power

Galera et al. (2005) Galera et al. (2005) Shen et al. (2012) Khabbazi et al. (2013) Nicholson and Bieniawski (1990) Mitri et al. (1994) Aydan and Kawamoto (2000) Galera et al. (2005) Sonmez et al. (2006) Shen et al. (2012) Hoek and Brown (1997) Hoek et al. (2002) Hoek et al. (2002) Hoek and Diederichs (2006) Hoek and Diederichs (2006)

Erm = (1 – D/2) 10((GSI-10)/40)   Erm = Ei 0:02 þ 1D=2 60þ15DGSI 11 1þe   1D=2 Erm = 100 75þ25DGSI 1þe

Sonmez et al. (2004) Barton (1983) Barton (1995) Palmström and Singh (2001) Kang et al. (2012) Barton (1996) Coon and Merritt (1970) Kayabasi et al. (2003) Gokceoglu et al. (2003) Zhang and Einstein (2004) Palmström and Singh (2001)

11

a 0.4

Erm = Ei (s ) ; s = exp.[(GSI − 100) / (9 − 3D)], a = 0.5 + 1 / 6(exp(−GSI / 15) + exp.(−20 / 3)) Erm = 25log10(Q) Erm = 10Q1/3 Erm = 8Q0.4 Erm = 10(0.32logQ + 0.585) Erm = 10Q1/3 c; Qc = Q(σci / 100) Erm = Ei(0.0231RQD − 1.32) Erm = 0.135[Ei (1 + RQD / 100) / WD]1.1811 Erm = 0.001[(Ei/σci) (1 + RQD / 100) / WD]1.5528 Erm = Ei 10(0.0186RQD – 1.91) Erm = 7RMi0.4

QN1 1 b Q b 30

Erm deformation modulus of rock mass (GPa), H depth of tunnel (m), α a parameter related to RMR, σci uniaxial compressive strength of intact rock (MPa), Ei elastic modulus of intact rock (GPa), WD weathering degree, β a constant to be determined using a minimization procedure for experimental values, RMi rock mass index.

particular, in application related to the estimation of the deformation modulus of rock masses. Model selection criteria can be used to identify the “best” model, although predictions conducted with such “best” models are sometimes associated to large uncertainties due to the high variability of data employed to build them. Within this context, further improvements could be achieved with a Bayesian updating approach, so that uncertainties can be reduced when new project-specific observations become available, hence helping designers to “maximize the value of new project-specific observations” (Zhang et al., 2004). This is because the Bayesian updating is expected to reduce the predictive uncertainty (i.e., the COV of the predictive distribution), which may be of interest to reliability-based designs, as the computed probability of failure of a design will be influenced by such predictive uncertainty. In this study, we present an approach to select the most appropriate model, among four commonly used candidate models, to estimate the deformation modulus of a rock mass, Erm, based on its RMR or GSI values. The proposed approach builds on the use of “information criteria” to rank models, considering both their fit and complexity. In addition, we propose a procedure, within the Bayesian framework for model updating, to systematically update the predictive distribution of Erm, and its associated predictive uncertainty, when new “project-specific” data are available. Finally, we use an example case of an assumed circular rock tunnel design to show that the reduction of predictive uncertainties achieved with the Bayesian updating has a significant influence on the computed reliability estimates of a rock engineering design.

2. Bayesian models: inference, selection and prediction 2.1. Introduction The relationship between observed (input or independent) and estimated (output or dependent) variables in real rock engineering problems, such as estimating the deformation modulus of rock masses, is often complex. But, in some cases, we may build a model, zi =Mi(x,Θ), to describe such complex behavior, where zi (i = 1 , 2 , . . . , r) denote the (observable) dependent variables, x = (x 1 , x 2 , ... , x m ) is the m-dimensional vector of the (observable) independent variables, and Θ = (θ1, θ2, ... , θp) is the p-dimensional vector of (unobservable) model parameters that need to be estimated based on observed data. (Geyskens et al., 1998). Mi(⋅, ⋅) represents the relationship between them. For instance, the correlation proposed by Serafim and Pereira (1983) to estimate the deformation modulus of a rock mass, i.e., Erm = 10(a⁎RMR+b), is a model with r = 1, in which the input variable x= RMR is one dimensional (m = 1), and in which the vector of parameters is given by Θ =(a = 1/40,b = ‐ 1/4) (Note, therefore, that p = 2). That is, in this model, RMR is the (observable) independent variable, Erm is the (observable) dependent variable, and Θ = (a, b) are the (unobservable) model parameters to be estimated. Since a model is only an approximation to the complex real-world phenomenon, models are usually imperfect and their predictions are uncertain (Geyskens et al., 1993). In addition, simplifications are unavoidable and, as it is clear from Table 1, there may be many different models for the same task. Here, we use a Bayesian approach to infer and update a set of model parameters, Θ, for

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27

a specific rock engineering problem (see Sections 2.2 and 2.3), to select the best model out of a set of candidate models (see Sections 2.4 and 2.5), and to make predictions based on such models (see Section 2.6). 2.2. Bayesian inference of model parameters Bayes' theorem can be used to ‘update’ our (uncertain) state of knowledge about the model parameters, Θ, after a set of observations of dependent and independent variables, y = (zi, x), become available. In the context of inference of parameters, Bayes' theorem can be expressed as (Ang and Tang, 1975; Lunn et al., 2012): f ðΘjyÞ ¼

LðΘ; yÞf ðΘÞ f ðyÞ

ð1Þ

where f(Θ) is the “prior” probability density function (PDF), which is related to the initial knowledge about Θ before the newly observed data is taken into account; L(Θ; y) ≡ f(y| Θ) is the “likelihood” function, which is the conditional probability of observing the additionally available data, y, given specific values of Θ; and f(Θ| y) is the “posterior” PDF that incorporates both the subjective judgment and the observed data. The normalizing constant f ðyÞ ¼ ∫LðΘ; yÞf ðΘÞdΘ ensures that f(Θ| y) integrates to one. 2.3. Posterior predictive distribution Sometimes, we would like to make predictions after having observed data, y. For example, in this study we aim to predict the deformation modulus of a rock mass, so that it can be used in reliability-based design. Assuming that observations are conditionally independent given Θ, we derive the posterior predictive distribution (Ntzoufras, 2009; Lunn et al., 2012) ~ jyÞ ¼ f ðy

Z Θ

~jΘÞf ðΘjyÞdΘ f ðy

ð2Þ

~ denotes our predicted quantity of interest. where y In general, analytical solutions to Eq. (2) are not available due to the ~ jΘÞ or f(Θ| y), but such integrals can be computed nucomplexity of f ðy merically. In this work, for instance, we employ WinBUGS (Lunn et al., 2000) to compute the predictive distributions given by Eq. (2). 2.4. The need for model selection criteria

^ = ln f ðyjΘÞ ^ is the maximum log-likelihood of the model, Θ ^ where lðΘÞ are the maximum likelihood estimates for model parameters an k is the number of the estimated parameters. The BIC is closely related to AIC, and it is defined as   ^ þ k  ln ðnÞ BIC ¼ −2l Θ

ð4Þ

where n is the number of observations employed in the likelihood. In general, BIC has a stronger penalty term (Gardoni et al., 2009). It is clear that models with lower AIC (or BIC) are preferred. The DIC can be regarded as a generalization of AIC (Lunn et al., 2012), with the advantage that it is particularly suitable within Bayesian frameworks. Following Spiegelhalter et al. (2002), for data y and model parameters Θ, a deviance is defined as −2 times the log-likelihood, i.e., DðΘÞ ¼ −2lðΘÞ

ð5Þ

and the DIC is defined as   DIC ¼ D þ pD ¼ D Θ þ 2pD

ð6Þ

whereD is the posterior mean deviance (i.e., the average of D(Θ)) that denotes the measure of fit, DðΘÞ is the deviance of posterior means (i.e., deviance evaluated at the average of Θ) and pD is effective number of parameters that represents a measure of complexity, i.e., pD = D−DðΘÞ. (For non-hierarchical models with little prior information, the effective number of parameters in DIC, pD, should be approximately equal to the true number of parameters in AIC or BIC, k (Lunn et al., 2012).) Therefore, DIC can also be considered as a Bayesian measure of model fit and complexity (Spiegelhalter et al., 2002). Although DIC is mathematically more complex than AIC or BIC, and in this case it could be justified to use AIC or BIC since they give similar results, we have mainly used DIC in this study because it has been reported as a more natural measure within Bayesian framework (Lunn et al., 2012). (Note also that “in simple models, DIC is asymptotically equal to the AIC, but is more general and can also be applied in hierarchical models, where the number of parameters is note clearly defined” (Luoma, 2008) and that there are available software, such as WinBUGS, that include functions to directly compute DIC values for different models.) 2.5. Model selection

As shown in Table 1, many different types of models, —such as linear, power, exponential etc.—, can be used to predict Erm. Therefore, an objective methodology is needed to select the “best” model, given some observations, among a set of candidate models. Simply comparing the values of R2 without penalizing model complexity is a poor model selection technique, as it is possible to improve the goodness of fit of a model by increasing its number of parameters. But too many parameters may result in overfitting and, therefore, too complex models should be penalized. Model selection techniques based on information criteria provide a convenient way to rank models considering both model fit and model complexity. Akaike information criterion (AIC) (Akaike, 1973) and Bayesian information criterion (BIC) (Schwarz, 1978) were two initial proposals in this direction; more recently, the deviance information criterion (DIC) (Spiegelhalter et al., 2002) was proposed with Bayesian applications in mind. Therefore, we will conduct Bayesian model selection based mainly on DIC, although results computed with the widelyused AIC and BIC (which are quite similar in this case) will also be provided for comparison. The AIC is defined as   ^ þ 2k AIC ¼ −2l Θ

21

ð3Þ

Many models to estimate Erm are possible (See Table 1). Due to their more significant influence on the deformability of the rock mass (Khabbazi et al., 2013), here we focus on models that use RMR (or GSI) to predict Erm. In particular, we consider four types of candidate models: liner, power, exponential, and logistic (see Table 2). The immediate purpose of model selection is to identify the most appropriate model for prediction; often, it also aims to “estimate model parameters that are of particular interest” (Johnson and Omland, 2004). Therefore, in addition to the DIC values employed for model selection, the estimates of model parameters will also be presented in this section. Table 2 The functional form of four commonly used models, with a, b and c being the corresponding model parameters. Model type

Reference

Formula

Linear Power Exponential Logistic (“S” shape)

Bieniawski (1978) Aydan et al. (1997) Serafim and Pereira (1983) Hoek and Diederichs (2006)

Erm = a∙RMR⁎ Erm = a∙RMRb Erm = a∙exp.(b∙RMR) Erm = a/[1 + b∙exp.(c∙RMR)]

⁎ The intercept of the linear model is set to zero to ensure that Erm equals zero at RMR = 0.

22

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27

As shown by Haldar and Babu (2008) and Miranda et al. (2009), the lognormal distribution is often suitable for modeling non-negative geotechnical parameters for which the assumption of normality is questionable. Therefore, we will assume that the Erm is lognormally distributed in this study. To illustrate the Bayesian model selection approach, we use a subset of the database of reliable in situ Erm measurements reported by Hoek and Diederichs (2006). (Because we were using data scanned from their original figure, it was only possible to identify 312 observations out of the total of 494 that Hoek and Diederichs (2006) reported to be included in their analysis.) It has been suggested that “if the 1989 version of Bieniawski's RMR classification is used, then GSI = RMR89 – 5” (Hoek and Brown, 1997) and if the 1976 version of RMR is used, then GSI = RMR76 (Hoek and Diederichs, 2006). Therefore, for simplicity, it will be assumed that RMR = GSI for the analyses below. Given the Hoek and Diederichs (2006) data described above, the four models listed in Table 2 are fitted, in a Bayesian model assessment framework, using WinBUGS. Then the DIC of each model, which can be easily computed in Markov chain Monte Carlo (MCMC) analysis (Spiegelhalter et al., 2002), is employed for model selection, as “the minimum DIC is intended to identify the model that would be expected to make the best short-term predictions” (Lunn et al., 2012).

2.5.2. Power model Similarly, we can linearize the power model (see Table 2) by considering the logarithm for the i-th observation, i.e.,

2.5.1. Linear model For the i-th observation and taking its logarithm, the linear model (see Table 2) can be expressed as:

with μρ = μb = 0, τρ = τb =1 × 10−6 and η = λ = 1× 10−3. The summary statistics for model parameters and the DIC computed with WinBUGS and MCMC for the power model are listed in Table 3.

Di ¼ ρ þ F i þ εi

i ¼ 1; 2; :::; n

ð7Þ

where Di = ln (Erm,i), ρ = ln (a), Fi = ln (RMRi), n is the number of observations and εi is the Gaussian random correction accounting for model error. It is assumed that εi is normally distributed with common mean μ = 0 and variance σ 2, i.e., εi ~ N(0, σ 2). Based on the above assumptions, it is easy to show that, for given ρ, Di is normally distributed, i.e., Di ~ N(ρ + Fi, σ 2). That is, the likelihood function for parameters ρ and τ (with τ = 1/σ 2) can be expressed as: (rffiffiffiffiffiffi ) h τ i τ  exp − ðDi −ρ− F i Þ2 2π 2 i¼1 n

Lðρ; τÞ ¼ ∏

ð8Þ

and when no detailed information is available on model parameter, ρ, it is common to use a normal distribution with a large variance, indicating a very flat shape with very similar values over a wide range of ρ (Ntzoufras, 2009). The basic idea, therefore, is to use a prior that carries little information on model parameter, ρ, and that is why we call it the “low-informative” prior. It is considered as f ðρÞ ¼

rffiffiffiffiffiffi  2  τρ τρ  exp − ρ−μ ρ 2π 2

ð9Þ

with μρ = 0, τρ = 1 × 10−6. Similarly, it is widely accepted that a gamma distribution can be applied as a prior on τ (Ntzoufras, 2009). The reason is that with small scale and rate parameters, the variance of the gamma distribution is very large (e.g. the mean = 1 and variance = 1000 for η = λ = 1 × 10− 3), therefore, it can also be considered as a low-informative prior, as shown in Eq. (10). f ðτ Þ ¼

ηλ τ λ−1 expð−ητÞ Γ ðλÞ

for

τ N0

ð10Þ

with η = λ = 1 × 10−3. The inference of model parameters is conducted in WinBUGS using MCMC methods. The summary statistics for all model parameters are listed in Table 3, and the DIC for the assumed linear model is also listed in Table 3.

Di ¼ ρ þ bF i þ εi

ð11Þ

i ¼ 1; 2; :::; n

where Di is also normally distributed, i.e., Di ~ N(b⋅ Fi +ρ, σ 2). Therefore, the likelihood function for parameters ρ, b and τ can be expressed as: (rffiffiffiffiffiffi ) h τ i τ 2  exp − ðDi −ρ−b  F i Þ 2π 2 i¼1 n

Lðρ; b; τ Þ ¼ ∏

ð12Þ

and low-informative priors on ρ, b and τ are considered as: f ðρÞ ¼

rffiffiffiffiffiffi  2  τρ τρ  exp − ρ−μ ρ 2π 2

ð13Þ

f ðbÞ ¼

rffiffiffiffiffiffi h τ i τb 2 exp − b ðb−μ b Þ 2π 2

ð14Þ

f ðτ Þ ¼

ηλ τ λ−1 expð−ητÞ Γ ðλÞ

for

τ N0

ð15Þ

2.5.3. Exponential model For the exponential model (see Table 2), we also consider its logarithm for the i-th observation, i.e., Di ¼ ρ þ bRMRi þ εi

i ¼ 1; 2; :::; n

ð16Þ

As before, it can be seen that Di is normally distributed, i.e., Di ~ N(b ⋅ RMRi +ρ,σ 2). Therefore, the likelihood function for parameters ρ, b and τ can be expressed as: (rffiffiffiffiffiffi ) h τ i τ 2  exp − ðDi −ρ−b  RMRi Þ 2π 2 i¼1 n

Lðρ; b; τ Þ ¼ ∏

ð17Þ

The prior functions for the exponential model are the same as those employed for the power model (see Eqs. (13)–(15)) and not reported herein. The summary statistics for all model parameters and the DIC computed with WinBUGS and MCMC for the exponential model are presented in Table 3. 2.5.4. Logistic model The logistic function is a common sigmoid function, similar to the one adopted by Hoek and Diederichs (2006) to “constrain the increase of modulus as the rock becomes more massive”. For the i-th observation, the logarithm of the logistic model (see Table 2) can be expressed as  Di ¼ ln

 a þ εi 1 þ b  expðc  RMRi Þ

ð18Þ

where, since Erm is assumed to be lognormally distributed, Di is normally  h i  a distributed, i.e., Di  N ln 1þb expðcRMR ; σ 2 . That is, the likelihood iÞ function can be written as: Lðφ1 ; φ( 2 ; φ3 ; τ C Þ " rffiffiffiffiffiffi   2 #) n τC τC a  exp − Di − ln ¼ ∏ 2π 2 1 þ b  expðc  RMRi Þ i¼1 with φ1 = ln (a), φ2 = ln (b+ 1) and φ3 = ln (−c).

ð19Þ

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27

23

Table 3 Summary statistics, including prediction uncertainties, for the parameters of the different posterior models computed with WinBUGS using MCMC. Percentiles Model

Parameter

Mean

Std.dev

MC error

2.5%

50%

97.5%

Linear

a ρ = ln(a) σ τ D DðΘÞ pD DIC a b σ τ θ = ln(a) D DðΘÞ pD DIC a b σ τ A = ln(a) D DðΘÞ pD DIC a b c σ τ φ1 = ln(a) φ2 = ln(b + 1) φ3 = ln(−c) D DðΘÞ pD DIC

0.1349 −2.005 1.089 0.8473 938.1 936.1 2.002 940.1 7.037E−06 3.459 0.8234 1.482 −12.09 763.9 760.9 3.021 766.9 0.1149 0.06839 0.8077 1.54 −2.183 751.3 748.3 3.026 754.4 54.4 957.4 −0.08408 0.7977 1.579 3.962 6.829 −2.478 743.4 739.6 3.734 747.1

0.008277 0.06126 0.04409 0.06836

7.937E−05 5.878E−04 4.344E−04 6.819E−04

0.1196 −2.124 1.006 0.7186

0.1346 −2.006 1.088 0.8453

0.1521 −1.883 1.18 0.9878

5.354E−06 0.1615 0.03347 0.12 0.6641

5.097E−08 1.404E−03 3.350E−04 1.195E−03 0.005689

1.540E−06 3.138 0.7618 1.255 −13.38

5.633E−06 3.459 0.8224 1.478 −12.09

2.102E−05 3.775 0.8927 1.723 −10.77

0.02295 0.003084 0.03291 0.1249 0.1982

2.372E−04 3.220E−05 3.081E−04 0.001164 0.002034

0.0763 0.06245 0.7468 1.304 −2.573

0.1129 0.06837 0.8066 1.537 −2.181

0.1658 0.07445 0.8757 1.794 −1.797

14.46 263.3 0.005359 0.03216 0.1269 0.2601 0.271 0.06386

1.27 21.39 4.59E−04 3.09E−04 0.001204 0.02273 0.02201 0.005495

33.24 534.9 −0.09395 0.7372 1.339 3.504 6.284 −2.601

52.4 918.7 −0.08381 0.7965 1.576 3.959 6.824 −2.479

88.18 1574 −0.07417 0.8641 1.84 4.479 7.362 −2.365

Power

Exponential

Logistic

Since little prior information on φj is known, flat priors on φj ( j = 1 , 2, 3) and τ will be used in this study. They are given as:    rffiffiffiffiffiffi 2  τj τj  f φj ¼ exp − φ j −μ j 2π 2 f ðτ Þ ¼

ηλ τ λ−1 expð−ητÞ Γ ðλÞ

for

τ N0

ð20Þ

Erm. The exponential model and the logistic model predict similar values when RMR is less than about 75–80, but differences become more substantial for higher RMR values. (Remember, for instance, that the logistic model constrains the increase of modulus as the rock mass becomes more massive.)

ð21Þ

The assumed parameters are: η = λ = 1 × 10−3, μ j = 0 and τj = 1 × 10−6 ( j = 1 , 2, 3). The summary statistics of model parameters computed for the logistic model, as well as its DIC computed with WinBUGS and MCMC are listed in Table 3. Table 3 shows that the logistic model results in the smallest DIC among the four candidate models, and the DIC differences between the logistic model and the other three models in DIC indicate that, with an assumed lognormal distribution of Erm, the logistic model would be the most appropriate model for Hoek and Diederichs (2006) data, with other models being significantly less adequate. (According to Lunn et al. (2012), “differences of more than 10 might definitely rule out the model with the higher DIC, and differences between 5 and 10 are substantial”.) Using the mean values of the parameters listed in Table 3, Erm-RMR relationships for the four candidate models are plotted in Fig. 1 together with Hoek and Diederichs (2006) data. Results suggest that the linear model is not able to capture the complexity of the data. They also show that the exponential model predicts the largest Erm at RMR = 100, while the other three models give significantly smaller values of

Fig. 1. In situ observations of rock mass deformation moduli reported by Hoek and Diederichs (2006) and predictions with the four candidate models using the mean values of their parameters.

24

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27

Table 4 The MLE of model parameters as well as their corresponding AIC and BIC values. Model type (see Table 1) Parameter

Linear

Power

Exponential

Logistic

a b c σ AIC BIC Rank

0.13463 / / 1.0846 2247.4 2254.9 4

5.53E−06 3.4634 / 0.81910 2074.2 2085.4 3

0.11268 0.068401 / 0.80274 2061.6 2072.9 2

53.110 900.98 −0.083253 0.79162 2054.8 2069.9 1

Note: The AIC and BIC values are computed from the maximum log-likelihood using Eq. (17). The models are ranked according to their AIC and BIC values.

Next, the model parameters are also estimated using the maximum likelihood estimation (MLE) method. AICs and BICs for the four models will be reported for comparison. As before, we assume that Erm is log-normally distributed. Thus, ln(Erm) is normally distributed, i.e., ln(Erm) = ln [f(RMR)] + ε with ε ~ N(0, σ 2 ); alternatively, Erm can be expressed as: Erm = f(RMR) ⋅ eε with ε ~ N(0, σ 2). The likelihood function for the lognormally-distributed Erm model can be expressed as (Wang and Liu, 2006): "

2 #

ln Erm;i − ln ðμ i Þ 1 pffiffiffiffiffiffi exp − 2σ 2 i¼1 Erm;i 2π σ n

LðyjΘÞ ¼ ∏

ð22Þ

where L(y| Θ) is the likelihood of the observed data, y, given the parameter vector Θ = (a, b, c, σ) (see Table 2), Erm,i is the i-th observed value of Erm, μ i is the i-th estimated value of the mean of Erm given RMRi, computed using the formulas listed in Table 2 for given model parameters Θ. The maximum likelihood estimates of the parameters for our subset of Hoek and Diederichs (2006) dataset are obtained using the “interiorpoint” method in the optimization toolbox of MATLAB (Myung, 2003; Wang and Liu, 2006). The values of the estimated model parameters, as well as the AIC and BIC values computed for the four candidate models, are listed in Table 4. Table 4 shows that the model ranking based on BIC is identical to that based on AIC. It also coincides with the ranking provides by DIC in Table 3. Therefore, we can safely assume that the logistic model is the best for this dataset. (It should be noted that AIC values here are quite different from DIC values shown in Table 3 because of the log transformation performed; otherwise they should be very similar as the priors are low-informative.) In addition, it should be noted that models fitted assuming lognormal distribution for Erm are better, as they have lower AICs (or BICs) Table 5 AICs, BICs and DICs for Serafim and Pereira (1983) data. Method MLE Bayesian model selection

AIC BIC DIC

Linear

Power

Exponential

Logistic

92.39 93.67 33.63

82.80 84.72 23.86

81.00 82.92 22.62

83.03 85.58 23.90

than those fitted assuming normal distributions (for the sake of space, the results are not reported herein). For the observations reported by Serafim and Pereira (1983), we also make the model selection using a similar approach and the values of AICs, BICs and DICs listed in Table 5. We can see that the exponential model is the most appropriate model for Serafim and Pereira (1983) data, which provides supporting evidence for the analysis conducted by Serafim and Pereira (1983). Also, it should be noted that the exponential model, instead of the most complex model, i.e., the logistic model, is selected as the most appropriate model, hence illustrating how model selection criteria can help avoid overfitting problems. 2.6. Predictive analyses In practical engineering, one wants to use the selected model to predict the quantity of interest, which in this case is Erm. We know that, for Hoek and Diederichs (2006) data, the logistic model appears as the most suitable one based on the DIC (AIC and BIC) results of Section 2.5, so that, given this specific model, the Bayesian framework can be employed to obtain probabilistic estimates of its parameters. They can then be employed to conduct predictive analyses in the context of uncertainty, and they are also useful for reliability-based design. As an example, the predictive distribution of Erm at a given value of RMR is of particular interest in this study. For instance, assuming a value of RMR = 60, we can obtain the predictive distribution of D = ln(Erm) (see Table 6). Based on the distribution of D, obtaining the distribution of Erm = exp(D) is immediate. (Remember that D is normally distributed; therefore, the predictive distribution of Erm is directly obtained by a transformation that is available in standard engineering probability texts; see, e.g., Ang and Tang (1975).) 3. Bayesian updating with new data 3.1. Motivation In engineering, one may want to derive a “general” model to predict the (distribution of) an engineering property of interest, given information about such property observed from various locations or projects. However, the variability of data sources usually implies that the predictive uncertainty associated with this “general” model (as indicated by its COV) is often large, hence reducing its usefulness in practice for a particular project at a specific site. An alternative, of course, would be to derive a project- or site-specific model using exclusively new data gathered from that project (or location) as the project proceeds. However, this would imply to completely neglect the information provided by the original (general) dataset, which might be still available, particularly when the new project-specific dataset is limited. As presented in Section 3.2 (see also Feng and Jimenez (2014)), the Bayesian updating framework provides a way to systematically and effectively combine the information from the “general” and the “project-specific” data, producing an “updated” predictive distribution with reduced uncertainty. As an example, we show how the fitted “general” logistic model (see Table 3), that was identified as more feasible based on the Bayesian model selection conducted using Hoek and Diederichs (2006) dataset,

Table 6 Comparison of predicted medians, means, standard deviations and coefficient of variations (COV) of D and Erm at RMR = 60.

D

Erm

Based on Hoek and Diederichs (2006) data Based on Serafim and Pereira (1983) data Bayesian updating Based on Hoek and Diederichs (2006) data Based on Serafim and Pereira (1983) data Bayesian updating

Median (MPa)

Mean (MPa)

Std.dev (MPa)

COV

2.01 2.75 2.53 7.49 15.58 12.51

2.01 2.75 2.53 10.36 18.22 13.82

0.81 0.54 0.45 10.06 11.49 6.50

0.40 0.20 0.18 0.97 0.63 0.47

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27

25

can be updated using new available data. (As “updating” dataset, we employ the dataset of Erm - RMR values employed by Serafim and Pereira (1983) to develop their well-known correlation.) In particular, we show how the predictive uncertainty about Erm for a given RMR value (such as the relatively large value of COV[Erm] = 0.97 reported in Table 6 for RMR = 60) can be reduced using Bayesian updating when a set of new observations from a specific “project” are available.

Serafim and Pereira (1983) dataset. One can see that such updated distribution has the smallest uncertainty as indicated by a “narrower” probability distribution function (PDF) and a minimum COV of Erm. This illustrates that Bayesian updating can effectively reduce the predictive uncertainty. Table 6 also shows that the Bayesian updating method predicts the smallest standard deviation of Erm at RMR = 60. (These conclusions are also valid for other RMR values, although results are not reported herein.)

3.2. Bayesian updating

4. Influence on reliability estimates

To illustrate how to conduct the Bayesian updating, this section presents an example in which, for simplicity, we will assume that the new “project-specific” data that become available are those from Serafim and Pereira (1983). To that end, let us compute the median, mean and standard deviation of D = ln(Erm) employing (i) the dataset by Hoek and Diederichs (2006) (assuming the sigmoid model) and (ii) the dataset by Serafim and Pereira (1983). (μ D1 and σD1 are used to denote, respectively, the mean and standard deviation of D obtained for the sigmoid model and Hoek and Diederichs (2006) dataset; similarly, μ D2 and σD2 denote the parameters that would have been obtained applying the same model selection procedure directly with Serafim and Pereira (1983) dataset. Example results of these values, as computed for RMR = 60, are listed in Table 6.) If we assume that D is normally distributed with a known variance, the updated mean, μD″, based on the combined information of both datasets can be computed, using the Bayesian approach, as (Ang and Tang, 1975; Zhang et al., 2004):

To illustrate the effect of the different predictive distributions on a reliability-based design, we consider an example case with a circular tunnel of radius R and subjected to a hydrostatic stress p0. For simplicity, we assume that the rock mass remains elastic and that no support is applied. The inward displacement of the tunnel wall, u, can be computed as (Hoek and Brown, 1980)

μ D″ ¼

μ D1 σ 2D2 þ μ D2 σ 2D1 σ 2D2 þ σ 2D1

ð23Þ

and the updated standard deviation, σD″, is computed as σ D″ ¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 2 u σD σD 2 1 t σ 2D2 þ σ 2D1

ð24Þ

Table 6 also lists the updated mean and standard deviation of D, for RMR = 60, as computed using Eqs. (23) and (24). Once the information about the parameters of the distribution of D is available, the updated mean and standard deviation of Erm can be obtained using a standard transformation from normal to lognormal distribution (Ang and Tang (1975); see Table 6). Fig. 2 compares the predictive distributions of D = ln(Erm), for RMR = 60, obtained based on Hoek and Diederichs (2006) dataset only, on Serafim and Pereira (1983) dataset only, and on the Bayesian updating of the original model based on Hoek and Diederichs (2006) dataset using the new information provided by the

u 1þν p ¼ R Erm 0

ð25Þ

where ν is the Poisson's ratio of the rock mass (herein assumed as ν = 0.2). Selecting a limiting convergence ratio of u/R = 1%, the limit state function (LSF) becomes g ðxÞ ¼

0:01 0:01Erm −1 −1 ¼ ð1 þ ν Þp0 u=R

ð26Þ

The reliability index, β, and the probability of failure, Pf, of such LSF for the lognormal distributions of Erm listed in Table 6 are computed using FORM (First Order Reliability Method) (see, e.g., JimenezRodriguez et al. (2006)). The results are shown in Table 7. Taking p0 = 12.5 MPa and RMR = 60 as an example, the probability of failure computed after the Bayesian updating based on the new dataset (Pf = 1.0480E −6) is much smaller than that based on Hoek and Diederichs (2006) data (Pf = 0.0248). At first sight, it could be argued that the reduction of the computed probability of failure after updating with Serafim and Pereira (1983) dataset is due to the “shift towards the right” of the Erm distribution due to the good quality of the rock mass at the sites studied by Serafim and Pereira (1983) (Note that these sites had been mainly considered as locations for dams). However, Table 7 shows that, in this case, the reduction of the predictive uncertainty of Erm (as indicated by its smaller COV) has an even stronger influence on the computed estimates of Pf and β. (Note, for instance, that the Erm model based only on Serafim and Pereira (1983) dataset has higher mean of Erm, but a higher Pf; whereas the Erm model based on Hoek and Diederichs (2006) dataset updated with Serafim and Pereira (1983) dataset has a lower Erm mean but a lower Pf.) As indicated by previous research (see Jimenez and Sitar (2009)), this observation illustrates the importance of the tails of the distributions of quantities of interest in a reliability analysis, which are of course controlled by their COV.

Table 7 Reliability index and probability of failure for different distributions of Erm for an assumed RMR = 60, and for different p0 values. p0(MPa)

Fig. 2. Comparison of predictive distributions of D = ln(Erm) for RMR = 60 based on Hoek and Diederichs (2006) data, Serafim and Pereira (1983) data and Bayesian updating.

25.0

12.5

6.25

Erm based on Hoek and Diederichs (2006) data: β 1.1132 1.9637 0.1328 0.0248 Pf

2.8142 0.0024

Erm based on Serafim and Pereira (1983) data: β 2.8282 4.0261 0.0023 2.8355E−5 Pf

5.2240 8.7553E−8

Erm based on combined information using Bayesian updating approach: β 3.1934 4.7439 6.2945 7.0301E−4 1.0480E−6 1.5423E−10 Pf

26

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27

And, since the predictive COV of Erm decreases significantly after the Bayesian updating in our example case, this justifies why model updating has such a significant influence on our computed estimates of the probability of failure. This is, of course, of significant interest to reliability-based design, in which estimates are employed to accept design, to choose among designs, or to compute risks. 5. Conclusions This paper proposes an approach to select the most appropriate model, among four commonly used candidate models —linear, exponential, power, and logistic— to estimate the deformation modulus of a rock mass, Erm, based on its RMR or GSI values. The proposed approach, based on “information criteria” (such as DIC, AIC and BIC), considers model fit while penalizing model complexity. Results obtained with the proposed model selection approach suggest that the logistic model is the most adequate model (among those considered) when the extensive dataset of Erm-GSI values published by Hoek and Diederichs (2006) is employed to fit and rank models. In addition, we propose a procedure to systematically update the predictive model, and its uncertainty, as new project-specific data become available. The advantage of such Bayesian updating scheme is that it often reduces the generally large predictive uncertainty that is commonly associated to “general” models fitted using data from a wide variety of rock mass types and of project locations. In other words, the predictive Erm distribution obtained based on the general dataset can be updated, and its predictive uncertainty reduced, when “project-specific” data are employed. Such updating of the predictive uncertainty associated to a model can have significant consequences in the context of a reliability-based design, in which the probabilities of failure are employed to decide about the feasibility of a design, to choose among design alternatives, or to compute risks. In particular, we use an example case to show that the reduction of predictive uncertainties about Erm achieved with the Bayesian updating associated to “project-specific” data can significantly affect the computed estimates of the probability of failure (or, equivalently, of its reliability index) of an assumed circular rock tunnel design. Acknowledgments The first author is supported by China Scholarship Council (CSC) and, for insurance coverage, by Fundación José Entrecanales Ibarra. Jiayi Shen from Zhejiang University kindly provided his database of Erm v.s. RMR. Their supports are greatly appreciated. References Akaike, H., 1973. Information theory and an extension of the maximum likelihood principle. In: Pertaran, B.N., Csaaki, F. (Eds.), The Second International Symposium on Information Theory, Budapest, Hugary. Acadeemiai Kiadi, pp. 267–281. Ang, A.S., Tang, W., 1975. Probability concepts in engineering planning and design. Volume I: Basic principles. John Wiley & Sons, Inc., New York. Aydan, Ö., Kawamoto, T., 2000. Assessing mechanical properties of rock masses through RMR rock classification system. The GeoEng 2000 symposium, Melbourne, Australia (on CD). Aydan, Ö., Ulusay, R., Kawamoto, T., 1997. Assessment of rock mass strength for underground excavations. The 36th US rock mechanics symposium, New York, pp. 777–786. Barton, N., 1983. Application of Q system and index tests to estimate shear strength and deformability of rock masses. International Symposium on Engineering Geology and Underground Construction, Lisbon, Portugal, pp. 51–70. Barton, N., 1995. The influence of joint properties in modelling jointed rock masses. The 8th International Congress on Rock Mechanics, Tokyo, Japan, pp. 1023–1032. Barton, N., 1996. Estimating rock mass deformation modulus for excavation disturbed zone studies. International Conference on Deep Geological Disposal of Radioactive Waste, Winnepeg. Canadian Nuclear Society, pp. 133–144. Bieniawski, Z.T., 1978. Determining rock mass deformability: experience from case histories. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 15 (5), 237–247. http://dx.doi.org/10. 1016/0148-9062(78)90956-7. Burnham, K.P., Anderson, D.R., 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. Springer-Verlag, New York.

Cao, Z., Wang, Y., 2013. J. Geotech. Geoenviron. 139 (2), 267–276. http://dx.doi.org/10. 1061/(ASCE)GT.1943–5606.0000765. Cao, Z., Wang, Y., 2014a. Bayesian model comparison and characterization of undrained shear strength. J. Geotech. Geoenviron. 140 (6), 04014018. http://dx.doi.org/10. 1061/(ASCE)GT.1943-5606.0001108. Cao, Z., Wang, Y., 2014b. Bayesian model comparison and selection of spatial correlation functions for soil parameters. Struct. Saf. 49, 10–17. http://dx.doi.org/10.1016/j. strusafe.2013.06.003. Coon, R.F., Merritt, A.H., 1970. Predicting in situ modulus of deformation using rock quality indices. Determination of the in-situ Modulus of Deformation of Rock, ASTM STP 477, Philadelphia, pp. 154–173. Diederichs, M.S., Kaiser, P.K., 1999. Tensile strength and abutment relaxation as failure control mechanisms in underground excavations. Int. J. Rock Mech. Min. Sci. 36, 69–96. Feng, X., Jimenez, R., 2014. Bayesian prediction of elastic modulus of intact rocks using their uniaxial compressive strength. Eng. Geol. 173 (0), 32–40. http://dx.doi.org/10. 1016/j.enggeo.2014.02.005. Galera, J.M., Álvarez, M., Bieniawski, Z.T., 2005. Evaluation of the deformation modulus of rock masses: comparison of pressuremeter and dilatometer tests with RMR prediction. The ISP5-PRESSIO 2005 International Symposium, Madrid, Spain, pp. 1–25. Gardoni, P., Trejo, D., Vannucci, M., Bhattacharjee, C., 2009. Probabilistic models for modulus of elasticity of self-consolidated concrete: Bayesian approach. J. Eng. Mech. (ASCE) 135 (4), 295–306. http://dx.doi.org/10.1061//asce/0733-9399/2009/135:4/ 295. Geyskens, P., Der Kiureghian, A., Monteiro, P.J.M., 1993. BUMP: Bayesian Updating of Model Parameters. Tech. Rep. UCB/SEMM-93/06. Department of Civil Engineering, University of California, Berkeley, Calif. Geyskens, P., Kiureghian, A.D., Monteiro, P., 1998. Bayesian prediction of elastic modulus of concrete. J. Struct. Eng. 124 (1), 89–95. Gokceoglu, C., Sonmez, H., Kayabasi, A., 2003. Predicting the deformation moduli of rock masses. Int. J. Rock Mech. Min. Sci. 40 (5), 701–710. http://dx.doi.org/10.1016/ s1365-1609(03)00062-5. Haldar, S., Babu, G.L.S., 2008. Effect of soil spatial variability on the response of laterally loaded pile in undrained clay. Comput. Geotech. 35 (4), 537–547. http://dx.doi.org/ 10.1016/j.compgeo.2007.10.004. Hoek, E., Brown, E.T., 1980. Underground Excavations in Rock. Spon Press, London. Hoek, E., Brown, E.T., 1997. Practical estimates of rock mass strength. Int. J. Rock Mech. Min. Sci. 34 (8), 1165–1186. Hoek, E., Diederichs, M.S., 2006. Empirical estimation of rock mass modulus. Int. J. Rock Mech. Min. Sci. 43 (2), 203–215. http://dx.doi.org/10.1016/j.ijrmms.2005.06.005. Hoek, E., Carranza-Torres, C., Corkum, B., 2002. Hoek-Brown criterion-2002 edition. NARMS-TAC Conference, Toronto, pp. 267–273. Honjo, Y., Wen-Tsung, L., Sakajo, S., 1994. Application of Akaike information criterion statistics to geotechnical inverse analysis: the extended Bayesian method. Struct. Saf. 14 (1–2), 5–29. http://dx.doi.org/10.1016/0167-4730(94)90004-3. ISRM, 1975. Report of the Commission on Terminology (Lisbon). Jašarević, I., Kovačević, M.S., 1996. Analyzing applicability of existing classification for hard carbonate rock in Mediterranean area. ISRM International Symposium — EUROCK’96, Turin, Italy, pp. 811–818. Jimenez, R., Sitar, N., 2009. The importance of distribution types on finite element analyses of foundation settlement. Comput. Geotech. 36 (3), 474–483. http://dx.doi.org/10. 1016/j.compgeo.2008.05.003. Jimenez-Rodriguez, R., Sitar, N., Chacón, J., 2006. System reliability approach to rock slope stability. Int. J. Rock Mech. Min. Sci. 43 (6), 847–859. http://dx.doi.org/10.1016/j. ijrmms.2005.11.011. Johnson, J.B., Omland, K.S., 2004. Model selection in ecology and evolution. Trends Ecol. Evol. 19 (2), 101–108. http://dx.doi.org/10.1016/j.tree.2003.10.013. Kang, S.-S., Kim, H.-Y., Jang, B.-A., 2012. Correlation of in situ modulus of deformation with degree of weathering, RMR and Q-system. Environ. Earth Sci. 1–8 http://dx.doi.org/ 10.1007/s12665-012-2088-y. Kayabasi, A., Gokceoglu, C., Ercanoglu, M., 2003. Estimating the deformation modulus of rock masses: a comparative study. Int. J. Rock Mech. Min. Sci. 40 (1), 55–63. http://dx.doi.org/10.1016/S1365-1609(02)00112-0. Khabbazi, A., Ghafoori, M., Lashkaripour, G.R., Cheshomi, A., 2013. Estimation of the rock mass deformation modulus using a rock classification system. Geomech. Geoeng. 8 (1), 46–52. http://dx.doi.org/10.1080/17486025.2012.695089. Kim, G., 1993. Revaluation of geomechanics classification of rock masses. The Korean Geotechnical Society of Spring National Conference, Seoul, pp. 33–40. Li, D.Q., Tang, X.S., Phoon, K.K., Chen, Y.F., Zhou, C.B., 2013. Bivariate simulation using copula and its application to probabilistic pile settlement analysis. Int. J. Numer. Anal. Methods Geomech. 37 (6), 597–617. Lunn, D., Thomas, A., Best, N., Spiegelhalter, D., 2000. WinBUGS — a Bayesian modelling framework: concepts, structure, and extensibility. Stat. Comput. 10 (4), 325–337. http://dx.doi.org/10.1023/A:1008929526011. Lunn, D., Jackson, C., Best, N., Thomas, A., Spiegelhalter, D., 2012. The BUGS Book: A Practical Introduction to Bayesian Analysis. CRC Press, Boca Raton. Luoma, A., 2008. Bayesian model selection. Proceedings of 2008 Workshop on Information Theoretic Methods in Science and Engineering, Tampere, Finland. Miranda, T., Gomes, C.A., Ribeiro e Sousa, L., 2009. Bayesian methodology for updating geomechanical parameters and uncertainty quantification. Int. J. Rock Mech. Min. Sci. 46 (7), 1144–1153. http://dx.doi.org/10.1016/j.ijrmms.2009.03.008. Mitri, H.S., Edrissi, R., Henning, J., 1994. Finite element modeling of cable bolted slopes in hard rock ground mines. The SME Annual Meeting, Albuquerque, New Mexico, pp. 94–116. Myung, I.J., 2003. Tutorial on maximum likelihood estimation. J. Math. Psychol. 47 (1), 90–100. http://dx.doi.org/10.1016/S0022-2496(02)00028-7.

X. Feng, R. Jimenez / Engineering Geology 199 (2015) 19–27 Nicholson, G.A., Bieniawski, Z.T., 1990. A nonlinear deformation modulus based on rock mass classification. Int. J. Min. Geol. Eng. 8 (3), 181–202. http://dx.doi.org/10.1007/ bf01554041. Ntzoufras, I., 2009. Bayesian Modeling Using WinBUGS. John Wiley & Sons, Inc., Hoboken, New Jersey. Palmström, A., Singh, R., 2001. The deformation modulus of rock masses — comparisons between in situ tests and indirect estimates. Tunn. Undergr. Space Technol. 16 (2), 115–131. http://dx.doi.org/10.1016/S0886-7798(01)00038-4. Read, S., Richards, L., Perrin, N., 1999. Applicability of the Hoek-Brown failure criterion to New Zealand greywacke rocks. In: Vouille, G., Berest, P. (Eds.), The 9th International Congress on Rock Mechanics, Paris, pp. 655–660. Schwarz, G., 1978. Estimating the dimension of a model. Ann. Stat. 6 (2), 461–464. Serafim, J.L., Pereira, J.P., 1983. Considerations on the geomechanical classification of Bieniawski. International Symposium on Engineering Geology and Underground Openings, Portugal, Lisbon. SPG, pp. 1133–1144. Shen, J., Karakus, M., Xu, C., 2012. A comparative study for empirical equations in estimating deformation modulus of rock masses. Tunn. Undergr. Space Technol. 32 (0), 245–250. http://dx.doi.org/10.1016/j.tust.2012.07.004. Sonmez, H., Gokceoglu, C., Ulusay, R., 2004. Indirect determination of the modulus of deformation of rock masses based on the GSI system. Int. J. Rock Mech. Min. Sci. 41 (5), 849–857. http://dx.doi.org/10.1016/j.ijrmms.2003.01.006. Sonmez, H., Gokceoglu, C., Nefeslioglu, H.A., Kayabasi, A., 2006. Estimation of rock modulus: for intact rocks with an artificial neural network and for rock masses with a new empirical equation. Int. J. Rock Mech. Min. Sci. 43 (2), 224–235. http://dx.doi.org/10. 1016/j.ijrmms.2005.06.007. Spiegelhalter, D.J., Best, N.G., Carlin, B.P., Linde, A.V.D., 2002. Bayesian measures of model complexity and fit. J. R. Stat. Soc. Ser. B 64 (part 4), 583–639. Tang, X.-S., Li, D.-Q., Rong, G., Phoon, K.-K., Zhou, C.-B., 2013a. Impact of copula selection on geotechnical reliability under incomplete probability information. Comput. Geotech. 49, 264–278.

27

Tang, X.-S., Li, D.-Q., Zhou, C.-B., Phoon, K.-K., Zhang, L.-M., 2013b. Impact of copulas for modeling bivariate distributions on system reliability. Struct. Saf. 44, 80–90. http://dx.doi.org/10.1016/j.strusafe.2013.06.004. Verman, M., Singh, B., Viladkar, M.N., Jethwa, J.L., 1997. Effect of tunnel depth on modulus of deformation of rock mass. Rock Mech. Rock. Eng. 30 (3), 121–127. http://dx.doi. org/10.1007/bf01047388. Wang, Y., Aladejare, A.E., 2015. Selection of site-specific regression model for characterization of uniaxial compressive strength of rock. Int. J. Rock Mech. Min. Sci. 75 (0), 73–81. http://dx.doi.org/10.1016/j.ijrmms.2015.01.008. Wang, Y., Liu, Q., 2006. Comparison of Akaike information criterion (AIC) and Bayesian information criterion (BIC) in selection of stock–recruitment relationships. Fish. Res. 77 (2), 220–225. http://dx.doi.org/10.1016/j.fishres.2005.08.011. Wang, Y., Huang, K., Cao, Z., 2013. Probabilistic identification of underground soil stratification using cone penetration tests. Can. Geotech. J. 50 (7), 766–776. http://dx.doi. org/10.1139/cgj-2013-0004. Wang, Y., Huang, K., Cao, Z., 2014. Bayesian identification of soil strata in London clay. Géotechnique 64 (3), 239–246. http://dx.doi.org/10.1680/geot.13.T.018. Zhang, L., 2010. Estimating the strength of jointed rock masses. Rock Mech. Rock. Eng. 43 (4), 391–402. http://dx.doi.org/10.1007/s00603-009-0065-x. Zhang, L., Einstein, H.H., 2004. Using RQD to estimate the deformation modulus of rock masses. Int. J. Rock Mech. Min. Sci. 41 (2), 337–341. http://dx.doi.org/10.1016/ S1365-1609(03)00100-X. Zhang, L., Tang, W., Zheng, J., 2004. Reducing uncertainty of prediction from empirical correlations. J. Geotech. Geoenviron. 130 (5), 526–534. http://dx.doi.org/10.1061/ (ASCE)1090–0241(2004)130:5(526).