Tunnelling and Underground Space Technology

Contents lists available at ScienceDirect

Tunnelling and Underground Space Technology

journal homepage: www.elsevier.com/locate/tust

CrossMark

Stochastic assessment of pillar stability at Laisvall mine using Artificial Neural Network

Musa Adebayo Idrisa'*, David Saianga,b, Erling Nordlunda

a Division of Mining and Geotechnical Engineering, Lulea University of Technology, SE-971 87 Lulea, Sweden b SRK Consulting (Sweden) AB, SE-931 31 Skelleftea, Sweden

ARTICLE INFO

Article history:

Received 2 May 2014

Received in revised form 29 April 2015

Accepted 8 May 2015

Available online 29 May 2015

Keywords:

Stochastic assessment Probability of failure Reliability index Artificial Neural Network Pillar stability

ABSTRACT

Stability analyses of any excavations within the rock mass require reliable geotechnical input parameters such as in situ stress field, rock mass strength and deformation modulus. These parameters are intrinsically uncertain and their precise values are never known, hence, their variability must be properly accounted for in the stability analyses. Traditional deterministic approaches do not quantitatively consider these uncertainties in the input parameters. To incorporate these uncertainties stochastic approaches are generally used. In this study, a stochastic assessment of pillar stability using Artificial Neural Network (ANN) is presented. The uncertainty in the rock mass properties at the Laisvall mine were quantified and the probability density function of the deformation modulus of the rock mass was determined using probabilistic approach. The variability of the in situ stress was also considered. The random values of the deformation modulus and the horizontal in situ stresses were used as input parameters in the FLAC3D numerical simulations to determine the axial strain in the pillar. ANN model was developed to approximate an implicit relationship between the deformation modulus, horizontal in situ stresses and the axial strain occurring in pillar due to mining activities. The closed-form relationship generated from the trained ANN model, together with the maximum strain that the pillar can withstand was used to assess the stability of the pillar in terms of reliability index and probability of failure. The results from this study indicate that, the thickness of the overburden and pillar dimension have a substantial effect on the probability of failure and reliability index. Also shown is the significant influence of coefficient of variation (COV) of the random variables on the pillar stability. The approach presented in this study can be used to determine the optimal pillar dimensions based on the minimum acceptable risk of pillar failure.

© 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CCBY-NC-ND license

(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

A pillar can be defined as the in-situ rock mass between two or more underground openings. It is the main support in room and pillar mines. The support provided by the pillars controls the rock mass displacement throughout the zone of influence of mining, while the mining proceeds.

The analysis and design of mine pillars generally seek to optimize the size of the pillars so as to maximize the extraction ratio (i.e. amount of ore extracted relative to the total amount of ore available) while maintaining the stability of the mine. Hence the design of pillars has both economic and safety implications. The knowledge of the pillar strength and the determination of the required safety factor for a given loading condition are the

* Corresponding author. Tel.: +46 736462465. E-mail addresses: idris.musa@ltu.se, abumariyam20@gmail.com (M.A. Idris).

most important aspects of pillar design. Conventional pillar design methods comprise the calculation of the mean pillar stress (e.g. the tributary area method and the method by Coates (1981) and the estimation of the pillar strength using empirical formulae (e.g. Obert and Duvall, 1967; Krauland and Söder, 1987; Sjöberg, 1992). Based on the stress and strength of the pillar the factor of safety can be calculated. The factor of safety is the ratio of the pillar strength to the induced stress in the pillar and the pillar fails when the ratio is less than 1.

Though the conventional methods are widely used for pillar design, Alber and Heiland (2001) have expressed some concerns about this conventional approach for pillar design at shallow depth. They observed that the pillar failure at shallow depth could not be properly explained by comparing pillar strength with stresses induced on the pillar by mining activities. They suggested amongst other approaches that pillar failure could be related to strain. Therefore, when considering the strain occurring in the

http://dx.doi.org/10.1016/j.tust.2015.05.003 0886-7798/© 2015 The Authors. Published by Elsevier Ltd.

This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

pillar the factor of safety can be determined as the ratio of the maximum strain that a pillar can withstand to the strain occurring in the pillar due to mining activities. Nevertheless, either ways of determining the factor of safety are largely deterministic and do not consider the inherent variability of the rock mass properties and that of the in situ stress field. Mean values of these input parameters are generally assumed. The results from the deterministic approach could be misleading depending on the distributive character of the rock property variation (Kim and Gao, 1995). Deng et al. (2003) have reported instances where pillars failed despite the fact that the failed pillars had been considered stable with factor of safety greater than 1.

Therefore, for a reliable design and analysis of construction elements such as mine pillars appropriate methods which incorporate the variability in the rock mass properties must be used. The methods which consider this variability are known as stochastic or probabilistic methods. With a stochastic approach, the stability analysis can be considered as a random system, where the occurrence of a pillar failure is a random event depending on the outcome of the random variables involved.

A number of stochastic approaches have been applied to various geotechnical problems, including underground excavation problems (e.g. Chen et al., 1997; Lilly and Li, 2000; Cai, 2011; Idris et al., 2011; Dohyun et al., 2012), tunnel support (e.g. Schweiger et al., 2001; Li and Low, 2010; Oreste, 2005), subsidence (e.g. Torano et al., 2000) and pillar stability (e.g., Pine, 1992; Joughin et al., 2000; Griffiths et al., 2002; Deng et al., 2003; Cauvin et al., 2009; Najafi et al., 2011; Recio-Gordo and Jimenez, 2012; Wattimena et al., 2013). Pine (1992) presented a probabilistic approach for pillar design whereby normal probabilistic distributions were assumed for the random variables and the safety margin. Joughin et al. (2000) employed the point estimate method (Rosenblueth, 1981) to account for rock strength variability in the probabilistic method they presented for the design of chromite pillars in South Africa. Griffiths et al. (2002) analysed the stability of underground pillar by using random field theory with elasto-plastic finite element algorithm in a Monte Carlo framework. Deng et al. (2003) presented a probabilistic mine design method which combines the finite element methods, neural network and reliability analysis. Cauvin et al. (2009) used probabilistic approach to assess the effect of uncertainty in mining pillar stability analysis. Najafi et al. (2011) utilized First Order Second Moment (FORM) and Advanced Second Moment (ASM) for the probabilistic stability analysis of chain pillar in a coal mine in Iran. A probabilistic model based on linear classifier theory to predict the behaviour of pillar in longwall and retreat room and pillar mining was presented by Recio-Gordo and Jimenez (2012). Wattimena et al. (2013) employed logistic regression to predict the probability of coal pillar stability for given pillar geometry and stress condition.

In general, stochastic assessment of pillar stability is performed by two procedures: the first step is to quantify the uncertainty in the rock mass properties in order to determine the basic statistical parameters (i.e. mean and variance) and probability density functions (PDFs) of the strength and deformation modulus of the rock mass using the Monte Carlo method. The Monte Carlo (MC) simulation technique is often adopted in the geotechnical stochastic analyses with implicit or explicit solutions but when the analysis is associated with numerical modelling then the MC simulation technique becomes time consuming and less appealing.

In the second step, the probability of failure is determined with respect to a specific failure criterion, which can either be the induced pillar stress exceeding the pillar strength or the strain occurring in the pillar exceeding the defined threshold strain value for the pillar. The onset of failure in the context of this study is defined as the limit state when the peak strength of the pillar is exceeded or the strain occurring in the pillar exceeds the peak

strain for the pillar. For underground excavations this limit state is not known explicitly, instead numerical analysis using the finite difference method (FDM) or the finite element methods (FEM) can be combined with function approximation tools to construct a closed-form expression for the limit state surface. Recently, many function approximation tools have been proposed such as the response surface method (RSM), the point estimate method (PEM), and the Artificial Neural Network (ANN) to model the relationship. ANN, due to its high performance, has been one of the tools used in geotechnical engineering to model the relationship between non-linear multivariate variables (Sonmez et al., 2006).

In this study, a stochastic approach was used to analyse the pillar stability at the Laisvall mine in Sweden while considering the variability in the rock mass properties and in the in-situ stresses. The uncertainty in the rock mass properties at the Laisvall mine were quantified and the probability density function of the deformation modulus of the rock mass was determined. Also the variability of the horizontal in situ stresses was considered. The random values of these parameters (i.e. deformation modulus and horizontal in situ stresses) were used as input parameters for the FLAC3D (Itasca, 2012) analyses to determine the axial strain in the pillar. The ANN model was developed to approximate an implicit relationship between the deformation modulus, horizontal stresses and the pillar axial strain within the range of possible values of the random input parameters. The closed-form relationship generated from the trained ANN model together with critical axial strain, which the pillar can withstand, was used to define a pillar performance function. The performance function was used to assess the stability of the pillar in terms of probability of failure and reliability index.

1.1. Artificial Neural Network (ANN)

The ANN, also referred to as neural network, is an information system that imitates the behaviour of the human brain by emulating the operation and connectivity of the brain to generate a general solution to a problem. ANN can be used to extract patterns and detect trends from problems where the relationship between the inputs and outputs are not sufficiently known. In recent years, ANN has been frequently used for functions approximation in different fields of science, including geotechnical engineering (Sahin et al., 2001). Basically, ANN consists of simple interconnected nodes or neurons as shown in Fig. 1 where p is the input, w is the weight, b is the bias, f is the transfer function and a is the output.

If the neuron has N number of inputs then the output a can be calculated as:

a = f WiPi + b) (1)

There are different types of transfer functions that can be used in ANN such as hard limit transfer function, linear transfer function, log-sigmoid transfer (Beale et al., 2012). The choice of the transfer function depends on the specification of a problem that the neuron is attempting to solve (Beale et al., 2012).

The architecture of ANN consists of the number of layers, the number of neurons in each layer and the neuron transfer functions. Two or more neurons can be combined in a layer and a network

/ fl }

Fig. 1. Simple structure of ANN model.

could contain one or two layers whereby each layer has different roles. There are output layers, input layers and intermediate layers or hidden layers. There is always one output layer and one input layer while there can be many hidden layers for an ANN. However, it is known that a network with one hidden layer can approximate any continuous function provided with sufficient connection weights (Hornik et al., 1989). The number of neurons in the input and output layer are determined by the number of the model input and output variables. There is no any exact guide for determining the number of neurons in the hidden layer, however some researchers (e.g. Aldrich and Reuter, 1994; Kanellopoulas and Wilkinson, 1997; Seibi and Al-Alawi, 1997) have proposed heuristic relations for determining the neuron size.

Before an ANN can be used to make projections or predictions it has to be trained. In the training or learning process, the network is presented with a pair of training datasets including input and corresponding target values. The network computes its own output using its initial weights and biases. Then, the weights and biases are adjusted iteratively to reduce the errors between the network output and the target output. Mean square error is used as error index during the training phase to improve the network performance (Tawadrous et al., 2009). One of the most commonly used learning algorithms in geotechnical engineering is the back-propagation algorithm (BPA) proposed by Rumelhart et al. (1986). In BPA there are two phases: forward prediction which calculates the output values of the ANN from training data and error back-propagation which adjusts the weight. There are many techniques to adjust the error but the steepest descent method is often used (Tawadrous et al., 2009). Once the training process of ANN is completed predictions can be made for a new input dataset.

2. Description of the case studied

The Laisvall mine in northern Sweden, although closed in 2001, was a lead-zinc mine owned and operated by Boliden Mineral AB. The orebody, hosted in flat-lying quartzitic sandstone interlayered with clayey sandstones, was mined using the room and pillar mining method. The annual production before the closure of the mine was 1.6 Mt. There were four main ore bodies in the mine, most of them situated in the lower sandstone, except the Nadok orebody which was entirely situated in the upper sandstone. The thickness of the overlying strata above the Nadok orebody varied between 110 and 300 m (Söder and Krauland, 1990). In most of the orebod-ies of the mine the roof consisted of the same type of sandstone as the ore itself.

Söder and Krauland (1990) conducted a full-scale pillar test between 1983 and 1988 in the Nadok orebody. This test comprised 9 pillars and the objective was to determine the stress level in the pillar at failure (i.e., the bearing capacity) to serve as input to pillar design and long term planning of the mine. The pillars were subjected to increasing stresses by decreasing the cross-sectional area of the pillars through blasting of slices of approximately 0.4 m thickness, thereby reducing the widths and lengths of the pillars in each of the mining steps. The pillars were initially 4.6 m high, 7.4 m wide and 8.1 m long. Pillar failure appeared to have occurred after six slices were blasted.

3. Stochastic estimation of rock mass properties using GSI system

Inputs for the rock mass parameters in the numerical simulation were determined stochastically from the GSI (Geological Strength Index) rock mass characterization system. The GSI system, developed by Hoek et al. (1995), is based on the description of the rock structure and the block surface condition. It is a system that

provides a complete set of mechanical properties for design purposes when used in the generalized Hoek-Brown criterion (Cai et al., 2007). The mechanical rock mass properties are the deformation modulus (Em), the Hoek-Brown strength parameters (mb, s and a), the tensile strength, and the equivalent shear strength parameters: cohesion (c) and the internal friction angle (/). When the GSI system is used for rock mass characterization, two groups of parameters need to be determined, which includes the uniaxial compressive strength (UCS) and the material constant (mi) of intact rock. The other group is the joint parameters which consist of joint geometry and strength parameters (Cai, 2011). These parameters varied considerably for the rock mass where the pillar test was conducted as reported by Söder and Krauland (1990), Palmström (1995), Exadaktylos and Stavropoulou (2008).

3.1. Evaluation of GSI from Vb and Jc

The determination of GSI based on visual observation may be subjective. Recently, Cai et al. (2004b) proposed a new approach to quantitatively determine GSI based on the block size and condition; block volume (Vb) and the joint condition factor (Jc). This approach is well suited for the stochastic approach and therefore adopted for this study.

The block size is determined from the joint spacing, joint orientation, number of joint sets and joint persistence. According to Cai et al. (2004b) the effect of the intersection angle between joint sets is relatively small compared to the joint spacing hence they suggested that for practical purposes the block volume for three or more joints can be approximated using

Vb = S1S2S3

where Si is the spacing of each joint set.

The joint surface condition (Jc) is defined by the roughness, weathering and infilling. It is similar to the factor used by Palmström (1995) to quantify the joint surface condition. The combination of these factors defines the strength of a joint or block surface. Cai et al., 2004b defined Jc as:

j _JwJs Jc ~ J Ja

where JW and JS are the large-scale waviness and small-scale smoothness, respectively (Barton and Bandis, 1990; Palmström, 1995) and JA is the joint alteration factor (Barton et al., 1974).

Once the value of Jc and Vb are known the GSI value can be determined using the equation proposed by Cai and Kaiser (2006) which is expressed in terms of Jc and Vb, thus

26.5 + 8.79 ln Jc + 0.9ln Vb 1 + 0.0151 lnJc - 0.0253 ln Vb

By extension, since the variability in the material properties is considered in this study, if the mean values and the coefficient of variations of Jc and Vb together with their probability density functions (PDF) are known then they could be used as input variables to determine the distribution of the GSI using the Monte Carlo simulation.

Palmström (1995) has conducted extensive evaluation of the joint characteristics of the pillars at the Laisvall mine and reported the Vb and Jc values for the rock mass at the mine. Exadaktylos and Stavropoulou (2008), based on the results of Palmström (1995), presented the range of values for the Vb and Jc, respectively, for the quarzitic sandstone pillar at the Laisvall mine. The mean values and the coefficient of variation (COV) for both Vb and Jc were estimated based on the ''Three-sigma rule'' described by Dai and Wang (1992). The three-sigma rule is based on the fact that 99.73% of all values of a normally distributed parameter fall within three

standard deviations of the mean. The normal distribution is assumed for both Vb and Jc and the distributions were truncated at their maximum and minimum values. The mean values and the COV for the Vb and Jc together with their respective truncated PDF were used as inputs in Eq. (4) to generate the distribution of the GSI using the Monte Carlo simulation technique available in the Excel add in program @RISK (Palisade, 2001). Table 1 shows the statistical parameters for Vb and Jc and statistical parameters for the GSI generated from the Monte Carlo simulation.

3.2. Determination of statistical parameters for UCS and, mi

The GSI value, the UCS and the mi of the intact rock are basic inputs for the generalized Hoek-Brown failure criterion (Hoek et al., 2002). The UCS and mi of the intact rock are determined by laboratory testing as described by Hoek and Brown (1997). In this study, the mean value and the range for UCS used for the rock types in the Nadok orebody have been reported by Soder and Krauland (1990), Palmstrom (1995). The information about the variation of the UCS is not reported. A normal distribution with a COV of 15% was assumed for the UCS. Generally, the COV of 1520% is applicable to most natural geotechnical materials (Rethati, 1988). A normal distribution with the mean value of 17 and COV of 10% was assumed for the mi. This was taken from the value suggested for sandstone by Marinos and Hoek (2000). Table 2 shows the statistical parameters for the UCS and mi used in this study.

3.3. Determination of the rock mass properties and their Probability Density Functions

The generalized Hoek-Brown criterion for jointed rock masses was employed to estimate the rock mass strength parameters and deformation modulus for the Nadok orebody. At failure the generalized Hoek-Brown criterion relates the major effective principal stress r to the minor effective principal stress r3 as follows

Table 2

Statistical parameters for the UCS and mi.

ri = Ö3 + rci[mb — -rci

where mb, s, a are Hoek-Brown strength parameters for the rock mass and rciis the UCS of the intact rock. Having determined the statistical parameters (i.e. mean and COV) of GSI and mi Monte Carlo simulation technique was used to determine the distributions of the Hoek-Brown strength parameters for the orebody using the following equations:

mb = mi exp

s = exp

/GSI - 100\ V28 - 14D)

/GSI - 100^ I 9 - 3D

a = 2 +1 (e-GSI/15 - e-20/3)

where D is a factor which depends on the degree of disturbance the rock mass has been subjected to by blast damage and stress relaxation. It varies from 0 for undisturbed in situ rock masses to 1 for very disturbed rock masses. For this study, blasting damage was assumed to be minimal, hence D = 0 was used.

Table 1

Statistical parameters for Vb, Jc and GSI.

Parameters Min. value Max. value Mean COV (%) Distribution

Jc 0.75 2.25 1.5 16.7 Normal

Vb (x103 cm3) 100 300 200 16.5 Normal

GSI 51 63.54 58.7 4.0 Normal

Parameters Mean value COV (%)

UCS (MPa) 210 15

mi 17 10

Different equations have been proposed to determine the deformation modulus of the rock mass based on the range of the UCS (Hoek et al., 2002). Hoek and Diederichs (2006) suggested that when reliable estimates of Young's modulus for intact rock are available, the deformation modulus of the rock mass (Em) can be calculated from:

1 - D=2

1 e((60+15D-GSI/11))

where Ei is the Young modulus and it was determined from the modulus ratio (MR) proposed by Deere (1968) thus:

Ei = MRrd(10)

The value of MR was determined from the mean values of UCS and Ei for the intact rock reported by Söder and Krauland (1990). Eqs. (9) and (10) were used to estimate the statistical parameters and the PDF for deformation modulus (Em) for the rock mass by using the already determined PDFs of the UCS and GSI. Similar to the previous parameters, the PDF and the statistical parameters were generated using the Monte Carlo simulation technique available in @RISK. Table 3 shows the probability density function (PDF), mean, standard deviation (STDEV) and the coefficient of variation (COV) for the deformation modulus (Em).

4. Variability in in-situ stress

The knowledge of in situ stresses is essential for the design of underground excavations and construction elements such as mine pillars. Due to spatial variation of the in-situ stress field at various locations within the rock mass, there is always an uncertainty in the results of in-situ stresses measurements (Cai, 2011).

In-situ stresses are usually reported as; vertical stress (rv), maximum horizontal stress (rH), and minimum horizontal stress (rh). The vertical stress is normally assumed to vary linearly with depth as:

rv = pgz

where q is the rock density, g is the gravity and z is the overburden. For this study, the rock density for the rock mass was assumed to be constant with a value of 2700 kg/m3 and the overburden varied from 110 to 300 m.

The maximum horizontal stress (rH) and minimum horizontal stress (rh) recorded at the mine varied between 20 MPa and 25 MPa (Marklund, 2013). Normal distribution was assumed for the horizontal in situ stresses with coefficient of variation (COV) of 3.7% using the three sigma rule. Only the variability in the horizontal stresses was considered in this study.

5. Numerical analysis

Since the objective of the numerical analysis was to determine the induced axial stress and the elastic axial strain occurring in the

Table 3

PDF, mean, standard deviation and coefficient of variation (COV) for Em.

Parameter Mean STDEV COV (%) PDF

Em (GPa) 23.47 4.04 17.17 Normal

Em — E

pillar due to the mining, linear elastic analyses were considered sufficient. Series of elastic numerical analyses were carried out using the FLAC3D code (Itasca, 2012) considering the possible variability in the rock mass properties and the in-situ stress conditions.

5.1. Model geometry

The layout of the test area and the geometry of the simulated pillars are shown in Fig. 2. Roller boundary conditions were used along all external boundaries. The mesh size used for the model was 0.4 m x 0.4 m x 0.4 m. It was chosen based on the results from the FLAC analysis where the induced axial stress in pillar was monitored for different mesh sizes until further reduction in the mesh size gave no significant difference in the model result.

5.2. Modelling sequence

The modelling sequence used in this study consists of the following steps:

(i) Generation of the model grid using a brick-shaped mesh.

(ii) Selection of appropriate boundary conditions and material model.

(iii) Initialization of in-situ stresses and material properties.

(iv) Stepping of the model to initial equilibrium prior to any excavations.

(v) Excavation of the roadways to create the rib pillar and later excavating the rooms across the rib pillar to create the panel pillar. These excavations were henceforth referred to as mining step 0 in the modelling sequence.

(vi) Sequential loading of the pillar induced by blasting the pillar in 6 steps, henceforth referred to as mining steps 1-6. Each step extracts one slice of approximately 0.4 m from the length of the pillar and one slice of approximately 0.4 m from the width of the pillar.

Fig. 3 shows the mining steps 1-6 and Table 4 shows the pillar dimensions corresponding to the mining steps. At each mining step the axial stress and maximum vertical displacement at both ends of the pillar were obtained. The mean pillar stress was estimated by averaging the vertical stress component computed by FLAC3D along the center line ab (see Fig. 4) of the pillar. The axial strain (e) was calculated as follows:

e = x 100 (12)

where Ua and Ub are the vertical displacement at points a and b, respectively while H is the height of the pillar.

5.3. Sensitivity analysis

In the room and pillar mining method, roadways are usually developed ahead of mining and are initially subjected to the in-situ stress state. As mining progresses, the pillars are subjected to an increased induced stress which consequently leads to

Table 4

Mining steps and corresponding pillar dimension.

Mining steps Height, H (m) Width, W (m) Length (m) W/H (approx.)

0 4.6 7.4 8.1 1.61

1 4.6 6.7 7.8 1.46

2 4.6 6.3 7.2 1.37

3 4.6 5.9 6.6 1.28

4 4.6 5.5 6.2 1.20

5 4.6 4.9 5.9 1.07

6 4.6 4.5 5.3 0.98

increase in the rock mass response in terms of displacement or strain. Traditionally, the stability of a pillar is assessed by using the criteria which compares the strength of the pillar to the maximum average stress in the pillar. The maximum average stress can be determined from two or three dimensional elastic numerical analyses while the strength is determined from empirical equations. The pillar stability can also be assessed by comparing the strain monitored in the pillar during the numerical analyses with the maximum strain that the pillar can withstand. Some of the input parameters for the elastic numerical analyses are the deformation modulus of the rock mass and horizontal in-situ stress. As discussed previously, these parameters exhibit inherent variability and uncertainty hence for this study there is need to find the appropriate criteria for assessing the stability of the pillar while considering the variability of these parameters. A sensitivity analysis was performed to determine the effect of the variability of these input parameters on the axial strain of and the mean axial stress in the pillar. Two cases were considered for the sensitivity analysis as follows:

(i) The deformation modulus was varied between 20 and 32 GPa while the other input parameters were kept constant and equal to their mean values.

(ii) The horizontal stresses were varied between 20 and 26 MPa while the other input parameters were kept constant and equal to their mean values.

The mean values of the input parameters used for the sensitivity analysis are shown in Table 5.

5.3.1. Results of the sensitivity analysis

The results of the sensitivity analysis are shown in Figs. 5 and 6. It is evident from Fig. 5 that the variation in the horizontal stress has significant effect on both axial strain and mean axial stress in the pillar. The strain and the mean stress increased when the width-to-height ratio of the pillar was reduced from 1.61 to 1.46 (i.e. mining step 0 to step 1). However, as shown in Fig. 6a the axial strain in the pillar decreases for an increase in the deformation modulus. Also the strain increases with reduction in the width-to-height ratio. Therefore, from Fig. 6a it is evident that the axial strain in the pillar is sensitive to the variations in the deformation modulus as well as in the pillar geometry. Conversely, the mean axial stress is not sensitive to any increase in the deformation modulus as shown in Fig. 6a. However the axial stress increased when the width-to-height ratio was reduced (i.e. from mining step 0 to step 1). These results show that the distribution of elastic stress in the pillar is a function of the pillar geometry and the in-situ stress. This has also been observed by Martin and Maybee (2000), Lunder and Pakalnis (1997).

The sensitivity analysis shows that given the same pillar geometry the stress in the pillar, under elastic conditions, will be the same for pillars with different rock mass properties if they are subjected to the same in situ stresses. Deformation (e.g. strain) is the response of the rock mass when subjected to stress/load hence this response (i.e. strain) can be used to analyse the stability of the pillar instead of using a stress/strength criterion. Therefore, in this

Table 5

Mean values of the input parameters for the sensitivity analysis.

Parameters Mean values

Deformation modulus Em, GPa 23.47

Maximum horizontal stress, oh, MPa 22.50

Maximum horizontal stress, rh, MPa 22.50

Vertical stress, ov MPa 5.43

Simulated

15.5 m (b)

16.2 m

Fig. 2. The layout of the test area (a) and the geometry of the simulated pillar (b).

Direction of mining steps Fig. 3. Plan view of the simulated pillar showing the mining steps (not to scale).

study strain-based criterion is used for assessing the pillar stability in this study.

5.4. Determination of the pillar performance Junction

In order to determine the probability of failure for the pillar as well as its reliability the relationship between the capacity and the demand on the pillar must be established. Concluding from the results of the sensitivity analysis presented above, the strain based criterion is adopted for this study. Therefore, the pillar capacity is expressed in terms of the critical strain (ecriticai) and the demand is the axial strain (e).

Sakurai (1981, 1997, 1999) has developed a concept of critical strain for rock masses. The critical strain is the ratio of uniaxial compressive strength to Young's modulus i.e.

^critical :

where ecritical is the critical strain, rcm is the uniaxial compressive strength (UCS) of the rock mass and Em is the deformation modulus of the rock mass. The UCS of the rock mass rcm is obtained by setting r3 = 0 and a1 = acm in Eq. (5), giving:

a cm = rcisa (14)

The rock mass constants a and s are obtained from Eqs. (7) and (8), respectively.

Since the critical strain may not be the strain at failure for most rocks, Sakurai (1981) proposed a relationship between strain at failure and critical strain as

^failure

^critical

where £faiiure is the failure strain and Rj is a parameter representing failure strength. According to Sakurai (1981) Rj ranges from 0.05 to 0.8. The range of Rj suggested by Sakurai does not have a general application as it includes both weak and hard rocks. However, based on the study conducted by Cai et al. (2004a) on the generalized crack initiation and propagation thresholds of brittle rock masses Cai (2011) proposed the Rj ranging from 0.1 to 0.3 for hard rock masses. For this study 0.2 was assumed to be an approximated value for Rj for the rock mass at the Laisvall mine based on a rough estimate using the available data from the mine.

Therefore, for this study the performance function for the pillar is the limit-state function which defines the relationship between the failure strain and the axial strain due to mining. The limit state function, g(X), and probability of failure Pjcan be defined as:

g(X)=R(X)-S(X) (16)

Pf = P[g(X) < 0]

where R is the failure strain and S is the axial strain. Both R(X) and S(X) are PDFs of the random variable X.

Aside from the limit state function, the safety index or reliability index (b) is another important quantity for the reliability of a pillar. The reliability index can be defined as:

b = 7T rm

where im and am are the mean and standard deviation of the PDF of g(X), respectively. The larger the reliability index, the smaller the probability of pillar failure. The limit state function, g(X), is assumed to follow a normal distribution hence b can be defined as

(13) b = -/-1(Pf ) = /-1(1 - Pf )

where / 1 is the inverse cumulative density function (CDF) of the standard normal variable at the probability level (1 - Pf).

5.4.1. Determination of closed-form equation for axial strain

In this study an approach that integrates ANN and finite difference analysis was adopted to obtain a closed-form equation of the

axial strain for each of the mining steps. The proceeding sections describe the processes involved in achieving this.

5.4.1.1. Determination of model training dataset. An important step in developing ANN models is to select the model input variables that have the most significant effect on the model performance (Faraway and Chatfield, 1998). For this study, the training datasets are the deformation modulus and horizontal in-situ stresses as input data and the corresponding axial strain as output or target data. In order to approximate the relationship between the axial strain and the random variables a training dataset is designed. Thirty-six (i.e. 62) combinations of random values of Em and OH were generated by considering equally spaced six sampling points of each random values within the range of ±3 standard deviations of their respective means. This represents almost the total population of their normal distributions.

The FLAC3D elastic simulations were carried out using each of the 36 combinations of Em and OH together with constant values for Poisson's ratio (0.2) and mean value for overburden (i.e. 205 m) to generate the corresponding 36 axial strains. The 36 datasets were divided into 3 sections randomly for the ANN model training, testing and validation. 70% of the datasets for the training, 15% for the testing and 15% for the validation.

The input and output data were normalized or scaled to be in the range of ±1 because the ANN algorithm especially the back-propagation, works best when the training datasets are scaled or normalized (Beale et al., 2012). The following equation was used for linear scaling of the input and target data (Goh and Kulhawy, 2003).

Xnorm = 2(* ' X°M - 1 (20)

\xmax xmin /

where xnorm is the normalized value of parameter x with maximum and minimum values of xmax and xmin, respectively. The ranges of values used for Em and OH as input parameters for the training datasets were 23.47 ± 4.04 GPa and 22.50 ± 0.83 MPa, respectively. The corresponding maximum and minimum values of the output values of the axial strain for the base case (i.e. at depth of 205 m) are shown in Table 6.

5.4.1.2. Designing the architecture of the ANN. The multilayer per-ceptron (MLP) feed-forward neural network was considered for this study. MLP is commonly used in geotechnical engineering (Sahin et al., 2001). For the development of the network, a commercial software package MATLAB (MathWorks Inc., 2012) was used to simulate the ANN operations. The network has one input layer, one output layer and one hidden layer. The input layer has two neurons and the output layer has seven neurons. The number of neurons in the hidden layer was chosen to be 5 utilizing a trial and error approach. The best configuration was achieved in this work by using the linear transfer function (purelin) in the output layer and the tansig transfer function in the hidden layer. The

Horizontal in-situ stress (MPa) Horizontal in-situ stress (MPa)

(a) (b)

Fig. 5. Effect of varying the horizontal stresses on: (a) axial strain and (b) mean axial stress in the pillar at mining steps 0 and 1.

Table 6

Maximum and minimum values of axial strain.

Mining step Axial strain (%) Max. Min.

0 0.083 0.037

1 0.090 0.040

2 0.095 0.042

3 0.100 0.044

4 0.102 0.045

5 0.111 0.049

6 0.116 0.051

network was trained with Levenberg-Marquardt back-propagation algorithm (trainlm).

5.4.1.3. Model training, testing and validation. The process of optimizing the connection weights is known as training (Sahin et al., 2001). The weights and biases were initialized to non-zero random values. Then the normalized training dataset of the input and output values generated from the elastic simulations with FLAC3D were presented to the network. The performance of the ANN model was measured in terms of an error criterion between the target output and the calculated output. The output calculated at the end of each feed-forward computation was compared with the target output to estimate the mean-squared error (MSE). Thereafter the back-propagation algorithm adjusted the weights and biases until the mean squared error was greatly minimized.

Fig. 7 shows the relationship between the output targets (axial strain values from FLAC3D simulations) and the predicted values obtained from the ANN training, testing and validation process for the base case. Each diagram in the figure represents plots for different percentages of the total dataset, for training (70%), validation (15%), test (15%) and all (100%). The model shows good correlation to the training, validation, test and of course all the datasets with R values greater than 0.9. R values is an indication of the relationship between the outputs and targets where R = 1 indicates exact linear relationship.

5.4.1.4. Implementation of the trained network. Having trained the ANN by optimizing the weights and biases, the network can be used to predict target values given the input pattern within the range of the input values used for the network training.

The optimal trained connection weights and biases were extracted from the network model and used to develop a mathematical expression relating the input values, Em and OH and the output variables, axial strain (e), for each mining step. The network model has two input neurons, five hidden neurons, and seven output neurons, each one for each mining step hence the closed-form equation relating the inputs variables and output variables is as follows:

Tn = f 2 № + LW (f 1 (B1 + IW.Pn))} (21)

where B1 is the bias vector for the input, IW is the vector for weight connection between neurons of the hidden layer and the single output, B2 is the vector for the bias at hidden layer neuron, LW is the matrix for the hidden layer weight. Pn is a matrix of normalized input vectors used as input values for the network training. Tn is the corresponding matrix of normalized output vectors and f1 and f2 are the transfer functions which are tansig and purelin, respectively. The normalized output vectors (Tn) could be converted to non-normalized output vectors (T) as follows:

T = 0.5(Tn + 1)(max T - min T) + min T (22)

where maxT and minT are the vectors containing the maximum and minimum values of e, respectively which were used for the network training and testing. These values for the base case have been presented in Table 6.

5.4.2. Pillar performance function

Using Eqs. (21) and (22), the pillar performance function can be expressed as:

g(X) = efailure - T. (23)

The PDFs and the statistical parameters of Oci, s, a and Em were used as input values to determine the random values of the critical strain (ecriticai) as shown in Eqs. (13) and (14). The strain at failure,

Training: R=0.99999

Validation: R=0.99999

> ■o

-1 -0.5 0 0.5

Target value of normalized axial strain

> ■o

-1 -0.5 0 0.5

Target value of normalized axial strain

Test: R=0.99998

All: R=0.99998

-1 -0.5 0 0.5 1

Target value of normalized axial strain

Target value of normalized axial strain

Fig. 7. Comparison of the predicted values of the normalized axial strain using ANN and the target values (obtained from the FLAC3D simulations) for the base case: (a) training dataset (70% of the total dataset); (b) test dataset (15% of the total dataset); (c) validation dataset (15% of the total dataset); and (d) all the data sets.

which is the peak strain, was then determined using Eq. (15). The PDF and the statistical parameters of Em and aH were used as the input values to generated T (i.e.e) in Eq. (22). The Monte Carlo (MC) technique was used to simulate the performance function (i.e. Eq. (23)) using the program @RISK. In each simulation, the values of the variables X (i.e. aci, s, a, Em and aH) were randomly generated according to their PDFs. The PDFs, means and standard deviations of the performance function (g(X)) were determined from the MC simulation after 105 simulations at each mining step. The probability of failure Pf for the pillar at each mining step was calculated as

Pf = P[g(X) < 0]= N

where Nj is the number of simulations with g(x) < 0 (i.e. number of pillar failures), and N is the total number of simulations. The reliability index (b) was also estimated from the statistical parameters of the performance function using Eq. (18).

6. Results and discussions

In this section the results of the stochastic assessment of the pillar behaviour at the Laisvall mine are presented and discussed below.

6.1. Effect of mining activities on the pillar stability

The effect of the mining steps or the gradual reduction of the pillar width-to-height ratio on the stability of the pillar is demonstrated with the base case (i.e. at depth 205 m) and presented in Fig. 8. It is obvious from the figure that the mining steps have significant effect on the pillar stability. By reducing the

cross-sectional area of the pillar the pillar strength is reduced thereby allowing more strain to occur in the pillar. Consequently the probability of pillar failure significantly increased as shown in Fig. 8, while the pillar reliability index decreased as the width-to-height ratio reduced.

6.2. Effect of mining depth

As mentioned previously, the overburden varied between 110 and 300 m. In order to determine the effect of the various mining

—•— Probability of pillar failure —♦— Pillar reliability index 1.0 -

° 0.2 -| Ph

0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 Width-height ratio (W/H)

Fig. 8. Effect of mining activities on the pillar stability in terms of failure probability and reliability index (base case).

depths additional probabilistic analyses following the same procedure for the base case were performed. The mining depths considered were 220, 240, 260 and 280 m. Fig. 9 shows the effect of the mining depth on the performance of the pillar. It can be seen that the probability of pillar failure increases as the depth increases. Likewise, the corresponding pillar reliability index decreases as the mining depth increases; the lower the reliability index, the larger the corresponding probability of pillar failure. The failure probabilities of the pillar at all the considered depths are very small up to step 4 but increased significantly after step 4 for all the depths.

6.3. Effect of changing the coefficient of variation (COV) of the random variables on pillar stability

The UCS of the intact rock is a factor in the determination of the deformation modulus (Em) and the critical strain as indicated in Eqs. (10) and (14). Therefore, any changes in the UCS will have effect on the pillar stability. The effect of changing the coefficient of variation (COV) of the UCS on the pillar stability was investigated individually using five different values of COV (10%, 15%, 20% 25% and 30%). The results for the pillar at a depth of 280 m are presented in Fig. 10 to illustrate the effect. The depth 280 m is chosen as its results contain both small and large values of the probability of failure hence it can represent the trend for the remaining depths.

Since COV is a non-dimensional measure of variability the smaller its value is the smaller the amount of uncertainty of the variable. It can be observed from Fig. 10 that the probability of failure increases as the COV of the UCS increases which shows that when the variability in the material properties of the rock mass increased then the probability of pillar failure will also increase. The variation in the COV has no effect on the mean factor of safety (Fig. 11). This shows that when the stability of a pillar is analysed in terms of reliability and probability of failure the conventional factor of safety is an inappropriate indicator of pillar stability. Therefore, in such a case the probability of failure or reliability index is the most useful criterion of pillar stability, because it combines information on both the mean values and the variability. Nevertheless, the factor of safety should not be abandoned in favour of reliability analyses but it has to be used as a complementary tool to reliability analyses (Dai and Wang, 1992).

& 15 -

•s -g

i— Step 0 Step 1 Step 2 Step 3 Step 4 I- Step 5 —►— Step 6

10 15 20 25 30

Coefficient of variation (COV) of intact rock (%)

Fig. 10. Effect of varying COV of intact rock UCS on the pillar probability of failure at different mining steps.

<£ 1.64

B 1.44 £ 1.40

S 1.25 1.17

—■— Step 0 —•— Step 1 -A- Step 2 —T— Step 3 —♦- Step 4 —*— Step 5 ► Step 6

10 15 20 25 30

Coefficient of variation (COV) of intact rock UCS (%)

Fig. 11. Effect of varying COV of intact rock UCS on the mean factor of safety at different mining steps.

6.4. Assessment of pillar performance

The performance of a structure can be described in terms of the specific limit state being reached. The limit state is the acceptance

criteria established by considering the nature, importance, and consequences of failure of the structure under consideration. Many performance criteria in terms of reliability index (b), probability of failure (Pf) and mean factor of safety (mFS) have been

220 230 240 250 260 270 280 Mining depth (m)

■«— Step 0 -•- Step 1 Step 2 — Step 3 -»- Step 4 <— Step 5 Step 6

220 230 240 250 260 270 Mining depth (m)

Fig. 9. Probability of failure (a) and reliability index (b) for different mining depths at different mining steps.

Table 7

Typical values for reliability index and probability of failure and corresponding performance level.

Reliability index 4 Probability of failure Performance level

5.0 2.9E-7 High

4.0 3.2E-5 Good

3.0 1.35E-3 Above average

2.5 6.21E-3 Below average

2.0 0.023 Poor

1.5 0.067 Unsatisfactory

1.0 0.159 Hazardous

proposed for geotechnical structures (Lunder and Pakalnis, 1997; US Army Corps of Engineer, 1997; Priest and Brown, 1983). The pillar performance for the different mining steps at different depth is hereby assessed based on the criterion proposed by US Army Corps of Engineer (1997). The criterion is summarized in Table 7. Fig. 12a-g show b and the corresponding performance levels for different depths and each mining step. The performance levels based on the US Army Corps of Engineer criterion are indicated in Fig. 12a-g. Targeted performance level of any geotechnical project depends on the nature of the project and the level of risk that can be accepted. Risk is a product of the probability of failure and the consequence of the failure. For instance an underground nuclear waste disposal will require high performance level because of the severity of the consequence of its failure. Priest and Brown (1983) proposed b = 2.3 for geotechnical structures if the consequences of failure are moderately serious. Hence for this study a performance level below average was chosen i.e. b = 2.5. Therefore, based on the chosen performance level for the pillar it can be seen from Fig. 12a-g that the targeted performance level, below average, can be achieved for the pillar at all depths at mining steps 0-2 (i.e. W/H = 1.61-1.37). The same performance level can also be achieved for mining depths 205, 220, and 240 m at mining step 4 (i.e. W/H = 1.20). At mining depth 205 m the targeted performance level can only be achieved at mining step 5. The performance level cannot be achieved for the pillar for all the depths at mining step 6.

0 1 2 3 4 5 6

Mining step

Fig. 13. Influence of Rf value on the pillar reliability index and performance level.

It is obvious that the optimum dimension of the pillar that can be achieved depends on the level of the risk that could be tolerated. Therefore, when a targeted performance level is chosen for an underground structure the procedure described in this paper can be used for optimal design of pillar dimension especially when the variability and uncertainty of the design input parameters are to be considered. However, when using this approach with the strain based criterion it should be noted that the value of Rf (failure strength parameter) has a significant influence on the probability of failure as well as the reliability index as illustrated in Fig. 13. Fig. 13 shows the pillar reliability index at mining depth 280 m for different Rf values. For Rf = 0.1 the performance level (i.e. below average) can only be achieved at mining step 0 while the same performance level can be achieved up to mining step 4 if Rf value is 0.3. Hence, the Rf value should be carefully determined through field monitoring programs for specific mine before using this approach.

(a) step 0

(b) step 1

73"" £

205 220 240 260 280 Depth (m)

(e) step 4

205 220 240 260 280 Depth (m)

205 220 240 260 280 Depth (m)

(f) step 5

-•3--4

205 220 240 260 280 Depth (m)

(c) step 2

(d) step 3

-ES2-.

"""id"'

...Good...

. Above average..

Below average

205 220 240 260 Depth (m)

(g) step 6

220 240 260 Depth (m)

V/A 280

Above average Below average

205 220 240 260 280 Depth (m)

Fig. 12. Pillar reliability index (b) and the corresponding performance level for different depths at each mining step.

7. Conclusion

A stochastic assessment of the pillar stability of the Laisvall mine has been presented in this paper. The uncertainty and variability in the rock mass parameters as well as that of the horizontal in-situ stresses were considered. Elastic 3D numerical analyses were performed to determine the axial strain in the pillar at different mining steps for different depths. Artificial Neural Network (ANN) was used to generate a closed-form expression for the relationship between the random variables (i.e. Em and OH) and the axial strain using the results of the 3D numerical analyses. The closed-form equation together with the strain at failure was used to evaluate the performance function. The performance function was determined from Monte Carlo simulations of 105 trials and the probability of failure and reliability index were determined from the simulation results. Based on the results and discussions presented in this study, the main conclusions are summarized as follows:

1. Rock mass properties as well as in-situ stresses are intrinsically variable such that using their mean values could have significant impact on the design performance as shown in this study. Therefore, this variability should be properly considered in the determination of the input parameters for the numerical analyses using a stochastic approach. This will enable the evaluation of the probability of failure during the planning stage, thereby reducing the risk to an acceptable level.

2. The strain-based criterion used in this study seems appropriate when the variability in the rock mass properties is considered in the analysis as pillars with different rock mass properties deform differently when subjected to the same stress condition. The strain-based criterion can also be used without having to determine the dimension of the pillar which is required to determine the strength of the pillar from empirical equation when using stress/strength criterion. Furthermore, strain is relatively cheaper and easier to measure in the field than the stress which is expensive and time consuming. Therefore using the approach presented in this study may offer some potential and advantages than the traditional strength/stress based criterion.

3. The trained ANN model is a good approximation tool to generate an implicit relationship between the random variables (Em and OH) and the axial strain as shown in this study, thereby reducing the number of FLAC3D numerical simulations. Also the ANN model was able to capture the tails of the distribution of the random variable because, the range of the value for the random variable which covers almost 99.7% (i.e. ±3 standard deviation from the mean) of the distribution was considered. It should be noted that the trained ANN model can only be used in the similar situation described in this study within the range of values of the training dataset.

4. The width-to-height ratio of a pillar is very significant to the stability of the pillar and as the ratio reduces through the mining steps the reliability index reduces and consequently the probability of failure of the pillar increases.

5. Increase in the in-situ stresses which could be linked to the increase in the depth affects the stability of the pillar. As the depth increases the probability of pillar failure increases and the reliability index of the pillar reduces.

6. The variation in the COV of the intact rock UCS has effect on the reliability index and the probability of failure. This study indicates that, an increase in the COV of the intact rock UCS results in an increase in the probability of failure and reduces the reliability index. However, the variation in the COV has no effect on the mean factor of safety. Therefore, the mean factor of safety

may not be an adequate parameter to assess the stability of pillar when considering the variability and uncertainties in the rock mass material properties. Instead, the reliability index and probability of failure are the most useful indicators for assessing the stability of pillar in such situation. 7. The assessment of the pillar performance depends on the nature, importance and the consequence of the failure. The management decision on the level of risk that can be taken may affect the design of the pillar with regard to its dimension during the planning stage.

The approach presented in this study shows that when considering the uncertainty and variability of rock mass properties, the stochastic approach can be a useful tool for analysing and assessing the pillar stability. This will allow setting acceptable criteria for pillar design based on the minimum acceptable risk or performance level.

Acknowledgements

The authors acknowledged the Centre of Advanced Mining and Metallurgy (CAMM) at the Lulea University of Technology, Lulea, Sweden for financial support. The authors thank Mr. Per-Ivar Marklund of Boliden Mineral AB for providing information about the Laisvall mine.

References

Alber, M., Heiland, J., 2001. Investigation of a limestone pillar failure. Part 2: Stress history and application of fracture mechanics approach. Rock Meck. Rock Engng. 34, 187-199.

Aldrich, C., Reuter, M.A., 1994. The application of neural nets in the metallurgical

industry. Miner. Eng. 7, 793-809. Barton, N.R., Bandis, S.C., 1990. Review of predictive capability of JRC-JCS model in engineering practice. In: Proceedings of International Symposium on Rock Joints, Balkema, Rotterdam, pp. 603-610. Barton, N.R., Lien, R., Lunde, J., 1974. Engineering classification of rock masses for

the design of tunnel support. Rock Mech. 6, 189-239. Beale, M.H., Hagan, M.T., Demuth, H.B., 2012. Neural Network Toolbox User's Guide.

The MathWorks Inc., Natick, MA. Cai, M., 2011. Rock mass characterization and rock property variability considerations for tunnel and caven design. Rock Mech. Rock Eng. 44, 379-399. Cai, M., Kaiser, P.K., 2006. Visualization of rock mass classification systems. Geotech.

Geol. Eng. 24, 1089-1102. Cai, M., Kaiser, P.K., Tasaka, Y., Maejima, T., Morioka, H., Minami, M., 2004a. Generalized crack initiation and crack damage stress thresholds of brittle rock masses near underground excavations. Int. J. Rock Mech. Min. Sci. 41,833-847. Cai, M., Kaiser, P.K., Uno, H., Tasaka, Y., Minami, M., 2004b. Estimation of rock mass strength and deformation modulus of jointed hard rock masses using the GSI system. Int. J. Rock Mech. Min. Sci. 41, 3-19. Cai, M., Kaiser, P.K., Tasaka, Y., Minami, M., 2007. Determination of residual strength parameters of jointed rock masses using GSI system. Int. J. Rock Mech. Min. Sci. 44, 247-265.

Cauvin, M., Verdel, T., Romuald, S., 2009. Modelling uncertainties in mining pillar

stability analysis. Risk Anal. 29,1371-1380. Chen, G., Jia, Z., Ke, J., 1997. Probabilistic analysis of underground excavation

stability. Int. J. Rock Mech. Min. Sci. 34, 51.e1-51.e16. Coates, D.F., 1981. Rock mechanics principles. CANMET, energy, mines and

resources Canada. Monograph 874. Dai, S.H., Wang, M.O., 1992. Reliability Analysis in Engineering Applications. Van

Nostrand Reinhold, New York. Deere, D.U., 1968. Geological consideration, rock mechanics in Engineering Practice, Division of Civil Engineering, School of Engineering, University of Wales, Swansea. Wiley, New York, pp. 1-20. Deng, J., Yue, Z.Q., Tham, L.G., Zhu, H.H., 2003. Pillar design by combining finite element methods, neural networks and reliability: a case study of Feng Huangshan copper mine, China. Int. J. Rock Mech. Min. Sci. 40, 585-599. Dohyun, P., Hyung-Mok, K., Dong-Woo, R., Won-Kyong, S., Choon, S., 2012. Application of a point estimate method to the probability limit-state design of underground structures. Int. J. Rock Mech. Min. Sci. 51, 97-104. Exadaktylos, G., Stavropoulou, M., 2008. A specific upscaling theory of rock mass parameters exhibiting spatial variability: analytical relations and computational scheme. Int. J. Rock Mech. Min. Sci. 45,1102-1125. Faraway, J., Chatfield, C., 1998. Time series forecasting with neural network: a

comparative study using the airline data. Appl. Stat. 47, 231-250. Goh, A.T.C., Kulhawy, F.H., 2003. Neural network approach to model the limit state surface for reliability analysis. Can. Geotech. J. 40,1235-1244.

Griffiths, D.V., Fenton, G.A., Lemons, C.B., 2002. Probabilistic analysis of underground pillar stability. Int. J. Numer. Anal. Methods Geomech. 26, 775791.

Hoek, E., Brown, E.T., 1997. Practical estimates of rock mass strength. Int. J. Rock Mech. Min. Sci. 34,1165-1186.

Hoek, E., Diederichs, M.S., 2006. Empirical estimation of rock mass modulus. Int. J. Rock Mech. Min. Sci. 43, 203-215.

Hoek, E., Kaiser, P.K., Bawden, W.F., 1995. Support of Underground Excavations in Hard Rock. A.A. Balkema, Rotterdam Brookfield.

Hoek, E., Carranza-Torres, C., Corkum, B., 2002. Hoek-Brown failure criterion - 2002 ed. In: Proceedings of the 5th North American Rock Mechanics Symposium and 17th Tunnelling Association of Canada Conference: NARMS-TAC 2002, July 710, University of Toronto, pp. 267-271. <www.rocscience.com> (2.10.02).

Hornik, K., Stinchcombe, M., White, H., 1989. Multilayer feed-forward networks are universal approximators. Neural Networks 2, 359-366.

Idris, M.A., Saiang, D., Nordlund, E., 2011. Probabilistic analysis of open stope stability using numerical modelling. Int. J. Min. Mineral Eng. 3, 194-219.

Itasca, 2012. FLAC3D - Fast Langrangian Analysis of Continua in Three - Dimension, Version 5.0. <www.itascacg.com>.

Joughin, W.C., Swart, A.H., Wesselo, J., 2000. Risk based chromite pillar design: Part II- Non-linear modelling. In: Proceedings SANIRE Symposium.

Kanellopoulas, I., Wilkinson, G.G., 1997. Strategies and best practice for neutral network image classification. Int. J. Remote Sensing 18, 711-725.

Kim, K., Gao, H., 1995. Probabilistic approaches to estimating variation in the mechanical properties of rock masses. Int. J. Rock Mech. Min. Sci. Geomech. Abst. 32, 111-120.

Krauland, N., Soder, P.-E., 1987. Determining pillar strength from pillar failure observation. Eng. Min. J. 8, 34-40.

Li, H.Z., Low, B.K., 2010. Reliability analysis of circular tunnel under hydrostatic stress field. Comput. Geotech. 37, 50-58.

Lilly, P.A., Li, J., 2000. Estimating excavation reliability from displacement modelling. Int. J. Rock Mech. Min. Sci. 37, 1261-1265.

Lunder, P.J., Pakalnis, R., 1997. Determination of the strength of hard-rock mine pillars. Bull. Can. Min. Metall. 90, 51-55.

Marinos, P., Hoek, E., 2000. GSI: a geologically friendly tool for rock mass strength estimation. In: Proceedings of the GeoEng2000 at the International Conference on Geotechnical and Geological Engineering, Melbourne, Technomic Publishers, Lancaster, pp. 1422-1446.

Marklund, P.I., 2013. Personal Communication.

Martin, C.D., Maybee, W.G., 2000. The strength of hard-rock pillars. Int. J. Rock Mech. Min. Sci. 37, 1239-1246.

MATLAB, 2012. Version 2012b. The MathWorks Inc. <www.mathworks.com>.

Najafi, M., Jalali, S.E., Bafghi, A.R., Sereshki, F., 2011. Prediction of the confidence interval for stability analysis of chain pillars in coal mines. Saf. Sci. 49,651-657.

Obert, L., Duvall, W.I., 1967. Rock Mechanics and the Design of Structures in Rock. John Wiley & Sons Inc., New York.

Oreste, P.A., 2005. Probabilistic design approach for tunnel supports. Comput. Geotech. 32, 520-534.

Palisade Corporation, 2001. @RISK version 4.0. <www.palisade.com>.

Palmstrom, A., 1995. RMi - a rock mass characterization system for rock engineering purposes. Ph.D. Thesis, University of Oslo, Norway.

Pine, RJ., 1992. Risk analysis design applications in mining geomechanics. Trans. Inst. Min. Metall, Sect. A, A149-A158.

Priest, S.D., Brown, E.T., 1983. Probabilistic stability analysis of variable stopes. Trans. IMM. Sect. A 92, A1-A12.

Recio-Gordo, D., Jimenez, R., 2012. A probabilistic extention to the empirical ALPS and ARMPS systems for coal pillar design. Int. J. Rock Mech. Min. Sci. 52,181187.

Rethati, L., 1988. Probabilistic solutions in geotechnics. Developments in Geotechnical Engineering, vol. 46. Elsevier, Amsterdam.

Rosenblueth, E., 1981. Two-point estimates in probabilities. Appl. Math. Modell. 5, 329-335.

Rumelhart, D.E., Hinton, G.E., Williams, R.J., 1986. Learning internal representation by error propagation. In: Rumelhart, D.E., McClelland, J.L. (Eds.), Parallel Distributed Processing. MIT press, Cambridge, pp. 318-362.

Sahin, M.A., Jaksa, M.B., Maier, H.R., 2001. Artificial neural network applications in geotechnical engineering. Aust. Geomech. 36, 49-62.

Sakurai, S., 1981. Direct strain evaluation technique in construction of underground opening. In: Proceedings of 22th US Symp. Rock Mech, MIT, Cambridge, pp. 278-282.

Sakurai, S., 1997. Strength parameters of rocks determined from back analysis of measured displacements. In: Proceedings of the First Asian Rock Mechanics Symp. ISRM, Seoul, pp. 95-99.

Sakurai, S., 1999. Interpretation of field measurements in tunnel practice. In: Proceedings of the Ninth International Congress of Rock Mechanics, Paris, pp. i-vii.

Schweiger, H.F., Thurner, R., Pöttler, R., 2001. Reliability analysis in geotechnics with deterministic finite elements. Int. J. Geomech. 1, 389-413.

Seibi, A., Al-Alawi, S.M., 1997. Prediction of fracture toughness using artificial neural networks (ANNs). Eng. Fract. Mech. 56, 311-319.

Sjöberg, J., 1992. Failure modes and pillar behaviour in the Zinkgruvan mine. In: Proc. 33rd U.S. Rock Mech. Symp (Santa Fe, June 8-10,1992). Rotterdam: A.A. Balkema, pp. 491-500.

Söder, P.E., Krauland, N., 1990. Determination of pillar strength by full scale tests in the Laisvall mine. In: Proceedings of 11th Plenary Scientific Session of the International Bureau Strata Mechanics, World Mining Congress, Novosibirsk, pp. 39-59.

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, 224-235.

Tawadrous, A.S., DeGagne, D., Pierce, M., Mas Ivars, D., 2009. Prediction of uniaxial compression PFC3D model micro-properties using artificial neural networks. Int. J. Numer. Anal. Meth. Geomech. 33, 1953-1962.

Torano, J., Rodriguez, R., Ramirez-Oyanguren, P., 2000. Probabilistic analysis of subsidence-induced strains at the surface above steep steam mining. Int. J. Rock Mech. Min. Sci. 37, 1161-1167.

U.S. Army Corps of Engineers, 1997. Engineering and design introduction to probability and reliability methods for use in geotechnical engineering. Engineering Technical Letter No. 1110-2-547, Department of the Army, Washington, DC.

Wattimena, R.K., Kramadibrata, S., Sidi, I.D., Aziz, M.A., 2013. Developing coal pillar stability chart using logistic regression. Int. J. Rock Mech. Min. Sci. 58, 55-60.