Geomechanical data are never sufficient in quantity or adequately precise and accurate for design purposes in mining and civil engineering. This is because rock masses are naturally complex and variable at all scales. Qu (2017) stated that, in theory, the geomechanical characteristics of rock masses are not completely random variables. This is because rocks were formed and continuously modified by a variety of complex processes, causing physical heterogeneity that results in variations in measured physical properties, even within one rock type. Moreover, the presence of natural fractures creates spatial and directional variations in rock mass properties, i.e., such fractures cause a rock mass to become inhomogeneous and anisotropic (Jing 2013). Therefore, in geomechanical investigations and laboratory tests, it is necessary to statistically analyze the parametric characteristics of rock as random variables (Yegulalp and Mahtab 1983; Uzielli 2008; Mayer et al. 2014).

The aim of statistical analysis is to find the most representative value of the parameter needed. This value can be selected based on main statistical indices and on normal distribution parameters. In geomechanics, such analyses are performed to assess the heterogeneity and anisotropy of rock masses and the corresponding variation in the laboratory test results of physical properties of rocks. Increasing the number of data usually increases this variability, and the samples become more representative for the entire population (Joughin 2017). Many researchers have used complex data analysis methods and complex statistical models (Briševac et al. 2016; Pandit et al. 2019; Babets et al. 2019), including fuzzy models (Gokceoglu and Zorlu 2004), Monte Carlo simulations (Fattachi et al. 2019), Bayesian models (Feng and Jimenez 2014; Wang and Aladejare 2016), and hierarchical cluster analysis (Mayer et al. 2014). All of these methods lead to the derivation of a specific parameter, such as uniaxial compressive strength or Young’s modulus, and allow us to assess the uncertainty of the data set.

There are no standardized requirements for the scope of site investigation prior to mining or tunneling in hard rocks. Therefore, the sampling of rock masses is usually very inadequate due to the costly, time consuming, and technically complicated drilling required and the often limited access to the work site. In addition, rock masses are often heavily fractured, which prevents the collection of samples suitable for laboratory testing (e.g., tensile strength TS or uniaxial compressive strength UCS) (Fig. 1). Under such conditions, even cores that are 8–10 m long may produce few samples. Accordingly, the information on rock beds is sparse and insufficient. In contrast, a reliable predictive model requires a large number of high-quality data. There are known cases when even extensive site investigations could not provide rock samples usable for laboratory testing. Gokceoglu and Zorlu (2004) mentioned that even though approximately 200 blocks were collected from the field, only 82 of the obtained sample sets were able to be used for rock mechanics tests.

Intensely crushed core from the floor of the work site

In some instances, there are no core pieces that are sufficiently long to prepare samples for laboratory testing (Fig. 2). The fewer boreholes and samples that exist, the lower the quality of the geomechanical research and statistical analysis, which are both indispensable in rock mass numerical modeling.

Entirely crushed core from the floor of the work site

Another issue in geomechanical research is that even for the same type of rock beds, at very close distances, the strength parameters and Young modulus commonly vary by a factor of 2–3 (Majcherczyk et al. 1999) and sometimes vary more considerably by factors up to 5–8 (Dychkovskyi et al. 2020).

The objective of this paper is to determine the variability in the rock properties at a selected point in the roof of a mining roadway and then illustrate how the statistical approach used for the available geomechanical data can influence the results of numerical modeling.

Four statistical approaches were applied in this study: average values (the most common approach in geomechanical data analyses and project applications), average minus standard deviation, median, and average value minus statistical error. The generated maps of stress, displacement and yield zone showed that the most substantial difference in the numerical results was related not to the variation in the available data but to the statistical approach used for the data set.

2.1 Derived physical parameters of rocks and their statistical parameters

The laboratory tests were carried out on sedimentary rocks that were sampled from the roof of Roadway B-7 in the Pniowek coal mine. The samples were cut from three different cores, T-87/07, T-88/07 and T-89/07, which were drilled at a chainage of 248 km in the mining roadway. The distance between the core holes was approximately 1.4 m, and the axes of both outer holes were tilted from the vertical toward the roadway walls at an approximate angle of 10°–15° (Fig. 3). This layout ensures that the coring through the roof rocks covers the whole 6.1 m width of Roadway B-7. Worldwide research shows (Tien and Kuo 2001; Tavallali and Vervoort 2010; Khanlari et al. 2014; Yao et al. 2019) that coring at a small inclination to the normal direction of rock beds—up to 10°–13°—should not have an effect on the mechanical parameters of samples and can be neglected.

Lithology of the boreholes drilled in Roadway B-7

Three different rock types were identified through geological logging: claystone, mudstone and fine-grained sandstone (Fig. 3). It is worth noting that even though all three core holes were drilled at the same chainage in Roadway B-7, the thicknesses of the encountered rock beds were substantially different from core hole to core hole.

The general objective of the tests was to verify the assumptions for the design of roadway support. Among others, bulk unit weight (γ), uniaxial compressive strength (UCS), Young’s modulus (E), and tensile strength according to the Brazilian test (TS) have been tested.

The preparation of samples and the measurement procedures fulfilled the International Society of Rock Mechanics (ISRM) standards (Ulusay and Hudson 2006). The sample diameter was 50 ± 1 mm (the same as that of the core).

Young’s modulus was derived as an E_tan in the range of 20%–80% of the ultimate stress (Eq. (1), Małkowski et al. 2018). This range corresponds to the linear part of the stress–strain curves of the tested rocks.

where F is the maximum force exerted on the sample, N; A is the original cross-sectional area upon which the force is applied, m²; L is the length of the sample, m; and ΔL is the sample length reduction, m.

2.2 Laboratory test results

The following results of unit weight, compressive strength, tensile strength and Young’s modulus for the rocks were obtained in accordance with the relevant ISRM standards. Statistical analysis was performed on the data set from each rock type.

The number of samples cut from cores T-87/07, T-88/07 and T-89/07 was 24, 30 and 26, respectively, which contained different amounts of specific rock types (Fig. 4). Due to this shortage of samples, there were only 5 or 6 tensile strength tests for each rock type identified in the core holes and only 4–5 Young modulus tests. A small number of tests affects the validity of the statistical analysis. Although the core holes were 8 m long, the obtained core was fractured, and it was not possible to cut more than 20–30 samples from one core for laboratory testing, which could be a desirable number for statistical analysis. This is a typical situation in research on sedimentary rocks in mining environments.

Different rock types used in the laboratory compressive strength tests

The statistical analysis of the laboratory test results confirmed the rock heterogeneity. The ranges of compressive strength for claystone, mudstone, and sandstone were 59–147, 48–211 and 14–187 MPa, respectively (Fig. 5); hence, there were substantial variations (3–4-fold or greater) between the compressive strength ranges of the different rocks. However, the average compressive strength values of claystone, mudstone and sandstone were quite similar at approximately 90 MPa; the different rocks also had a similar standard deviations of approximately 31–38 MPa. Due to the high range of values, the rocks had high coefficients of variation of 37%–39%. The average tensile strength for the specific rock types varied from 5.5 MPa up to 8.6 MPa (Fig. 6). Testing 4–8 samples for tensile strength immediately influenced the variations in the standard deviation and the coefficient of variation. The same difference was observed when studying the Young’s moduli of the different rocks (Fig. 7). Here, the small number of samples was due to the difficulty in finding core sections sufficiently long to satisfy the standard sample specification for this test, i.e., the required diameter-to-length ratio was 1:2. The most stable statistical results were obtained for unit weight even though the number of samples was different (Fig. 8). This stability was verified by the corresponding coefficient of variation, which was only 3%–5%.

Comparison of chosen statistical parameters for compressive strength data

Comparison of chosen statistical parameters for tensile strength data

Comparison of chosen statistical parameters for Young's modulus data

Comparison of chosen statistical parameters for unit weight data

In summary, the results show that the geomechanical parameters were highly variable, especially in the case of interlayers inside the rocks. In this case, the difference between the values was as great as 200%, but the difference could potentially be greater than tenfold for compressive strength. Such variation in results is commonly observed for heterogeneous sedimentary rocks (Palchik 2011). For this reason, advanced statistical techniques are employed for geomechanical data analysis (Yegulalp and Mahtab 1983; Gokceoglu and Zorlu 2004; Wang and Aladejare 2016).

2.3 Case studies for different statistical approaches

The obtained results of the geomechanical parameters were typical of design projects conducted in sedimentary rocks. If we consider the support design and roadway (or tunnel) stability analysis, the main concern was the statistical preparation of the input data. The conservative approach, wherein unfavorable values are assumed for the rock parameters, leads to “overdesign” and cost increases, whereas assuming excessively favorable values can be risky for the designed structure. Day et al. (2017) claimed that in some instances, changes in material parameters led to an increase in the global factor of safety, whereas the reliability index for these parameters decreased. In engineering practice, advanced statistical analyses are very seldom conducted; rather, simple statistical parameters are used, such as standard deviation or statistical error.

When choosing the most suitable scenarios for numerical analysis, we assumed that the most common approach for geomechanical data was the average value (Case 1). Then, we also assumed that the most common confidence level used in statistical investigations was 95%. Therefore, for this value, we calculated the statistical error. Furthermore, we analyzed three other cases: average minus standard deviation (Case 2), average minus statistical error (Case 3), and median (Case 4). We chose these variants because they are quite easy to investigate in engineering practice—as recommended by Hammah and Curran (2009)—and because they can consider the uncertainty in geomechanical data.

After performing the analysis above, we determined the percentage difference between Case 1 (average value) against the other three cases (Fig. 9).

Difference between the geomechanical parameters of Case 1 (average) and the other three cases: a compressive strength, b tensile strength, c unit weight, and d Young’s modulus

Figure 9 shows significant differences between Case 1 (average) and the other three cases. The greatest difference appears between Case 2 (average minus standard deviation) and Case 1. The difference can reach nearly 40% for the compressive or tensile strength. The lowest difference appears for unit weight, which is approximately 1%. Note that Case 4 (median) is usually close to Case 1, but the reason for this agreement could be the small number of samples.

The geomechanical parameters obtained from the four selected statistical approaches (Cases 1–4), which were further used in the numerical modeling, are shown in Table 1.

The variations in the measured geomechanical parameters of the rocks reveals how substantially the statistical approach can affect the state of stress and strain when the data set is input to the numerical model. Numerical methods are widely used for solving a variety of technical problems (Studeny and Scior 2009; Jing 2013; Sun et al. 2015; Smith et al. 2019; Li et al. 2020). However, Feng and Hudson (2010) emphasized that selecting the correct methods for obtaining the data required for the selected modeling method was one of the steps of appropriate design procedures in geomechanics.

The problem defined above has been checked with the use of Phase 2 software. In the model, the roadway cross-section (Fig. 10) was based on the roof lithostratigraphic profile logged at 248 km in Roadway B-7 and on the floor rock profile obtained from other site investigations. The other assumptions used in the numerical model were as follows: the field of stress used in the modeling was at a depth of 790 m; the size of roadway B-7 was as: w = 5.5 m and h = 3.8 m; the vertical boundaries were constrained from movement in the horizontal direction, and the bottom and top boundaries of the model were fixed; the rock mass was an elastic–plastic material, which is a fundamental representation of rock masses (Jing 2013); the generalized Hoek–Brown failure criterion was applied (Jing 2013); the dilation parameter was one-third of the m_b parameter for weaker rocks and two-thirds of the m_b parameter for stronger rocks (Phase 2 Tutorial); the convergence criterion used was absolute energy; the residual parameters were reduced by 20%.

Numerical model

There are many options and details in numerical modeling (Jing 2013, Malkowski et al. 2018) that need to be selected appropriately, including the geomechanical properties of the rocks, a proper physical model describing the behavior of the rock mass, and the adopted failure criterion.

However, the objective of this particular modeling was only to show the difference in the results generated from the same geomechanical data set—used for the input to the model—after processing with different statistical approaches. Four kinds of maps were generated and studied: major principal stress σ₁ and damage zone, vertical stress σ_z, total displacement, and yielded elements.

Analysis of the major principal stress σ₁ and damage zones in Case 1 (average) revealed a gradual increase in claystone and coal around the roof, ribs and floor of the excavation (Fig. 11). The shear and tension zones occurred particularly in coal and claystone above the roof of the excavation and in the floor. The claystone in the roof was sheared at a distance up to several meters from the roadway contour. The rocks around the roadway were entirely distributed. Looking further at Case 3 (Fig. 13), an increase of approximately 1 m in the affected regions around the excavation was observed, particularly around the ribs, roof and floor. The shear zones in mudstone were wider, and in general, the distressed zone was larger. In Case 4 (Fig. 14), no additional layers—upward or downward—were affected, but the lateral spread of the affected zone increased along the layers (Fig. 12).

Major principal stress σ₁ with yielded elements of Case 1

Major principal stress σ₁ with yielded elements of Case 2

Major principal stress σ₁ with yielded elements of Case 3

Major principal stress σ₁ with yielded elements of Case 4

The strongest effect was observed in Case 2 (Fig. 12), wherein a substantial increase in the distressed zone occurred in the fine-grained sandstone layer, and the whole coal layer was sheared. In this case, the stress concentration on the contact of sandstone and claystone in the roof was the highest among the four cases (Fig. 14).

The total displacement of the roadway contour was up to 0.28 m at the ribs in Case 1 (Fig. 15). In Case 2 (Fig. 16), although the total displacement was reduced by approximately 0.1 m, it was larger at the roof and floor part of the excavation contour. In Case 3 (Fig. 17), a further increase in the total displacement up to 0.32 m was observed. The highest displacement appeared in Case 4, reaching up to 0.35 m in the ribs (Fig. 18). Importantly, the lowest displacement in the floor was observed in Case 1.

Total displacement of Case 1

Total displacement of Case 2

Total displacement of Case 3

Total displacement of Case 4

Overall, the application of the four different cases to the numerical model verified how important the statistical approach was to the geomechanical data set. The results derived from these cases had substantial differences (Table 2) in the value of stress and the shape and extent of the damage zone.

The largest differences were observed in the total displacement of the excavation contour and in the range of the yielded zone. From studying the stress and displacement distribution around the roadway at specific points in the model, a large difference in the values was found between the analyzed cases. In some instances, this difference was negligible, whereas in other instances, the difference was greater than 100%, which would dramatically change the support design. For example, the displacement in the roof can vary from 5 cm up to 13 cm, and in the rib section, it can vary from 0.24 m up to 0.34 m, so the predictable horizontal convergence can vary from 48 to 68 cm depending on the case used for geomechanical data processing. This immediately suggests different designs for rock mass supports.

It is important to note that no support was applied in the analyzed model, which could considerably reduce the yielded zone range. However, the rationale of this approach was that the initial modeling without any support can provide basic information regarding what type of support was needed. This tells, for example, what length of bolts should be applied, what amount of shotcrete can be used, and whether cable bolts or steel sets are needed. The final decision depends on the needs and engineering solutions used at the mine and can be further tested in the numerical model.

The authors of this study realize that laboratory geomechanical data without any reductions were applied in the numerical analysis, so these data do not reflect the heterogeneity of the rock mass (Bello 1988). However, the aim of this study was to show the difference in the results from numerical models when different statistical approaches were used to process the data set. The issue of how to reduce the geomechanical properties of rocks derived in a laboratory is a problem that must be addressed in another study (Kulatilake et al. 2004). Moreover, further recalculations of the rock mass parameters would obscure the results of the conducted numerical study on the selected cases.

1. (1)

   Usually, the number of samples that can be prepared from a single core is limited to 20–30 from an 8-m-long core. This is a typical situation in engineering practice. Accordingly, more than one borehole must be drilled in the same location, and more cores must be investigated.

2. (2)

   The variation in the laboratory data showed that the same rock type sampled from different locations can be significantly different geomechanical parameters. The high level of variation also highlighted the sample heterogeneity. For sedimentary carboniferous rocks, the variation in the obtained data can be as high as tenfold, even for a single rock type. This variability is especially problematic for the rock strength.

3. (3)

   The approach used for statistical processing of the geomechanical input data had a substantial effect on the numerical modeling results of Roadway B-7 in Pniowek coal mine. The analysis showed that the most visible difference was in the extent of the destressed zone, which was also reflected in the yielded part of the rock mass and in the displacements of the roadway contour. Case 2 (average minus standard deviation) was the most conservative and least risky option for modeling. Using Case 1 (average) was the least conservative approach, and the cost of support will be low. This is probably the reason why it is typically used in design procedures. The differences between the above cases can reach up to fourfold in displacement and in the yielded element zone. Although, if we consider the safety of the working environment as a priority, we may opt to choose Case 2 as the basis for the support design; however, the cost of support will be rather high when selecting this approach.

4. (4)

   According to the study, taking the average values minus the statistical error (Case 3) for the numerical analysis seems to be the best solution because the median value (Case 4) always depends on the number of tests and the symmetry of the distribution. Using Case 3 for numerical modeling, the confidence level can be adjusted to the meaning of the designed object.

5. (5)

   The results obtained in this study should be compared with the real conditions at the mine site. This is the only way to calibrate the numerical model, which will help indicate the best statistical approach for the model input data.