Application of two-level factorial design to sensitivity analysis of keyblock statistics from fracture geometry

International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

Application of two-level factorial design to sensitivity analysis of keyblock statistics from fracture geometry P. Starzeca,*, J. Anderssonb a

Swedish Geotechnical Institute,Chalmers Vasa, Hugo Grauers gata 5B, SE-412 96 Gothenburg, Sweden b JA Streamflow, Stockholm, Sweden Accepted 5 March 2002

Abstract A stochastically modeled fracture network offers potential for more realistic assessment of stability status in underground excavations than predictions based entirely on deterministic features. The reliability of probabilistic models, however, depends strongly on an accurate estimation of the model’s variables, i.e., the fracture network properties from the field and laboratory observations. In this study, predictions of keyblocks by implementing stochastically generated fractures in the Central Storage Facility for Spent Nuclear Fuel (CLAB 2 Centralt Lager Anv.ant Br.ansle) located in southeast Sweden are presented. The fracture network model is built by using fracture mapping in the floor of the facility and incorporates fracture size, shape, orientations, termination mode, spatial arrangement and fracture mechanical properties. The predicted volume of individual keyblocks is best-fitted with the Pareto probability distribution function. Subsequently, a statistical two-level factorial analysis is performed to examine the impact of both single fracture properties and their interactions on the predictions made. In the factorial experiments, the block predictions are made for eight different fracture models where three factors: fracture radii, orientation and termination are each assigned two levels intentionally departing from the best estimates found for the CLAB 2 site. This allows us to express the experiment results as the degree to which each of the eight computed block statistics deviated from the most likely prediction. It is found that fracture orientation is the only statistically significant factor influencing the keyblock statistics while the input from other variables/fracture properties and their interactions is less significant. The results of our study yield a prospective approach for improving the effectiveness of the stochastic model variable estimation and for more optimal field mapping strategies. r 2002 Elsevier Science Ltd. All rights reserved.

1. Introduction The volume, shape and amount of unstable rock blocks formed by intersections between joints/ fractures and contours of underground chambers/ tunnels depend on both the dimensions and the geometry of the excavation itself as well as on the geometry and other properties of fractures/joints intersecting the excavation. The commonly applied Block Theory of Goodman and Shi [1] focuses on the potentially largest keyblocks and orientations of major joint systems in relation to the orientation of the underground object. Several authors [2–5] attempted to employ a probabilistic approach for *Corresponding author. E-mail address: [email protected] (P. Starzec).

keyblock predictions based on the stochastic representation of fracture geometries and their locations. Many other researchers addressed the predictions of unstable rock blocks both for underground facilities and for rock slopes, and the numbers of the numerical codes for keyblock predictions and rock support design have been developed; SATIRN [6], UNWEDGE [7,8], DRKBA [9], MSB [10], KBTUNNEL [11], to mention only a few. These codes are all based on similar principles but differ in terms of the type of input data, model assumptions and potential outcome from the simulation analysis. Other well-known codes/methodologies offer less potential for modeling more complex fracture networks but on the other hand permit the examination of interactions between block triggering, displacement, deformation and stress field in relation to the time factor: for example, 3DEC [12] and DDA [13].

1365-1609/02/$ - see front matter r 2002 Elsevier Science Ltd. All rights reserved. PII: S 1 3 6 5 - 1 6 0 9 ( 0 2 ) 0 0 0 2 6 - 6

244

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

One important aspect of keyblock analysis refers to the estimation process of fracture geometries and mechanical properties based on field and/or laboratory measurements. Depending on the type of estimated fracture variable/property, its uncertainty will have a different impact on the final block predictions. Further, in some circumstances the interactions between several variables may influence the predictions to a higher degree than the individual variables do. Knowing which fracture properties significantly affect the amount of unstable blocks, their size or total removable volume along an underground chamber and which properties are of less importance could result in the optimization of field sampling strategy and of more cost-effective tunnel support design. The majority of these kinds of studies have been focused on how predicted block statistics depends on varying one property at the time e.g., by varying fracture size, fracture intensity or tunnel dimensions [14–16]. Therefore, a deeper analysis of the possible interactions among properties of the fracture network and their quantitative bearing on block predictions could help us make more reliable geological risk estimates. This paper presents: (i) probabilistic three-dimensional keyblock predictions based on two-dimensional fracture mapping from the floor of the CLAB 2 underground chamber hosted in crystalline rock unit, southeast Sweden, and (ii) sensitivity analysis of keyblock predictions to fracture size, orientations and terminations based on statistical two-level factorial design.

2. Stochastic model of fracture network at the CLAB 2 site 2.1. Fracture mapping data at the CLAB 2 site The probabilistic keyblock predictions imply that simulation of potentially unstable rock blocks relies on a stochastic fracture field generated from observed/measured fracture properties. In the current work, the FracMan fracture generator was applied [17] allowing us to represent such properties as fracture orientations, size, shape, termination mode, intensity and spatial arrangement with the best-fit continuous statistical distributions or mathematical/statistical processes estimated from measured quantities at the site. The process of model generation should be preceded by studies on structural/geological homogeneity of the site to determine whether or not fracture properties are spatially independent. Relying on results from such study, the modeled rock portion is represented either by stationary models where fracture properties exhibit the same statistical moments within the entire volume or by nonstationary or multidomain models, i.e., an assembly of more than one stationary model where each model

represents a different tectonic/structural or lithological unit. Although there are no strict procedures designed entirely for homogeneity studies on fracture geometries, a good start is to apply non-parametric equivalents to classical p- and F -statistics [18] for testing mean and variance of a variable in question within different sample populations (different lateral and/or vertical locations at the site) and gradually switch to more advanced statistical procedures. Kulatilake et al. [19] make inferences about the statistical homogeneity of jointed rock mass by studying the spatial variation of box fractal dimension for trace lengths. Dershowitz and Ushida [20], give an overview of several methods for spatial analysis of fracture field data including trend analysis. In a previous study [21], the geological/structural homogeneity at the CLAB site was investigated by using non-parametric Kolmogorov–Smirnov, Mann–Whitney and Kruskal–Wallis tests on fracture intensity obtained from a number of boreholes. Also, the moving average for fracture intensity vs. depth was calculated to account for possible trends in vertical direction. However, the results did not show any clear signs of anisotropic or inhomogeneous fracturing pattern and a one-domain, stationary stochastic model was adopted for the site. The excavated CLAB 2 chamber is intended to serve as an interim storage facility for spent nuclear fuel. It is located about 25 m underground and is 115 m in length, 21 m in width and 27 m in height. Fig. 1 shows fractures and weakness zones mapped on the floor of the cavern divided in the figure into two sections 57.5 m long each (Fig. 1a depicts the northern half of the floor and Fig. 1b the southern half). While in the previous work on keyblock predictions for the site the veins and the largeaperture fractures were selected for fracture network generation [21], in this study the access to the mapped fracture trace data made it possible to apply a more geomechanically reliable criterion for the selection of the proper type of discontinuities for setting up a stochastic fracture model. According to geological expertise at the site [22], the open fractures and the fractures filled with chlorite were the most significant for keyblock prediction due to their low shear strength compared to tight fractures or to those filled with quartz. For the purpose of stochastic fracture network modeling, we selected 411 fractures as a representative hard data-input for variable estimation. The shortest trace length selected was 1.1 m long and the longest trace was 25 m long, but longer traces were also expected since a certain number of the mapped traces were censored by the floor’s boundaries. The approach adopted in the FracMan code is that the observed statistics for fracture geometric properties obtained from 1D and/or 2D fracture mapping are approximated with theoretical probability density functions and, consequently, each generated fracture in the

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

245

Fig. 1. Fracture traces (light line) and weakness zones (chain-patterned lines) with depicted dip angle mapped in the CLAB-2 chamber floor. The upper figure represents northern part of the floor and the lower figure depicts the southern part. Each part of the floor is 57.5 m long and 21 m wide.

model is a product of one Monte Carlo sampling from each probability function representing a certain fracture property. Hence, before the stochastic fracture model reflecting the closest-to-nature conditions for the CLAB 2 was set up, the all modeled fracture variables (orientation, size, shape, spatial pattern, termination, intensity) had to be estimated. The remainder of this chapter presents the estimation methodology for frac-

ture variables from the observed field data and their match to approximate continuous theoretical statistical probability functions. 2.2. Discontinuity orientation In order to account for natural variation in discontinuity orientations a statistical model is built that

246

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

represents the discontinuity orientation characteristics of the entire rock mass studied. Usually, hemispherical projection methods are utilized as a graphical approach to the analysis of orientation data. Sets of preferred orientation can be identified by eye. The orientation limits for each set can then be specified either in terms of a range of trend and plunge angles or by the manual delineation of a range of orientations. This method has the advantage of allowing an individual’s expertise and familiarity with a particular site to play a part in the identification of sets but has the disadvantage of being subjective and therefore susceptible to personal bias and inconsistency. A less subjective approach is based on statistical cluster identification algorithms with the assumption that a discontinuity set will exhibit a significantly greater degree of clustering than a totally random distribution of orientations [23,24]. For the CLAB 2 site, the mapped fractures were first corrected for orientation sampling bias with the upper limit for the correction factor CF ¼ 6: The applied correction compensated for subhorizontal fractures having low probability to be mapped on the tunnel floor as they run subparallel to the sampling plane. The data were corrected by multiplying each fracture by the integer of the reciprocal of sin a; where a is the fracture dip angle (the floor of CLAB 2 is a horizontal plane). The upper limit of the correction factor was set up to avert increasing to infinity the amount of low-dipping fractures as sin a approaches 0 for fractures dipping at very low angles in respect to the tunnel floor. Establishing the upper limit for CF is usually not explicit and in other studies, it has been reported to be between 5 and 10 depending on the author’s choice [25,26]. In this study, the upper limit for correction factor was set up to CF ¼ 6: To limit a possible personal misjudgment in terms of orientation cluster recognition, we initially computed eigenvalues and spherical variance [18] on the data corrected for the orientation sampling bias, i.e., for 564 fractures as represented by poles to fracture on the lower hemisphere projection in Fig. 2. The eigenvalue analysis showed that l1 > l2 > l3 (where l1 is the 1st eigenvalue, l2 is the second and l3 is the third), which means that the data are not evenly distributed over sphere and are scattered around great circle but with preferred direction [27]. Also, the ratio found between the logarithm of the first and second and the logarithm between the second and third eigenvalue were both close to 1.0 indicates that the orientations of pole vectors are not randomly distributed on the sphere [28]. Hence we could reject the hypothesis of uniformity on the sphere; and also supported by visual inspection, we concluded that there was at least one dominant fracture orientation set, which could be specified by the eigenvector associated with the first eigenvalue, i.e., pole

Fig. 2. Lower hemisphere projection of poles to 564 fractures mapped and corrected in the CLAB 2 floor. Fractures corrected for orientation sampling bias with maximum correction factor CF ¼ 6:

azimuth about 3371 and pole plunge about 301(cf. Fig. 2). The rejection of the uniformity by statistical methods is roughly confirmed by the geological interpretation from the tectonic record of the site. The area ( was formed through magmatic intrusions (Smaland and Varmland intrusions) following after Svekokarelian . orogenes (1800 million years ago). The region has been subjected to several contracting tectonic movements of different magnitudes. It is expected that such tectonic ‘‘diversity’’ would have resulted in rather non-uniform distribution of fracture orientations and formation of distinguishable fracture sets. Next, the interactive set identification system (ISIS) developed by Dershowitz and his co-workers [17] was applied in order to separate fractures into distinct orientation clusters and to find the best statistical fit between theoretical probability distributions and the observed data. We examined several different clusters’ configurations and several different probability distributions to find the best-fit between the observed data and continuous statistical distributions. Finally, the best-fit to the observed data was obtained after separating the fractures into two fracture sets and approximating them with Fisher’s probability density function [18,29]. However, the Kolmogorov–Smirnov goodness-of-fit statistics for this pdf was rather poor and not statistically significant. Table 1 summarizes the output from the ISIS analysis for the orientation data in Fig. 2. The concentration coefficient k in Table 1 is a measure of the degree of clustering, or preferred orientation within the

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

247

Table 1 Summary of orientation set identification for fractures selected from floor mapping in the CLAB 2 and corrected for orientation sampling bias Orientation set

Mean pole azimuth and plunge (1)

Concentration coefficient for best-fit probability density function, k

Fracture percentage of the observed total (%)

Set 1 Set 2

336.6 247.6

7.63 9.81

66.1 33.9

29.9 14.5

population [29]. High values of k imply less scatter around the resultant orientation vector. For more details, see Refs. [18,25] where the authors demonstrate a simple way to calculate k for fractures and other directional geological data types.

Table 2 The most probable fracture radius distribution obtained from simulated sampling with FracSize procedure for fractures mapped in the CLAB 2 floor Orientation set

Distributional form

Distribution parameters

2.3. Fracture size

Set 1 (336.6, 29.9) Set 2 (247.6, 14.5)

Lognormal Lognormal

Mean=4.8 m, std. dev.=1.7 m Mean=4.5 m, std.dev.=1.4 m

In the current work, fracture shapes were approximated with ellipsoidal discs with aspect ratio, i.e., the ratio between longer and shorter disc axis 2:1. The decision to assume that fractures mapped in the CLAB 2 floor were elongated rather than equidimensional was based on the observations of fracture traces in the tunnel [22]. Even when trace length data are accessible, the derivation of the discontinuity size is still an intricate task. Priest [25] presented an overview of the most common difficulties in a proper estimation of fracture dimensions with focus on an analytical approach originating from the work done by Warburton [30]. Another method was proposed by La Pointe et al. [31] who advocated the use of fracture size derivation through simulated sampling using forward stochastic modeling. In order to estimate the representative fracture radius for the CLAB 2 site we used FracSize module implemented in the FracMan fracture simulator. The FracSize module determines the distribution of fracture radius that best matches the observed trace length data. In brief, the procedure starts with generating a fracture field with a user-assumed fracture radius distribution. In the next step, the FracSize simulates the sampling process on a synthetic plane or a borehole and compares the sampled trace length statistics with the observed one. To automate the search of radius distribution, the optimization algorithm based on simulated annealing is applied [32], and the coherence between the simulated and observed trace lengths are measured by Kolmogorov– Smirnov or Chi-Squared statistics. The results from the fracture radii estimation for the CLAB 2 data are given in Table 2. 2.4. Fracture locations in space Knowledge about the spatial fracture pattern as inferred from either one-dimensional borehole sampling

or two-dimensional outcrops/tunnel mapping is crucial for setting up a realistic three-dimensional stochastic spatial fracture pattern. Fracture organization in space has been treated by many authors [33–36] and is still a challenging task for structural geologists and modelers. Although a majority of field studies have shown that fractures tend to be randomly distributed in space [25], some workers have provided evidence for some more systematic fracture location patterns that could be described by processes other than pure Poisson [35,37,38]. Previous studies on fracture spacing as found from the borehole camera and core mapping at the CLAB site by applying power density spectra and variogram analysis did not show any clear tendency for fracture clustering, and a stochastic fracture model based on a random location pattern was adopted as representative for the site [21]. Likewise, in the current study, the Baecher algorithm revised terminations (BART) fracture location model was implemented [17] where fractures were located uniformly in space, following the Poisson process, while fracture termination modes, resulting in non-uniform fracture locations, were allowed. Fracture termination was assigned by termination fraction T (Eq. (1)). T¼

tf ; tf þ trm

ð1Þ

where tf is the number of terminations against another fracture and trm the number of terminations in intact rock mass. For fractures mapped in the CLAB 2 floor, the termination fraction was found to be T ¼ 13%: In the BART model, at first ð1  TÞ * N fractures were generated based upon uniformly located fracture surface

248

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

points (N is the total number of fractures). T fractures terminating at intersections were then generated from locations uniformly positioned on the surface of existing fractures. These positions were not the fracture centers, but rather the locations at which the fracture termination occurred. The fractures generated from these locations were defined with centers at a distance away from the termination location such that the fracture surface area distribution for the generated fracture was preserved. Even if the BART model proved to correctly reflect the fracture pattern for the CLAB 2 as the model was verified through synthetic trace sampling (see next section), this was not a test of clustering evidence and a more detailed investigation on possible fracture clustering by means of 2D fractal box counting method might be done in further studies at the site.

2.5. Fracture intensity In the FracMan code, fracture intensity is expressed either as the total number of fractures in a given volume or as the volumetric fracture intensity, which is the ratio between the total area of all fractures and the volume of the model. Since the volumetric fracture intensity and the two-dimensional intensity measure (ratio between the total fracture trace lengths and the area within which fractures were mapped) are linearly related [39], it is usually easier to verify a stochastic model based on the volumetric fracture intensity rather than based on the total fracture number. A simple reason for adopting volumetric intensity is that after the model has been generated a synthetic fracture trace mapping on a plane reflecting dimensions and orientation of real rock exposure/outcrop can be done; and, accordingly, statistics for synthetically sampled data and the traces observed in the field can be compared. Subsequently, the originally assumed volumetric intensity is adjusted until the smallest allowed discrepancy between the twodimensional synthetic and the observed fracture intensity is reached. The observed two-dimensional fracture intensity in the CLAB 2 floor was found to be 1.32 m/m2. The stochastic fracture network with the BART model was generated within a slab 130 m long, 60 m wide and 60 m high. For 20 Monte Carlo realizations with fracture geometries represented by probability density distributions as described in Sections 2.1–2.3 (see also Table 2), the synthetic sampling of fracture traces was done and compared to the observed data to estimate the most probable volumetric fracture intensity. Finally, the best estimate was obtained for the volumetric intensity of 1.55 m2/m3.

3. Keyblock predictions 3.1. Methodology Computation of the number of keyblocks along the CLAB 2 facility having a horseshoe shape (115 m long, 21 m wide and 27 m high), and with the long axes striking N12E was accomplished with RockBlock PC-based code [40]. The main advantage of using RockBlock was its full compatibility with FracMan output fracture models. Keyblock predictions were made for 10 Monte Carlo runs of the stochastic fracture network with its estimated properties as depicted in Table 3. The horizontal stress vector was not considered, thus the keyblock predictions here were regarded as the worst possible scenario where no forces other than friction on the fracture planes have a stabilizing effect on block fall/sliding. The mechanical properties of fractures at the site estimated in earlier studies [41] were assumed to be also valid for the type of discontinuities considered in this work. Factor of safety SF was based on limit equilibrium assumptions and the Barton–Bandis model [42]. Free falling blocks had FS ¼ 0; stable blocks had an effectively infinite FS; while for one-plane and twoplane sliding the factor of safety was calculated according to Eqs. (2) and (3), respectively: FS ¼

jN1 j tanðfÞ ; S

ð2Þ



jN1 j tan f1 þ jN2 j tan f2 FS ¼ ; S12

ð3Þ

Table 3 Components of the stochastic fracture model built from fracture mapping in the CLAB 2 floor Model variable

Properties

Type of discontinuities

Open fractures and fractures filled with chlorite BART model (modified random location model) Two sets: Set 1: 336 (pole trend), 29.9 (plunge), Set 2: 247 (pole trend), 14.7 (plunge), Fisher’s pdf Ten-sided polygons, aspect ratio 2:1 Set 1: radii log N distributed, mean 4.8 m, std. dev. 1.7 m Set 2: radii log N distributed, mean 4.5 m, std. dev. 1.4 m 1.55 m2/m3

Fracture location pattern Fracture orientations

Fracture shapes Fracture size

Volumetric fracture intensity Fracture terminations Model size Friction angle Safety factor

13% 130 * 60 * 60 m3 301 (applies to both open and filled fractures) Barton–Bandis model

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

where N is the normal force to face, f the friction angle and S the magnitude of shear force. The additional subscripts in Eq. (3) were introduced to account for twoblock face sliding. 3.2. Results Any probabilistic prediction of keyblock statistics raises the usual question: What keyblock properties and which distributional parameters do constitute a reliable basis for rock support design? As pointed out by Chan and Goodman [43], the number of keyblocks alone is not a good indicator of the support required, but could yield useful information about the average size of keyblocks in combination with the other measures. While the matter of the correct statistical representation of keyblock predictions is still a topic of scientific dispute in many academic and professional communities, the authors suggest what they consider the most

critical keyblock statistical parameters and present them in Table 4. The reported predictions are valid after all blocks below 1000 kg in weight have been arbitrarily eliminated from the analysis since their detrimental impact on tunnel stability was considered of little significance. The lower weight cut-off corresponds to the rock volume of 0.38 m3. We also point out that the keyblock statistics presented in Table 4 consider only free falling and sliding keyblocks with a factor of safety FSo1: Fig. 3 presents the best-fit pdf to the keyblock predictions resulting from combining all falling and sliding blocks with FSo1 found after completing all 10 Monte Carlo simulations. The Pareto pdf (Fig. 3) is characterized by two parameters: location L; which in this particular case physically corresponds to the minimum block volume, and shape b; which as the name implies describes the curvature of the probability function. The formula for Pareto pdf computation implemented in the Crystal Ball software package [44] used for matching our data with continuous pdfs is defined as follows: f ðxÞ ¼

Table 4 Statistical parameters of the predicted keyblocks for the CLAB 2 facility. Only free falling blocks and sliding blocks with FSo1 included. The predictions are based on 10 Monte Carlo runs Distributional parameter

Predictions

Median keyblock volume Vm

0.8 m3

Average of the largest keyblock volume Vmax and standard deviation db

Vmax ¼ 12 m3 db ¼ 10:2 m3

Average of total unstable volume Vt and standard deviation dt

Vt ¼ 57 m3 dt ¼ 15 m3

Best-fit pdf to block volume

Pareto distribution: location L ¼ 0:39 m3, shape b ¼ 0:96

Average number of keyblocks Nb and standard deviation db

Nb ¼ 38 db ¼ 7

249

bLb ; xðbþ1Þ

ð4Þ

where x is the block volume here.

4. Sensitivity analysis for keyblock predictions As mentioned in the introductory chapter, several studies have been reported that focus on the relationship between keyblock predictions, fracture geometry and excavation shape, size, and orientation [14–16]. Most of this work, however, addresses the impact of one fracture variable on the prediction at the time and does not investigate possible interactions between variables and how such interactions might affect the predictions made. In the following sections, we present a brief description of factorial analysis [45] and demonstrate how this methodology is applied to study the interrelation among

Fig. 3. The simulated keyblock volume (vertical bars) in m3 and the best-fit Pareto probability density function (solid line) found for the 115 m long CLAB 2 tunnel. Only free falling blocks and sliding blocks with FSo1 with volume larger than 0.38 m3 (B1000 kg) included.

250

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

fracture geometric properties and its impact on keyblock predictions for the CLAB 2 facility. 4.1. Basic principles of factorial analysis Factorial analysis or, more precisely, factorial designs gained its recognition in the late 1960s, and the concept was first applied for cost-effective optimization of industrial processes [46], increasing efficiency and productivity [47], and improving computer performance [48]. The method offers a relatively simple way to study the impact of single variables constituting a mathematical or physical model for a given phenomenon on the model performance as well as the impact from the interactions among several variables. To perform a general factorial design, an investigator selects a number of levels for each variable (factor) and then runs experiments with all possible combinations of levels and variables. If there are l levels for each variable and there are n variables, the complete arrangement of all experimental runs is called an l n factorial design. Table 5 illustrates a design matrix for 22 factorial experiment involving 22=4 experiments in which there are two variables A and B and the response from the experiments is a measurable quantity Y : In the experiment, the variables A and B; which can be both quantitative and qualitative are assigned two different values/levels coded as plus and minus signs. The notation AB in Table 5 represents the interaction between the main factors A and B and its sign is the product of their sign multiplication. Since the experimental results y1y: y4 are valid for different configurations of variables-levels, the next step of the procedure is to extract the main effects, i.e., the contribution from the main factors as well as the interactive effect(s) and to quantify them. For instance, the main effects from factors A and B are after Box et al. [45]: MA ¼ ðy1 þ y2  y3 þ y4 Þ=2

ð5Þ

MB ¼ ðy1  y2 þ y3 þ y4 Þ=2

ð6Þ

Once all the effects are quantified, a number of tests are performed to investigate which effects are statistically significant and which of them are the products of chance. When there are no replicates available from the experiments, which, due to the experiment’s cost and time is not uncommon, the ordinary testing procedure is to first standardize all the effects, plot them on normal probability paper and see how they conform to the normal probability density function. Those effects that depart from normal pdf are considered to be significant contributors to the experiment’s results and their influence should be further investigated. 4.2. Two-level factorial design for the keyblock prediction model at the CLAB 2 While the methodology for factorial analysis as presented in the preceding section is rather strict, the selection of proper variables and their levels is case specific and entails a common sense principle regarding the expected results. In other words, the user can freely apply any levels to the specified variables, but at the same time the assumptions made ought to be realistic and restricted by the scale of the studied phenomenon. For example, setting up a minus level to some few centimeters for the fracture radii and its plus level beyond 1 km when in fact the observed fracture size at the site is around 1 m will shift the focus from real fracture dimensions towards extremities and end up with misleading results. Considering the computational resources available (PC: 650 MHz, 128 Mb RAM) at the time this study was pursued, and to guarantee the transparency of the analysis itself, it was resolved to construct a 23 factorial analysis, that means an experimental design for testing three variables on two levels. Relying on a number of pilot-simulations with FracMan and RockBlock prior to the design of the experiment as well as on studies reported elsewhere, we selected fracture size SZ, orientations OR and termination fraction T as the variables we expected to have the major impact on keyblock predictions.

and the interactive effect from AB: MAB ¼ ðþy1  y2  y3 þ y4 Þ=2:

ð7Þ

Table 5 Design matrix for 22 factorial experiment. A and B are model variables and Y is the result from the experiment Runs

A

B

AB

Y

1 2 3 4

 +  +

  + +

+   +

y1 y2 y3 y4

4.2.1. Devising levels for factorial experiments The plus and minus level for each variable should always be defined in relation to the experiment’s resulting parameter and its practical meaning. For the case of the CLAB site, one possible option for the resulting variable Y (see Table 5) would be to measure one or several keyblock distributional parameters and to see how they relate to the variation in fracture geometries, i.e., OR, SZ and T. However, a more interesting approach in our belief was to investigate how a two-level alteration of fracture size, orientations and termination fraction affect the discrepancy between predictions made for ‘‘the most likely model’’ (MLM)

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

for the site and the predictions made for a number of modified models, i.e., the models where variables SZ, OR and T are entered on their plus and minus levels. In more practical terms, the aim of the sensitivity analysis became to examine how inaccurate estimates of size, orientation and termination as related to their estimates in the MLM results in deviation from the most likely keyblock predictions. The properties of the MLM are summarized in Tables 1–3. While there are still some uncertainties in our estimates, for the purpose of the sensitivity analysis we assume the MLM as completely verified and perfectly reflecting the natural variability in fracture network geometry at the CLAB site. Thus, the plus and minus level in this work mimic either over or underestimation of a specific variable as compared to its estimate in the MLM. Starting with OR, while studying the stereo plot of the fracture poles in Fig. 2 there might be at least two additional interpretations as an alternative to the two orientation clusters found from ISIS analysis as depicted in Table 1. Depending on the interpreter, one could state that there are in fact three distinct orientation clusters or alternatively no distinct clusters at all and for the latter the data should be bootstrapped instead. A nonparametric bootstrap distribution [49–50] implies that generated fractures in the stochastic model are assigned the same orientations as the observed fractures but with some estimated dispersion reflecting the natural variation. In order to imitate these two possible ‘‘misinterpretations’’, i.e., (i) three orientation sets and (ii) no preferred orientation clusters, the minus level for orientation was defined as bootstrap with a very high concentration coefficient k ¼ 100 suggesting low spread from the observed orientations, and the plus level was defined as three distinct orientation sets where each set was found as the best-fit to the theoretical statistical distribution following the principle of ISIS algorithm. Setting up minus and plus levels for SZ was done assuming that no more sophisticated fracture size analysis could be done for any factorial design and either simplistic estimations of fracture radii directly from the fracture trace length data were made or infinite fracture size (the postulation originating from classical Block Theory [1]) was assumed. For the minus level, the fracture radii were assumed to be equal to the double of the mean fracture trace length observed in the tunnel floor and that they follow lognormal distribution with standard deviation equal to double the trace length’s standard deviation. Even if the last assumption was rather conservative, we considered it plausible as it was applied in another study done by Kicker [51] who modeled fracture network at Yucca Mountain tunnel. Since the modeling tool used in our work could not handle infinite fracture size, the plus level for radii was

251

Table 6 Variables and their levels for 23 factorial designs for the keyblock predictions for the CLAB 2 facility Variable

‘‘’’ level

‘‘+’’ level

OR (orientation)

Nonparametric bootstrap, k ¼ 100

SZ (fracture radii)

log N pdf, r=15.5 m, s ¼ 10:7 m 0%

Set 1: 357, 45. Fisher’s k ¼ 9:5 Set 2: 323, 13.4. Fisher’s k ¼ 23:9 Set 3: 248, 15.3. Fisher’s k ¼ 10:2 Uniform pdf, r ¼ 300 m 50%

T (termination fraction)

set up to 300 m which was deemed to be practically infinite as regards to the CLAB 2 facility dimensions and the tunnel’s location depth. Lastly, the levels for T were arranged. The minus level was set up to 0%, i.e., no terminations at fractures at all, and the plus level was set up to 50%. Table 6 summarizes the input data for the design matrix for 23 factorial analysis for the CLAB 2 tunnel. The volumetric fracture intensity was kept constant (1.55 m2/m3) for all the factorial designs. 4.2.2. Results from factorial experiments To estimate how all the combinations between two levels and among three variables affect keyblock predictions in relation to the predictions obtained from the MLM, we decided that the results from the factorial designs (column Y in Table 5) would be based on differences between the total unstable volume for the MLM and for each model from all factorial designs. This parameter was chosen for the reason that it is of practical importance for rock support requirements. Arithmetic averages for the unstable block volume from ten Monte Carlo simulations of the stochastic fracture network for all models, i.e., for the MLM and the additional eight models from 23 factorial designs were computed and are presented in Table 7. Based on the results presented in Table 7, the complete design matrix for 23 factorial experiments where the resulting variable Y is the difference between the average of total unstable block volume for MLM and each of the eight other models was derived (Table 8). The three main factors (OR; SZ and T), the three twofactor interactions (OR * SZ; OR * T and T * SZ) and one three-factor interaction (OR * T * SZ), which were derived according to Eqs, (5)–(7) are displayed in Table 9. In order to see what factors and factor-interactions significantly affect the variable Y, the calculated factors

252

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

Table 7 The average of the total unstable block volume and standard deviations for MLM and eight other models. The results for each model are based on 10 Monte Carlo runs Statistical parameters

MLM

Model 1 OR T SZ

Model 2 OR+ T SZ

Model 3 OR T+ SZ

Model 4 OR+ T+ SZ

Model 5 OR T SZ+

Model 6 OR+ T SZ+

Model 7 ORT+ SZ+

Model 8 OR+ T+ SZ+

Average (m3) Standard deviation (m3)

57 15

144 95

118 34

154 59

132 54

142 47

125 56

157 48

179 51

Table 8 The complete design matrix for 23 factorial design for the CLAB 2 with fracture orientation OR, fracture radii SZ, termination fraction T. The test result Y is the difference between total the average of unstable block volume for MLM and each of the eight models Experiment

OR

T

SZ

OR * T

OR * SZ

T * SZ

OR * T * SZ

Y (m3)

1 2 3 4 5 6 7 8

 +  +  +  +

  + +   + +

    + + + +

+   + +   +

+  +   +  +

+ +     + +

+   + +   +

87 62 97 75 85 68 100 122

Table 9 The main factors and the factor-interactions resulting from the factorial designs for the CLAB 2 site OR

T

SZ

OR * T

OR * SZ

T * SZ

OR * T * SZ

10.6

23

13.9

10.4

12.9

11.4

8.8

and interactions in Table 9 were fitted with normal pdf by implementing the Herd–Johnson method [52] with the 95% confidence interval. Fig. 4 illustrates that only one data point (factor OR) significantly departs from the normal probability plot. The diverging data point corresponded to impact of fracture orientation OR on the difference between total unstable block volume for MLM and other models. It was concluded that all other single factors and factor interactions did not have any significant impact on the measured quantity Y, and were the product of chance. As a result, the influence of OR on Y was investigated in more detail (see next section). In addition to the factorial analysis for the difference between average unstable block volume, we performed the same analysis for the median of the total unstable volume, however, no statistically significant impact of single factors and interactions on the quantity Y were observed. Also, a more sophisticated interpretation of factorial analysis was done where the difference between total unstable block volume for MLM and the other models was expressed as Kolmogorov–Smirnov and

Fig. 4. Normal probability plot (with 95% confidence interval) for main factors and interactions for 23 factor designs for the CLAB 2 regarding total unstable block volume.

Mann–Whitney statistics [18] based on the total unstable volume-data from 10 Monte Carlo realizations. Neither for the analysis based on Kolmogorov–Smirnov

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255

253

nor Mann–Whitney statistics there was any significant impact of single factors or interactions on Y. 4.2.3. The effect of fracture orientations on the total unstable block volume The statistically significant impact on the difference between total unstable block volume for MLM and for the other models encouraged further analysis. Since the OR impact on the quantity Y proved to be negative (10.6 according to Table 9), the shift from the OR(i.e., bootstrap with no clusters) to OR+ (three sets, each represented with Fisher’s pdf) implies lower Y and consequently lower discrepancy between MLM and the models with three orientation clusters. Thus, the total unstable block volume predicted from models where fracture orientations were approximated with three clusters was more accurate as regards to MLM than the predictions based on models with bootstrapped orientations.

5. Conclusions Stochastic fracture network built from 2D fracture mapping in the CLAB 2 facility floor was derived. Orientation of fractures was approximated with two major clusters each following Fisher’s pdf, fracture size, i.e., the representative radius was best fitted with lognormal distribution and fracture locations followed Poisson process. The proposed model was used for simulation of unstable blocks along the full length of the facility. The simulation included 10 Monte Carlo realizations and resulted in the average of 38 unstable blocks where their volume was best-fitted with the Pareto distribution. As the distribution of the block volume and the average block amount are of less practical value the estimation of the total reinforcement required should rely on other predicted parameters such as the total unstable volume Vt ¼ 57 m3 or the largest unstable block Vmax ¼ 12 m3. Since the predicted largest block volume showed by far more spread around the average than the total unstable volume, the last parameter was considered to be a more correct measure of instability of the cavern. The presented factorial analysis illustrated that incorrect estimation of fracture orientation pattern can have most severe consequences in terms of correct prediction of the stability conditions. The prediction of total unstable block volume based on models with three orientation clusters approximated with Fisher’s pdf, i.e., plus level in factorial design was clearly closer to the prediction based on Most Likely Model than prediction derived from models where orientations were represented with non-parametric bootstrap distribution (minus level) as Fig. 5 indicates. Fig. 5 shows that also OR * T * SZ; OR * T and SZ depart from the normal

Fig. 5. The impact of fracture orientations OR on measured quantity Y (deviation of the predicted total unstable volume from the most likely model for the CLAB 2 facility).

probability plot, however, these deviation is considerably smaller than OR and rather product of chance. The impact of three factor-interactions on the block prediction is rather a complex case, nonetheless, the quick look at Table 8 leads to the conclusion that Model 2 with OR+, T and SZ- gives the closest match to the MLM. The obtained results are case specific, however, we conclude that in general, the accuracy in the prediction of total unstable volume is most sensitive to correct estimation of fracture orientations. Factorial analysis has a potential for more stringent sensitivity analysis and may be used for an optimal variable estimation for stochastic fracture models.

Acknowledgements This study was supported with Grant 821(621) from the Swedish Rock Engineering Research (SveBeFo). Golder Assoc. Seattle is appreciated for permission to use RockBlock code. The Swedish Nuclear Fuel and Waste Management Co. (SKB) provided us with dataset and all written documentation on the CLAB site. We thank Editor JA. Hudson and two anonymous reviewers for useful suggestions.

References [1] Goodman RE, Shi G. Block theory and its application to rock engineering. Englewood Cliffs, NJ: Prentice-Hall, 1985. [2] Shapiro A, Delport JL. Statistical analysis of jointed rock data. Int J Rock Mech Min Sci Geomech Abstr 1991;28(5):375–82. [3] Mauldon M. Relative probabilities of joint intersections. In: Tillerson JR, Wawersik WR, editors. Proceedings of the 33rd US Symposium on Rock Mechanics, vol. 33. Santa Fe, June 3–5, 1992. p. 767–74. [4] Stone CA, Kuszmaul JS, Boontun A, Young D. Comparision of an analytical and a numerical approach to probabilistic keyblock analysis. In: Aubertin M, Hassani F, Mitri HS. Proceedings of the

254

[5]

[6]

[7] [8]

[9] [10]

[11] [12]

[13]

[14]

[15]

[16]

[17]

[18] [19]

[20]

[21]

[22] [23] [24]

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255 Second North American Rock Mechanics Symposium, vol. 2, Montreal, June 19–20, 1996. p. 1769–75. Kuszmaul JS. Estimating keyblock sizes in underground excavations: accounting for joint set spacing. Int J Rock Mech Min Sci 1999;36:217–32. Priest SD, Samaniego JA. The statistical analysis of rigid block stability in jointed rock masses. Proceedings of the Fifth Australia–New Zealand Conference on Geomechanics, The Institutions of Engineers Australia, Sydney, 1998. p. 398–403. Hoek E, Carvalho J, Li B. Unwedge users manual. Toronto: Rocscience Inc., 1991. Li B. The stability of wedges formed by three intersecting discontinuities in the rock surrounding underground excavations. Ph.D. dissertation, University of Toronto, 1991. Stone CA. A matrix approach to probabilistic keyblock analysis. Ph.D. dissertation, Michigan Technological University, 1994. Jakubowski J. Prediction of the load of tunnel support in the rock mass of blocky structure by statistical methods. Ph.D. disserta! tion, University of Mining and Metallurgy: Krakow, Poland, 1995. PanTechnica Corporation, PT workshop, KBTunnel module, version 2.0 users manual. MN: Eden Prairie, 2000. Cundall PA. Formulation of a three-dimensional distinct element model-Part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. Int J Rock Mech Min Sci Geomech Abstr 1988;25(3):107–16. Shi G-H. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block system. Ph.D. dissertation, University of California, Berkeley, 1988. Chan L-Y. Application of block theory and simulation techniques to optimum design of rock excavations. Ph.D. dissertation, University of California, Berkeley, 1986. Hoerger SF, Young DS. Probabilistic prediction of keyblock occurrences. In: Hustrulid, Jhonson, editors. Rock mechanics contributions and challenges. Rotterdam: Balkema, 1990. p. 229–36. Jakubowski J, Tajdus A. The 3D Monte-Carlo simulation of rigid blocks around a tunnel. In: Rossmanith HP, editor. Proceedings of the Second International Conference on the Mechanics of Jointed and Faulted rock, Vol. 2. Rotterdam: Balkema, 1995. p. 551–6. Dershowitz W, Lee G, Geier JE, Foxford T, La Pointe P, Thomas A. FracMan. Interactive discrete feature data analysis geometric modeling and exploration simulation. User documentation version 2.6. Seattle: Golder Assoc. Inc., 1998. Davis JC. Statistics and data analysis in geology. New York: Wiley, 1986. Kulatilake PHSW, Fiedler R, Panda BB. Box fractal dimension as a measure of statistical homogeneity of jointed rock masses. Eng Geol 1997;48:217–29. . o. Hard Rock Laboratory. TRUE Dershowitz W, Ushida M. Asp Block Scale project. Spatial data analysis. Swedish Nuclear Fuel and Waste Management Co., ITD-99-30. Stockholm, 1999. Starzec P, Andersson J. Stochastic fracture modeling for stability assessment of underground facilities in crystalline rocks. In: Elsworth D, Tinucci JP, Heasley KA, editors. Proceedings of the 38th U.S. Rock Mechanics Symposium, DC Rocks 2001, vol. 2. Washington DC, 7–10 July 2001, p. 1483–90. Berglund J. personal communication. Shanley RJ, Mathab MA. Delineation and analysis of clusters in orienattion data. J Math Geol 1976;8(3):9–23. Mathab MA, Yegulalp TM. A rejection criterion for definition of clusters in orientation data. In: Goodman RE, Heuze FE. editors. Issues in Rock Mechanics, Proceedings of the 22nd Symposium on Rock Mechanics, American Institute of Mining Metallurgy and Petroleum, Berkeley, p. 116–23.

[25] Priest SD. Discontinuity analysis for rock engineering. London: Chapman & Hall, 1993. [26] La Pointe P, Wallman PW, Follin S. Estimation of effective block conductivities based on discrete fracture network analysis using . o. site. Swedish Nuclear Fuel and Waste data from the Asp Management Co., TR 95-15. Stockholm, 1995. [27] Swan ARH, Sandilands M. Introduction to geological data analysis. London: Blackwell Science, 1995. [28] RockWare Inc. RockWorks98 manual. Golden CO: RockWare Inc., 1998. [29] Fisher R. Dispersion on a sphere. Proc Roy Soc London 1953;A217:295–305. [30] Warburton PM. A stereological interpretation of joint trace data. Int J Rock Mech Min Sci Geomech Abstr 1980;17:181–90. [31] La Pointe P, Wallman P, Dershowitz W. Stochastic estimation of fracture size through simulated sampling. Int J Rock Mech Min Sci Geomech Abstr 1993;30:1611–7. [32] Press W, Teukolski S, Vetterling W, Flannery B. Numerical recipes in C, the art of scientific computing. Cambridge UK: Cambridge University Press, 1992. [33] Priest SD, Hudson JA. Discontinuity spacings in rock. Int J Rock Mech Min Sci Geomech Abstr 1976;13:135–48. [34] Baecher GB, Lanney NA, Eistein HH. Statistical description of rock properties and sampling. In: Wang FD, Clark GB, editors. Energy resources and excavation technology. Golden CO: Colorado School of Mines, 1977. p. 5C11. [35] Barton C, La Pointe P. Fractals in the Earth Sciences. New York: Plenum Press, 1995. [36] Chil!es JP. Fractal and geostatistical methods for modelling of a fracture network. Math Geol 1988;20(6):631–54. . o. site: [37] Geier JE, Thomas AL. Discrete-feature modeling of the Asp discrete fracture network models for the repository scale, vol. 5. Stockholm: Swedish Nuclear Power Inspectorate, 1996. [38] Tansi C, Sorriso-Valvo M, Greco R. Relationship between joint separation and faulting: an initial numerical appraisal. Eng Geol 2000;52:225–30. [39] Dershowitz W, Herda H. Interpretation of fracture spacing and intensity. In: Tillerson JR, Wawersik WR, editors. Rock mechanics; Proceedings of the 33rd US Symposium, vol. 33. 1992. p. 757–66. [40] Dershowitz W, Carvahlo J, Foxford T. Fracman/rockblock. discrete fracture rock block stability analysis. User documentation. Version 1.0. Seattle: Golder Assoc. Inc., 1995. [41] Stille H, Fredriksson A. CLAB-Etapp 2: Bergmekanisk utredning . allanden och inverkan p(a befintligt bergrum. av stabilitetsforh( Swedish Nuclear Fuel and Waste Management Co., PR 96-06, Stockholm, 1996. [42] Barton N, Bandis SC. Review of predictive capabilities of JRCJCS model in engineering practice. In: Barton N, Stephansson O, editors. Proceedings of the International Symposium on Rock Joints, Loen, Norway. Rotterdam: Balkema, 1990. p. 603–10. [43] Chan L-Y, Goodman RE. Predicting the number and dimensions of key blocks of an excavation using Block Theory and joint statistics. In: Farmer IW, Daemen JJK, Desai CS, Glass CE, Neuman SP, editors. Proceeding of the 28th US Symposium on Rock Mechanics, University of Arizona, Tuscon, 1987. p. 81-87. [44] Wainwright E, Yastrebov A, James D. Crystal Ball 2000 professional edition. Denver CO: Decisioneering, 2000. [45] Box GEP, Hunter WG, Hunter JS. Statistics for experimenters. New York: Wiley, 1978. [46] Daniel C. Applications of statistics to industrial experimentation. New York: Wiley, 1976. [47] Hunter WG, Chacko E. Increasing industrial productivity in developing countries. Int Devt Rev 1971;13(3):11. [48] Hunter JS, Naylor TH. Experimental designs for computer simulation experiments. Manage Sci 1971;16:422.

P. Starzec, J. Andersson / International Journal of Rock Mechanics & Mining Sciences 39 (2002) 243–255 [49] Efron B. The Jacknife the bootstrap and other resampling plans. SIAM Monograph, vol. 38. Penn, PI: Society of Industrial and Applied Mathematics, 1982. [50] Hjort U. Computer intensive statistical methods: Validation model selection and bootstrap. London: Chapman & Hall, 1994.

255

[51] Kicker D. Drift degradation analysis. ANL-EBS-MD-000027, Office of Civilian Radioactive Waste Management, Las Vegas, NV, 2000. [52] Ryan BF, Joiner BL. Minitabt handbook 4th ed. Canada: Duxbury Thomson Learning, 2001.