1 Introduction

Joint roughness affects the mechanical and hydraulic behavior of fractured rock masses (Bae et al. 2011; Chen et al. 2021; Lê et al. 2018; Yong et al. 2018d; Zhang et al. 2017; Zhao et al. 2018). During the last 4 decades, considerable efforts have been devoted to the quantitative description of the surface roughness of rock joints (Barton 1978; Fardin 2007; Grasselli et al. 2002; Kulatilake and Um 1999; Liu et al. 2017; Maerz et al. 1990; Morelli 2014; Tse and Cruden 1979; Zheng and Qi 2016; Ban et al. 2021). However, given the variability and limited accessibility of natural features, achieving an exact characterization of in situ fractures is generally difficult (Vogler et al. 2017). Natural discontinuities often show spatial variations of roughness (Alameda-Hernández et al. 2014; Wu et al. 2020). The uncertainty in the joint roughness measurements greatly affects the estimation of the rock joint shear strength and further influences the evaluation of the rock mass stability.

Based on rock joint roughness investigations in the field, Du (1998) summarized the following three characteristics of joint roughness: (1) heterogeneity, (2) anisotropy, and (3) nonuniformity. Heterogeneity indicates that the joint roughness measurement varies with position, even when measured on the same joint surface and along the same measurement direction, because the roughness varies considerably along different cross-sections of a joint surface. Anisotropy indicates that the joint roughness behaves differently in different directions, even when measured on the same joint surface, which is attributed to the formation process of rock joints. Nonuniformity indicates that the joint roughness is dependent on the type of rock, mineral grain size, and the weathering conditions. Joint roughness estimations are generally performed on a set of joint profiles measured parallel to the direction of the shear displacement for a given rock joint, and the average value is commonly used to represent the roughness of the entire joint surface (Diaz et al. 2017; Du 1999; Liu et al. 2017; Lopez et al. 2003; Tang et al. 2015; Xia et al. 2014; Zhang et al. 2016). The International Society for Rock Mechanics (ISRM) advised obtaining surface roughness from rock joint samples located at different positions and of varied length scales (Ulusay and Hudson 2006). These studies have already considered the influence of heterogeneity on joint roughness estimations and are helpful for obtaining representative evaluation results. However, the areal extent of the exposed joint surface in situ is generally limited, so researchers often have to evaluate the joint roughness based on insufficient samples. Sampling bias is an important issue in statistics and likely leads to incorrect predictions (Spiegel and Stephens 2017). Thus, researchers need to know how many individuals they must survey to draw conclusions about the population. In this study, the population refers to all samples extracted from different locations of a given rock joint, and an individual refers to one of the potential samples. In previous studies, investigators have not paid enough attention to the influence of the sample number, and the problem of estimating the joint roughness based on an inadequate number of samples remains unsolved.

Many investigations have been conducted on the scale effect of shear strength through laboratory and in situ tests (Bandis et al. 1981; Fardin et al. 2001; Hencher and Richards 2014; Tatone and Grasselli 2012; Ye et al. 2016). It was found that the selection of testing samples affects the observation of scale effects, and sampling bias may lead to incorrect conclusions about the scale effect (Barton 1990). For example, representative joint samples of different sizes are needed to accurately describe the scale-dependent behavior of shear strength. In addition, the average value of the samples subdividing a large joint is used to represent the shear strength of each size, and the relationship between the average value and sample size is used to reveal the scale effect on the shear strength of rock joints (Bahaaddini et al. 2014; Bandis et al. 1981; Morelli 2014). However, the number of joint samples used in previous studies sometimes appears to be insufficient for statistically investigating the scale effect on shear strength (Yong et al. 2018b), which may lead to the inability to estimate the representative characteristics of the scale effect. In practice, the evaluation of structurally controlled slope instability requires careful investigations of the shear behavior of the potential sliding plane (Huang et al. 2019; Tang et al. 2017; Yang et al. 2019a, 2019b). Thus, the representative roughness of the potential sliding plane should be estimated. Systematically characterizing the heterogeneity of joint roughness is expected to provide insights into the exploration of representative samples. Unfortunately, methods for handling the heterogeneity of the surface roughness of rock joints are still lacking.

To address these concerns, the global search method was introduced to characterize the heterogeneity based on a statistical analysis of the roughness of all samples extracted from different locations of a given rock joint. The roughness heterogeneities of both two-dimensional (2D) profiles and three-dimensional (3D) surface topographies were quantitatively evaluated, and then this proposed method was applied for searching representative samples.

2 Heterogeneity of Joint Roughness and Its Investigation in the Field

The accuracy of heterogeneity investigations depends on several factors, such as the accessibility of the plane of interest, the areal extent of the exposed plane, and human error. However, generally, very limited portions of potential sliding planes are exposed in the field. Thus, it is practical to study the heterogeneity of the rock joint roughness based on the exposed part of the potential sliding plane. In rock engineering practice, three situations of exposed potential sliding planes are commonly encountered in field surveys (labeled situations A, B, and C in Fig. 1).

Fig. 1
figure 1

Three typical situations of the exposed sliding plane commonly encountered in the field. a Situation A, b situation B, c situation C

Situation A reflects the case in which a very limited area of a rock joint surface is exposed, and the direction of the long edge of the exposed area is approximately the same as the sliding direction (Fig. 1a and Fig. 2). In this situation, measuring large joint profiles is feasible, and the scale effect on the joint roughness can be investigated. However, it is difficult to obtain many profiles parallel to the sliding direction. To study the heterogeneity, Du (1994, 1999) suggested that the joint roughness should be determined based on the segments (S1, S2, S3, ···) obtained by evenly dividing the measured large joint profile (Fig. 2). The roughness of the entire joint profile is represented on the basis of the average joint roughness coefficient (JRC) of each segment.

Fig. 2
figure 2

Arrangement of the measured joint profiles under different situations

Situation B also describes the case where a small area of the rock joint surface is exposed, but the direction of the short edge of the exposed area is the same as the sliding direction (Figs. 1b and 2). In this situation, sufficient profiles (P1, P2, P3, ···) can easily be obtained along the sliding direction (Fig. 2). However, it is not feasible to measure the large joint profiles along the sliding direction in this case. The obtained profiles are too short to investigate the scale effect on the joint roughness. In this situation, the heterogeneity reflects the roughness diversity of small joint profiles (P1, P2, P3, ···) measured at different positions parallel to the sliding direction, and the average roughness of these profiles was used to represent the overall roughness of the studied rock joint.

Situation C describes the case in which a relatively large area of the rock joint surface is exposed (Figs. 1c and 2). In contrast with situations A and B, situation C offers a better condition for measurement, and therefore, obtaining sufficient joint profiles is easier to achieve. In this situation, the heterogeneity of the joint roughness is relatively comprehensively revealed. As shown in Fig. 2, a number of closely spaced long profiles can be easily measured, and the variation of joint roughness along the parallel profiles can be obtained to characterize the roughness heterogeneity.

For samples with small sizes, acquiring many cross-sections along the potential shear direction is difficult to achieve in situation A. This study focuses on the sampling problem in situation A, and a new method for acquiring sufficient samples is offered to exploit detailed roughness information about the heterogeneity of the joint roughness.

3 Global Search Method for Investigating the Heterogeneity of Joint Roughness

The conventional sampling methods include the following three categories: (1) the equal‑partition sampling method, (2) the simple random sampling method, and (3) the processive magnifying sampling method. Their definitions and limitations are summarized as follows:

The equal‑partition sampling method refers to the sampling process by which smaller sized samples are obtained using equal partitioning of a larger sample. As shown in Fig. 3, we can obtain 100 mm, 200 mm, and 500 mm samples using equal partitioning of the entire joint profile. The obtained samples of each sample size cover the entire length of the original joint. However, for the samples with lengths of 300 mm, 400 mm, 600 mm, 700 mm, 800 mm and 900 mm, the end parts of the entire joint must be abandoned, because the whole length of the joint profile cannot be equally divided by these sample sizes. For example, when a 600-mm-long sample is taken from the entire joint profile, the length of the abandoned part reaches 40% of the overall length of the original joint. Thus, the roughness of the obtained 600 mm sample can only represent the roughness characteristics of the entire joint in the first 600 mm of length. The roughness characteristics of the abandoned part are neglected. Thus, sampling bias occurs as a result of the samples obtained unable to represent all the samples at different locations on the joint.

Fig. 3
figure 3

taken from the original surface by the partition sampling method

Samples of different sizes

The simple random sampling method refers to the sampling process in which samples of different sizes are arbitrarily taken from the original surface. Compared with the equal‑partition sampling method, the locations of different-sized samples obtained using the simple random sampling method are random and irregular, and only one sample is selected for each sample size. Figure 4 shows three examples of the sampling schemes obtained based on the simple random sampling method. The locations of the samples obtained in schemes A, B, and C are significantly different, and the morphology of the samples of the same size are also different. The representativeness of samples primarily depends on the personal judgment and choice of researchers, so the sampling results are generally difficult to be reproduced.

Fig. 4
figure 4

taken from the original surface by the simple random sampling method

Samples of different sizes

The processive magnifying sampling method refers to the sampling process in which smaller sized samples are obtained by an edge cutout from a side section or by cutting a sample from the middle section. Three sampling results are obtained based on the processive magnifying sampling method, as shown in Fig. 5, where scheme A acquires samples starting from the left edge, scheme B acquires samples starting from the right edge, and scheme C acquires samples starting from the middle. All smaller sized samples are part of larger samples. Compared to the simple random sampling method, the processive magnifying sampling method has explicit sampling rules by which sampling results can be reproduced by different researchers. However, the representativeness of the obtained samples is unclear.

Fig. 5
figure 5

taken from the original surface by the processive magnifying sampling method

Samples of different sizes

For the simple random sampling method, the sampling results primarily depend on personal judgment, the representativeness of different-sized samples are not quantitatively evaluated. For the equal‑partition sampling method and the processive magnifying sampling method, the heterogeneity of joint roughness is not fully considered, and the representativeness of joint samples is not quantitatively evaluated. To solve the sampling problem, the global search method was introduced to characterize the heterogeneity. The roughness profile of a natural rock joint surface is continuous, while the digitized roughness profile data obtained by measurement devices are available only at a certain interval of horizontal spacing (Kulatilake and Um 1999; Li and Zhang 2015; Yong et al. 2018d). Here, the sampling interval of a digitized profile was labeled SI, the length of the original profile was labeled L, and the lengths of the samples located within the profile (l ≤ L) were labeled l (Fig. 6). The process of evaluating the overall roughness based on the roughness of different samples is, in essence, a process of representing the characteristics of a population by evaluating the selected individual samples. As shown in Fig. 6, there were many joint samples (Sample No.1 to Sample No. (N)) taken from different positions. The total number of individual samples N depends on SI, which can be given as follows:

$$N = \frac{L - l}{{SI}} + 1.$$
(1)
Fig. 6
figure 6

Schematic diagram of the global search method

The joint roughness population of samples taken from different positions can be expressed as follows:

$${\mathbb{R}} = \left\{ {R_{1} ,R_{2} , \cdots ,R_{N} } \right\}.$$
(2)

For instance, for an original profile with a length of L = 1000 mm, the number N of the samples with a length of l = 100 mm can be calculated by Eq. (1). Based on the above sampling process, we can obtain N = 1801 samples under a sampling interval of SI = 0.5 mm. However, as shown in Fig. 3, only ten individual samples can be obtained using the partition sampling method.

As shown in Fig. 6, the samples were extracted from different positions on the original profile by moving an SI distance each time. All the individual samples in the population can be obtained based on this sampling process. The change in the starting point (SP) and the end point (EP) makes each sample independent from the other samples. Sufficient samples can be obtained using the global search method without abandoning any side sections of the original profile. Furthermore, the same sampling result can be obtained from the forward and backward directions of the rock joint profile (Fig. 6a, b).

Generally, heterogeneity refers to a phenomenon in which individual trials have results that are different from each other. Here, the characterization of the heterogeneity is based on the analysis of the roughness differences of all individual samples. Heterogeneity indicates that rock joint roughness measurements vary with position, which is inversely proportional to the similarity between individual samples. The smaller similarity between any two joint samples taken from different positions indicates higher heterogeneity, and vice versa.

Let T = (t1, t2, …, tm) and R = (r1, r2, …, rm) be two m-dimensional vectors. Then, the similarity measures between the two vectors T and R are calculated by the Jaccard similarity measure (Ye 2014) as follows:

$$J(T,R) = \frac{T \cdot R}{{\left\| T \right\|_{2}^{2} + \left\| R \right\|_{2}^{2} - T \cdot R}} = \frac{{\sum\limits_{i = 1}^{m} {t_{i} r_{i} } }}{{\sum\limits_{i = 1}^{m} {t_{i}^{2} + \sum\limits_{i = 1}^{m} {r_{i}^{2} - \sum\limits_{i = 1}^{m} {t_{i} r_{i} } } } }},$$
(3)

where \(T \cdot R = \sum\limits_{i = 1}^{m} {t_{i} r_{i} }\) is the inner product of the vectors \(T\) and \(R\), and \(\left\| T \right\|_{2} = \sqrt {\sum\limits_{i = 1}^{m} {t_{i}^{2} } }\) and \(\left\| R \right\|_{2} = \sqrt {\sum\nolimits_{i = 1}^{m} {r_{i}^{2} } }\) are the Euclidean norms of \(T\) and \(R\), respectively.

Regarding the anisotropy of joint surface roughness, many researchers have suggested that the roughness parameters should be calculated in different directions for full characterization of the 3D surface roughness (Grasselli et al. 2002, Tatone and Grasselli et al. 2012). Let Ri = (ri1, ri2, …, rim) and Rj = (rj1, rj2, …, rjm) (Ri, Rj ∈ \({\mathbb{R}}\), i ≠ j) be two m-dimensional vectors of joint roughness. Here, m is the number of orientations. Ri and Rj represent the roughness individual samples in different orientations of any two samples.

Take the directional roughness \(\theta_{\max }^{*} /[C + 1]_{3D}\) of a sample as an example. The roughness was calculated in directions between 0° and 360° in 5° increments. Thus, the number of the orientations m is 73. Then, the joint roughness vectors can be expressed by Ri = (ri1, ri2, …, ri73). Here, ri1 is the obtained roughness in direction of 0°, ri2 is the obtained roughness in direction of 5°, …, and ri73 is the obtained roughness in direction of 360°.

The joint roughness heterogeneity can be quantified as follows:

$$H = \frac{1}{h}\sum\limits_{{R_{i} ,R_{j} \in {\mathbb{R}}}} {\frac{{\left\| {R_{i} } \right\|_{2}^{2} + \left\| {R_{j} } \right\|_{2}^{2} - R_{i} \cdot R_{j} }}{{R_{i} \cdot R_{j} }}} = \frac{1}{h}\sum\limits_{{R_{i} ,R_{j} \in {\mathbb{R}}}} {\frac{{\sum\limits_{k = 1}^{m} {r_{ik}^{2} + \sum\limits_{k = 1}^{m} {r_{jk}^{2} - \sum\limits_{k = 1}^{m} {r_{ik} r_{jk} } } } }}{{\sum\limits_{i = 1}^{m} {r_{ik} r_{jk} } }}} ,$$
(4)

where \(h = N\left( {N - 1} \right)/2\) is the number of all the combinations of any two individual samples in \({\mathbb{R}}\); \({\mathbb{R}}\) is the population of joint roughness vectors; and rik and rjk indicate the roughness parameters of Ri and Rj in orientation k (1 ≤ k ≤ m). The schematic figure is shown in Fig. 7.

Fig. 7
figure 7

Schematic figure of the quantification of the joint roughness heterogeneity

The determination of the sample number is a significant issue in the sampling survey (Murthy 1967). We suppose that n individual samples are selected from a population of N. Let the joint roughness at different positions be denoted by JRC1, JRC2, and JRCi. On the basis of the partition sampling method, the joint roughness heterogeneity is indicated by the expected value \(\hat{\theta }\), which is given as follows:

$$\hat{\theta }\sim N\left( {\theta ,V(\hat{\theta })} \right).$$
(5)

The population parameter \(\theta\) is given as follows:

$$\theta = \frac{1}{N}\sum\limits_{i = 1}^{N} {JRC_{i} } .$$
(6)

The accuracy of the estimation of population parameter \(\theta\) is usually defined by the absolute error \(d\) or the margin of error r. The margin of error statistically expresses the random sampling error in a survey’s results, and it shows how close the sample’s results are to those when the entire population has been sampled. To ensure the accuracy of the estimation with confidence 1-α, the difference between the expected value \(\hat{\theta }\) and the population parameter \(\theta\) should be within the error, which is expressed as follows:

$$P\left( {\left| {\hat{\theta } - \theta } \right| \le d} \right) = 1 - \alpha ,$$
(7)
$$P\left( {\frac{{\left| {\hat{\theta } - \theta } \right|}}{\theta } \le r} \right) = 1 - \alpha .$$
(8)

Based on the statistical measurement results in previous studies (Du 1999; Tanyas and Ulusay 2013; Ye et al. 2016), the result of the JRC values is an approximate normal distribution when the sample number exceeds 30. Thereafter, we can obtain the following:

$$P\left( {\frac{{\left| {\hat{\theta } - \theta } \right|}}{{\sqrt {V\left( {\hat{\theta }} \right)} }} \le t} \right) = 1 - \alpha ,$$
(9)

where \(V\left( {\hat{\theta }} \right)\) is the sampling variance and \(t\) is the bilateral quantile of the standard normal distribution under a confidence level of 1-α. For instance, \(t\) equals 1.96 when α = 0.05.

The absolute error \(d\) and \(V\left( {\hat{\theta }} \right)\) can be expressed as follows:

$$d = t\sqrt {V\left( {\hat{\theta }} \right)} ,$$
(10)
$$V\left( {\hat{\theta }} \right) = \left( {1 - \frac{n}{N}} \right)\frac{{S^{2} }}{n},$$
(11)

where S is the standard deviation of the population, N is the number of all samples extracted from different locations of a given rock joint, and n is the number of samples obtained by the partition sampling method.

Based on Eqs. (10) and (11), the effective sample number n0 under a confidence level of 1-α is given as follows:

$$n_{0} = \frac{{Nt^{2} S^{2} }}{{Nd^{2} + t^{2} S^{2} }} = \frac{{Nt^{2} S^{2} }}{{Nr^{2} \theta^{2} + t^{2} S^{2} }}.$$
(12)

For instance, we can obtain all individual samples of 400 mm-long samples on the original profile using the global search method. The number of samples can be calculated by Eq. (1), as follows:

$$N = \frac{1000 - 400}{{0.5}} + 1 = 1201.$$

Then, statistical analysis was performed based on the JRC values of the obtained 400-mm-long samples, and the population parameter (\(\theta\) = 7.49) and standard deviation (S = 1.42) were determined. If the confidence level is 95% and the margin of error r is 5%, then the effective sample number is calculated as follows:

$$n_{0} = \frac{{Nt^{2} S^{2} }}{{Nr^{2} \theta^{2} + t^{2} S^{2} }} = \frac{{1201 \times 1.96^{2} \times 1.42^{2} }}{{1201 \times 0.05^{2} \times 7.49^{2} + 1.96^{2} \times 1.42^{2} }} = 52.8.$$

Thus, the sample number needed to statistically describe the roughness should exceed 53. This finding challenges the number of samples obtained by the partition sampling method.

Based on the partition sampling method, n (n = L/l) individual samples can be attained by the equal partition of the joint profile. On the basis of Eq. (12), the back-calculated margins of error attained using different methods are tabulated in Table 1. The margins of error obtained using the simple random sampling method or the processive magnifying sampling method range from 31.3 to 48.9% with a confidence level from 90 to 99%. The partition sampling method has a margin of error in the range from 9.8 to 15.4%. The samples obtained using the partition sampling method are more representative than the samples selected using the simple random sampling method or the processive magnifying sampling method, but estimation errors consistently arise due to inferences about the population based on observations of part of it. The margin of error is commonly taken as 5% for engineering practice (Yong et al. 2018c). The obtained margin of error using the partition sampling method is comparatively larger. Thus, it is inadequate to use the expected value \(\hat{\theta }\) obtained using the partition sampling method to represent the population parameter \(\theta\). The global search method can be applied to analyze the roughness heterogeneity in the population to solve the sampling problem.

Table 1 Back-calculated margins of error based on different sampling methods

4 Rock Joint Sample and Roughness Measurements

4.1 Joint Selection and Description

The Heshangnong quarry is located at Qingshi Town, southeast of Changshan County, Zhejiang Province, China, approximately 213.5 km from Hangzhou city as shown in Fig. 8a. The exploitation of this quarry requires a pit with a length of 87 m, a width of 59 m and a maximum height of 79 m. In this pit, the overburden mainly consists of calcareous slate (Fig. 8b), which originated from Ordovician argillaceous limestone under the condition of light metamorphism. The stability of this open pit is controlled by the slate foliation, which generally dips approximately 55° to the NW. The grayish-green slate rock wall is foliated, very fine grained, and formed by the metamorphosis of intermediate tuff. The very distinct, continuous foliation planes in the overburden rock are oriented with strikes approximately parallel to the pit walls and dips towards the bottom of the pit (Fig. 8b).

Fig. 8
figure 8

Sites selected for the investigation. a Location map of the study site; b view of the structurally controlled open-pit slope

The roughness parameters of rock joints are usually scale-dependent, but the scale dependency of joint roughness is limited to a certain size, defined as the stationarity threshold. The change in the roughness parameter is not visible for sample sizes greater than the stationarity threshold. To determine the morphological and mechanical properties of rock joints at laboratory and field scales, the size of the samples should be equal to or larger than the stationary threshold (Fardin 2007). Following the stationary threshold size mentioned in previous studies (Fardin et al. 2004; Tatone and Grasselli 2013), a sample with an overall area of 1100 × 1100 mm2 (Fig. 9a) was sawed from the slate rock and transported to the laboratory. A study area with a size of 1000 × 1000 mm2 was obtained from the center to avoid damage had occurred while in transit to the laboratory.

Fig. 9
figure 9

Generation of the joint surface digitization. a Scanning of the joint surface. b Digitized joint surface. c Locations of the digitized profiles

4.2 Digitization of the Joint Surface

A 3D laser scanning system, MetraScan 750 (Fig. 9a), with a maximum accuracy of 0.030 mm, was used to measure the geometry of the joint surface, and its main components include a scanner, C-Track cameras, a controller, and a computer. The surface acquisition was solved by observing the laser lines projected on the rock joint surface. As the laser swept over the surface by the scanner, the scanner measured the 3D coordinates of the sample surface using seven laser crosses, and the data were registered depending on the triangulated position. The final 3D surface model of the study area was obtained following point cloud data processing using the scanner software (Fig. 9b). We obtained 51 cross-sectional profiles along the potential shearing direction (Fig. 9c), namely, P01, P02, …, and P51. These profiles were digitized at a sampling interval (SI) of 0.5 mm, as this value is often applied in previous studies (Li and Zhang 2015; Tatone and Grasselli 2010; Yong et al. 2018c).

4.3 Quantification of Joint Roughness

In this study, empirical, statistical, and fractal roughness parameters were used to describe the geometric irregularities of rock joints.

First, JRC, the most commonly used empirical roughness parameter, was taken to quantify the joint roughness. For field estimation purposes, the JRC values can be approximately determined based on the relations between the JRC and the maximum profile amplitude (Amax) measured over a sample length (L) (Barton 1981; Palmström 1995). The relation between the JRC and the undulation factor u (u = Amax/L) can be used to evaluate the scale effects on the JRC. Yong et al. (2018a) proposed a programmed method for determining the Amax of digitized joint profiles. As shown in Fig. 10, the main process of determining Amax is introduced as follows:

Fig. 10
figure 10

Flowchart illustrating the steps involved in determining Amax

Step 1: Arbitrarily select two points from the digitized profile (xp,yp), (xq,yq).

Step 2: Obtain the linear equation of the line passing through (xp,yp), (xq,yq) as follows:

$$\overline{y} = \frac{{y_{p} - y_{q} }}{{x_{p} - x_{q} }}(x - x_{p} ) + y_{p} .$$
(13)

Step 3: Select another point (xi,yi) on the digitized profile, where xi ∈(xp, xq). Then, determine the asperity amplitude by calculating the normal distance from this point to the line as follows:

$$A = \frac{{\left| {\frac{{y_{p} - y_{q} }}{{x_{p} - x_{q} }}(x_{i} - x_{p} ) + y_{p} - y_{i} } \right|}}{{\sqrt {1 + \left( {\frac{{y_{p} - y_{q} }}{{x_{p} - x_{q} }}} \right)^{2} } }}.$$
(14)

Step 4: According to steps (1) to (3), an iterative process was used to determine the asperity amplitudes, by changing the points (xp, yp), (xq, yq) and (x0, y0), which took into account all possible combinations of these three points from the digitized profile.

Step 5: The maximum asperity amplitude Amax was determined by finding the maximum value of A.

Based on the results, the joint roughness coefficient based on Barton’s straight edge method (Barton 1981; Yong et al. 2018a) is as follows:

$$JRC_{n} = 400\frac{{A_{\max } }}{{L_{n} }}.$$
(15)

In addition, four widely accepted statistical parameters were selected to investigate the joint roughness heterogeneity. The following three parameters are used to quantify of 2D joint roughness: the first-derivative root-mean-square \(Z_{2}\) (Tse and Cruden 1979), the roughness profile index or profile sinuosity \(R_{p}\) (Maerz et al. 1990), and the modified root-mean-square Z2 (Zhang et al. 2017). The first-derivative root-mean-square \(Z_{2}\) and the roughness profile index or profile sinuosity \(R_{p}\) are common roughness parameters, and they can be determined by the following equations:

$$Z_{2} = \sqrt {\frac{1}{L}\int_{x = 0}^{x = L} {\left( {\frac{dy}{{dx}}} \right)^{2} dx} } \approx \sqrt {\frac{1}{N}\sum\limits_{i = 1}^{N - 1} {\left[ {\frac{{\left( {y_{i + 1} - y_{i} } \right)}}{{(x_{i + 1} - x_{i} )}}} \right]}^{2} }$$
(16)
$$R_{p} = \frac{{L_{t} }}{L} = \frac{{\sum {_{i = 1}^{N - 1} \sqrt {\left( {x_{i + 1} - x_{i} } \right)^{2} + \left( {y_{i + 1} - y_{i} } \right)^{2} } } }}{L},$$
(17)

where x and y are the horizontal and vertical coordinates of the points along a profile; (xi+1, yi+1) and (xi, yi) represent the adjacent coordinates of a profile; N is the number of measured points over a profile; Lt is the true length of the profile, which refers to the sum of the lengths of the lines between two adjacent points; and L is the length of the profile projected on the x-axis (see Fig. 11a).

Fig. 11
figure 11

Variation of inclination angles along a rock joint profile (Zhang et al. 2014)

Zhang et al. (2014) found that the results of shear box tests in forward and reverse directions were different, which indicated that roughness was also different in forward and reverse directions. To account for the directional dependencies, they classified the dilation angles into positive and negative angles (see Fig. 11b) and suggested that only positive angles should be considered in the process of determining joint roughness. Based on Eq. (16), Zhang et al. (2014) defined the modified root-mean-square Z2 as follows:

$$Z^{\prime}_{2} = \sqrt {\frac{1}{L}\int_{x = 0}^{x = L} {\left( {\max \left( {{0,}\frac{dy}{{dx}}} \right)} \right)^{2} dx} } = \left[ {\sum\limits_{i = 1}^{N - 1} {\frac{{(\max (0,y_{i + 1} - y_{i} ))^{2} }}{{(x_{i + 1} - x_{i} )L}}} } \right]^{1/2} .$$
(18)

In this study, the roughness metric \(\theta_{\max }^{*} /[C + 1]_{3D}\) (Grasselli et al. 2002) was used to characterize the 3D joint roughness. Grasselli (2001) found that only the parts of the joint surface that face the shear direction and are steeper than a threshold inclination provide shear resistance. Then, he developed the roughness metric \(\theta_{\max }^{*} /[C + 1]_{3D}\) based on an angular threshold concept, and this parameter was initially developed to identify potential contact areas during direct shear testing of rock joints.

As shown in Fig. 12a, the joint surface is discretized into adjacent triangles. Each triangle orientation is uniquely identified by its azimuth angle α and dip angle θ. Azimuth angle α is the angle between the true dip vector d projected on the shear plane and the shear vector t, measured clockwise from t. Dip angle θ is the angle between the shear plane and the triangle. As presented in detail by Grasselli et al. (2002) and Grasselli (2006), the apparent dip angle θ* describes the apparent inclination of each triangle with respect to the shear direction, and it is obtained by projecting the dip angle along the vertical plane which contains the shear direction, as follows:

$$\tan \theta^{*} {\text{ = - tan}}\theta \cdot {\text{cos}}\alpha .$$
(19)
Fig. 12
figure 12

Schematic figure showing the definition of the roughness metrics \(\theta_{\max }^{*} /[C + 1]_{3D}\) (Tatone and Grasselli 2009). a Sketch of the apparent dip angle θ*. b Schematic diagram illustrating the use of the angular threshold concept to determine Aθ. c Cumulative distribution of Aθ as a function of θ*

For each surface, it is possible to calculate the area of the joint that has an apparent dip angle equal to or greater than a chosen threshold dip angle, which represents the area in contact or damaged during shearing (Grasselli, 2006). As shown in Fig. 12b, the areas that are steeper than a certain value of θ* are highlighted on the entire surface. The ratio of the areas steeper than a threshold dip angle to the area of the entire surface is denoted by Aθ∗. For example, the areas that are steeper than a threshold dip angle (θ* = 0) are highlighted in blue color, and the ratio \(A_{\theta *}\) of the highlighted area to the entire surface is 0.540. As shown in Fig. 12b, when the values of θ* equal to 0°, 5°, 10°, 20°, and 30° are considered, the corresponding values of \(A_{\theta *}\) are 0.540, 0.399, 0.271, 0.092, and 0.022, respectively.

By varying the θ* from 0° to the maximum apparent dip angle \(\theta_{\max }^{*}\), the ratios \(A_{\theta *}\) were obtained. Figure 12c shows the relationship between \(A_{\theta *}\) and θ*. It is possible to plot the variation of θ* as a function of the threshold dip angle, and their relationship can be expressed in terms of Tatone and Grasselli (2010)’s equation:

$$A_{\theta *} = A_{0} \left( {\frac{{\theta_{\max }^{*} - \theta^{*} }}{{\theta_{\max }^{*} }}} \right)^{C} ,$$
(20)

where A0 is the total potential contact area ratio when θ* = 0°; C is the dimensionless empirical fitting parameter calculated via non-linear least-squares regression.

Then, the roughness metric \(\theta_{\max }^{*} /[C + 1]_{3D}\) can be calculated based on the obtained parameters \(\theta_{\max }^{*}\) and \(C\). In addition, Tatone and Grasselli (2010) expand the application of the roughness metric and proposed a new parameter \(\theta_{\max }^{*} /[C + 1]_{2D}\) for the characterization of 2D roughness. In this work, this parameter was also used to quantify the roughness of joint profiles.

Furthermore, the fractal dimensions D of a rock joint profile were calculated using the compass-walking method, which is also a widely accepted parameter for quantifying the roughness of a natural rock joint profile (Li and Zhang 2015).

5 Measurement Results and Analysis

5.1 Heterogeneity of Joint Roughness by Different Roughness Parameters

The original joint profile P01 is applied here as an example (Fig. 13a), the global search method was applied to investigate the roughness heterogeneity in situation A. Profile P01 is 1000 mm long. Based on the partition sampling method, this profile can be equally separated into ten 100 mm subsections (Fig. 13b). The number of the samples obtained is n = 10. Based on the global search method, the samples were taken from different positions on the original profile by moving an SI distance each time (Fig. 13c). The total number N of subsections is 1801 based on Eq. (1).

Fig. 13
figure 13

Schematic diagrams of the sampling results. a The original joint profile P01. b Partition sampling method. c Global search method

The joint roughness parameters of all the individual samples in the roughness population obtained by the global search method are displayed in the histograms in Fig. 14. The frequency histograms of the roughness parameters indicate the joint roughness heterogeneity. We can evaluate the joint roughness heterogeneity by Eq. (4). In this study, the roughness of the joint was quantified by empirical, statistical, and fractal roughness parameters. To ease the comparison between the roughness heterogeneity by different parameters, all the roughness parameters were normalized into a range between 1 and 2 before calculating the roughness heterogeneity. The normalized values \(r_{i}\) of the roughness parameters are obtained using the following equation:

$$r_{i} = \frac{{J_{i}^{ * } - J_{{{\text{min}}}}^{*} }}{{J_{\max }^{*} - J_{{{\text{min}}}}^{*} }} + 1,$$
(21)

where \(J_{i}^{ * }\) is a roughness parameter which can be JRC, Z2, Rp, Z2’, \(\theta_{\max }^{*} /[C + 1]_{2D}\), and D. \(J_{\min }^{*}\) and \(J_{\max }^{*}\) are the minimum and maximum values of the roughness parameters, respectively.

Fig. 14
figure 14

Histograms of the heterogeneity of joint roughness using different roughness parameters. a JRC, b Z2, c Rp, d Z2’, e \(\theta_{\max }^{*} /[C + 1]_{2D}\), f D

The ratio of heterogeneity (H ≥ 1) indicates the overall difference between the roughness of the samples taken from different positions. When the ratio equals 1, it reflects that the joint surface roughness is homogeneous. Larger ratios indicate more noticeable differences in the roughness of samples collected from various positions. According to the obtained normalized values of different roughness parameters, the heterogeneity ratios H were calculated using Eq. (4). The heterogeneity ratios H for JRC, Z2, Rp, Z2’, \(\theta_{\max }^{*} /[C + 1]_{2D}\), and D were 1.051, 1.050, 1.050, 1.059, 1.047, and 1.046, respectively. This result indicates that different roughness estimation methods produce close evaluations of the heterogeneity.

5.2 Error in Joint Roughness Heterogeneity Estimated by the Partition Sampling Method

In previous studies, JRC is the most commonly used parameter for quantifying the joint roughness, because it provides insights into the shear strength and deformation behavior of rock joints (Morelli 2014). Here, it was taken to study the error in the joint roughness heterogeneity estimate by the partition sampling method.

Figure 15a is the frequency histogram of the JRC values of the samples obtained using the global search method. A normal distribution is shown here with a mean value \(\theta\) of 13.34 and a standard deviation S of 3.40. Taking into account the fact that all samples in the population can be obtained by the global search method, \(\theta\) is the population parameter. However, the JRC values of the samples obtained using the partition sampling method do not follow a normal distribution (Fig. 15b). Compared with the statistical results obtained by the global search method, the expected value \(\hat{\theta }\) obtained using the partition sampling method is close to the population parameter \(\theta\), but the standard deviation \(\hat{S}\) obtained using the partition sampling method is distinctly different from S. In addition, the roughness distribution obtained using the partition sampling method (Fig. 15b) is inconsistent with the result obtained using the global search method (Fig. 15a). For example, the JRC values in the ranges from 12 to 14 and from 18 to 20 are missing in Fig. 15b. However, the JRC values range from 12 to 14 with highest relative frequency (Fig. 15a). Therefore, the distinct differences in the standard deviation and roughness distribution indicate that the statistical results based on the partition sampling method may introduce estimation errors.

Fig. 15
figure 15

Comparison of the JRC values of 100 mm joint samples. a Histograms that show the distribution of the JRC values using the global search method. b Histograms that show the distribution of the JRC values using partition sampling method

In addition, we can evaluate whether the sample number is sufficient based on Eq. (12). Under a confidence level of 95% and a margin of error of 10%, the effective sample number is calculated using Eq. (12) as follows:

$$n_{0} = \frac{{Nt^{2} S^{2} }}{{Nr^{2} \theta^{2} + t^{2} S^{2} }} = \frac{{1801 \times 1.96^{2} \times 3.40^{2} }}{{1801 \times 0.1^{2} \times 13.34^{2} + 1.96^{2} \times 3.40^{2} }} \approx 25.$$

Only ten samples can be obtained by the partition sampling method. This number is less than n0. Thus, the sample number suggested by the partition sampling method is insufficient. Furthermore, only one sample could be obtained based on the simple random sampling method and the processive magnifying sampling method. It is smaller than the sample number by the partition sampling method, thus it is also insufficient. In statistics, the margin of error is a statistic expressing the amount of random sampling error in the results of a survey. It indicates how far the results deviate from the real population value. We can calculate the margin of error by Eq. (12) for a confidence level of 95%. The back-calculated margins for the partition sampling method, the simple random sampling method and the processive magnifying sampling method are 15.75%, 49.9% and 49.9%, respectively. The margin of error obtained using these traditional sampling methods is greater than 10%. Therefore, the roughness parameters of the samples obtained by these conventional sampling methods are not representative of the overall roughness.

For each joint profile with length of 1000 mm, we can obtain ten joint samples with the length of 100 mm or two samples with the length of 500 mm using the partition sampling method. As previously mentioned, these ten samples obtained by the partition sampling method are insufficient for the statistical analysis of the roughness of the 100 mm joint samples. The 500 mm samples were obtained from the profiles (P01, P02, …, and P51) in Fig. 9 using the partition sampling method. First, the joint samples were obtained from the profiles (P01, P02, …, and P51) in Fig. 9 using the partition sampling method. Then, the expected value \(\hat{\theta }\) and standard deviations \(\hat{S}\) of the samples taken from each profile were calculated. Figure 16 illustrates that the mean values and standard deviations present apparent differences between the values obtained by the global search method and the partition sampling method. For each profile, different samples can be extracted from different locations using the global search method, and the total number of the samples is 1001 based on Eq. (1). Based on the global search method, the heterogeneity of rock joints was characterized by the roughness of all individual samples. However, only two samples can be extracted from each profile using the partition sampling method. The heterogeneity estimated by the partition sampling method may be biased because of the insufficient samples. As shown in Fig. 16a, the expected values (\(\hat{\theta }\)) of the samples obtained from the profiles P01 to P31 are underestimated, but the \(\hat{\theta }\) values of the samples obtained from the profiles P32 to P51 are overestimated. The largest relative error between \(\hat{\theta }\) and \(\theta\) is 32.59%. Figure 16b shows a comparison of the standard deviations of the JRC values based on the partition sampling method and the global search method. For the partition sampling method, the standard deviation \(\hat{S}\) is a measure of the variation between the roughness of the two samples, which cannot reflect the variations between the roughness of samples taken from various locations. The differences between \(\hat{S}\) and \(S\) are greater than the difference between \(\hat{\theta }\) and \(\theta\), and the largest relative error between \(\hat{S}\) and \(S\) is 236.34%. Therefore, the expected values \(\hat{\theta }\) of 500 mm joint samples by the partition sampling method also failed to represent the overall roughness.

Fig. 16
figure 16

Comparison of the JRC values of 500-mm-long joint samples based on partition sampling method and global search method. a The mean values and b the standard deviations

The difference between the population value \(\theta\) by the global search method and the expected value \(\hat{\theta }\) by the partition sampling method can be quantified by the absolute values of the relative errors r’, as follows:

$$r^{\prime} = \frac{{\left| {\theta - \hat{\theta }} \right|}}{\theta } \times 100\% .$$
(22)

As shown in Fig. 17a, the relative errors r’ of the joint profiles (P01, P02,…, and P51) were calculated. It was validated that the differences between the results by the global search method and the partition sampling method not only vary from fracture to fracture but also vary with scale. The distribution of the relative error r’ can be divided into two sections. For the first section, the length of the samples does not exceed 500 mm, the maximum value of the relative error r’ increases with the sample length, and the overall maximum value is 32.59%. For the second section, the sample length is larger than 500 mm, and the maximum value of the relative error r’ decreases with the sample length. In the first section, the increasing trend of the maximum value of the relative error r’ is due to the reduction in the number of samples of each length (Fig. 17b). For the 100-mm-long samples, 10 samples can be obtained using the partition sampling method, while only 2 samples can be obtained for the 400-mm- and 500-mm-long samples. In the second section, the decreasing trend of the relative error r’ is due to the changes in the abandoned part generated by the partition sampling method. For a 600 mm sample, the length of the abandoned part is 400 mm, so the percentage of the abandoned length λ is 40%. The roughness of the obtained 600 mm sample can only represent the roughness characteristics of the joint in the first 600 mm of length. Then, the percentage of the abandoned length λ decreases as the sample length increases (Fig. 17c), which contributes to the decreasing trend of the relative error r’.

Fig. 17
figure 17

Comparison between the population value \(\theta\) by the global search method and the expected value \(\hat{\theta }\) by the partition sampling method. a Variation in the relative errors with the length of the joint samples. b The change in the sample number. c The change in the percentage of the neglected length

5.3 Applications for Searching Representative Samples

Surveying the shear behavior requires an understanding of all the potential sampling errors that can bias the result. Roughness heterogeneity as an efficient indicator of shear behavior may be used to explore representative test samples. Yong et al. (2018b) suggested that a representative sample can be obtained by finding the sample whose roughness is close to the maximum-likelihood estimation of the roughness probability distributions. Here, the representative samples were obtained based on the roughness probability distributions of the joint samples with the same size. The representativeness of each sample can be evaluated by the following:

$$\eta_{i} = \frac{{\left| {JRC_{i} - \theta } \right|}}{\theta } \times 100\% ,$$
(23)

where ηi is the representativeness coefficient; JRCi denotes the JRC of joint sample i; and θ is the mean value of the roughness of all samples obtained from different positions using the global search method.

The most representative sample of each size can be determined when ηi reaches the minimum value. Thereafter, the position of the sample can be determined on the profile, which can be selected for the scale effect study.

By taking the profile P01 as an example, the JRC values of the samples with the lengths of 100 mm to 1000 mm were obtained using the global search method. For the samples of 100 mm length, the mean value \(\theta\) of the roughness of all samples obtained from different positions is 13.34. We calculated the representativeness of each sample with the length of 100 mm using Eq. (23) and found that η of joint sample taken from 779.5 to 879.5 mm on P01 reached the minimum value. That is, the location of the representative sample of a length of 100 mm ranges from 779.5 to 879.5 mm in the x-axial direction. Following the same process, we can determine the location of the representative samples at the lengths of 200 mm to 1000 mm. The positions of the representative samples of different sizes are shown in Fig. 18.

Fig. 18
figure 18

Sampling locations of representative samples based on the global search method

6 Heterogeneity of 3D Joint Roughness

Heterogeneity exists in both 2D profiles and 3D surface topographies (Du 1994). The heterogeneity of joint roughness was analyzed on the basis of the joint profiles in the aforementioned sections. Here, the heterogeneity of the 3D joint roughness was investigated based on joint samples with sizes of 100 × 100 mm2 obtained from the digitized joint surface in Fig. 9.

Based on the global search method, independent roughness samples can be obtained by changing their positions, and a population of the roughness of all samples with a number of (N = 1801)2 can be obtained. Although sampling error exists when using the partition sampling method, the error usually decreases with the sample number. A maximum error of 11.8% was observed in the heterogeneity of 2D profiles with sizes of 100 mm. The digitized joint surface can be divided into 100 samples using the partition sampling method (Fig. 19). The heterogeneity of 3D surfaces was approximately estimated based on the partition sampling method. As shown in Fig. 19, there is a distinct difference in the surface roughness of two adjacent samples S1-1 and S1-2, as the topography of sample S1-2 fluctuates more than that of sample S1-1. Samples S1-1 and S1-10 are obtained from the corners of the joint surface. The surface of sample S1–10 is comparatively smoother than that of sample S1-1. The surface roughness of sample S5-5 is slightly rough, but the areas that face the shear direction are completely different from those of samples S1–1 and S1–10. It is well known that surface roughness is one of the main factors affecting the shear strength of rock joints. Thus, the shear behaviors of samples S1–1, S1–10, and S5–5 are probably different.

Fig. 19
figure 19

3D joint roughness samples with sizes of 100 × 100 mm2

As mentioned in Sect. 2, there are three situations of exposed potential sliding planes encountered in field surveys. The following two groups of 3D joint roughness samples were used to show the difference of heterogeneity in different situations: samples S1-1, S2-1… S10-1 obtained in the direction perpendicular to the shear direction were taken as an example to show the heterogeneity of the samples in situation A, and samples S1-1, S1-2…S1-10 obtained in the direction perpendicular to the shear direction were used to illustrate the heterogeneity in situation B. The former reflects the roughness diversity spread along the sliding direction, and the latter reflects the diversity of surface roughness in the direction perpendicular to the sliding direction. The surface topographies of these samples were displayed in a contour plot to visualize the anisotropy of surface roughness. The surface topographies of the samples in Fig. 20a had greater variance. The differences in the topographies of each sample in Fig. 20b are less apparent. The values of the directional roughness \(\theta_{\max }^{*} /[C + 1]_{3D}\) change with the locations of the joint samples in field situations A and B. In field situation A, the directions of the maximum and minimum values of \(\theta_{\max }^{*} /[C + 1]_{3D}\) change with the locations of the joint samples (Fig. 20c), but the directions of the maximum and minimum values are almost constant in field situation B (Fig. 20d). Using Eq. (4), the heterogeneity ratio H in situation A is 1.046, and the result for situation B is 1.039. Thus, the heterogeneity of the 3D joint roughness is more apparent in situation A than in situation B.

Fig. 20
figure 20

Comparisons between the heterogeneities of the samples. a Contour plots of the surface topographies of the samples under field condition A. b Contour plots of surface topographies of the samples under field condition B. c Polar plots of the 3D directional roughness surface of the samples in situation A. d Polar plots of the 3D directional roughness surface of the samples in situation B

The values of \(\theta_{\max }^{*} /[C + 1]_{3D}\) along the sliding direction of all samples are shown in Fig. 21, which vary from 5.24 to 10.69. These data follow a normal distribution with a mean value of 7.98 and a standard deviation of 1.12. Here, the representative samples were achieved based on the \(\theta_{\max }^{*} /[C + 1]_{3D}\) probability distributions of the joint samples. Sample S4-3 was found to be a representative sample, whose roughness metric \(\theta_{\max }^{*} /[C + 1]_{3D}\) is 8.08.

Fig. 21
figure 21

Histogram of the 3D joint roughness along the sliding direction

7 Conclusions

Roughness heterogeneity refers to the fact that joint roughness measurements vary with position, and it exists in both two-dimensional (2D) profiles and three-dimensional (3D) surface topographies. Little attention has been given to the influence of heterogeneity on joint roughness estimation in previous studies, and the problem regarding the inaccuracy of joint roughness estimation based on inadequate samples remains unsolved. According to the statistical analysis of the effective sample number, the roughness samples in conventional methods are normally insufficient to represent the overall roughness and the estimation error consistently arises in the characterization of the joint roughness heterogeneity.

A new method, namely the global search method, was proposed to characterize the heterogeneity of rock joints based on the analysis of the roughness of all individual samples. This method allows all independent samples within the populations to be obtained by changing the starting and end points, without abandoning any side sections of the original profile. This method is convenient for acquiring a sufficient number of samples of various sizes and helpful for improving the determination of the scale effect. It was verified that the expected value obtained from conventional methods cannot represent the overall roughness of the joint profile. The proposed method is capable of systematically characterizing the heterogeneity of the roughness of rock joints of various sizes.

In the case study on the natural slate joint, the joint roughness heterogeneity was systematically investigated based on both 2D profiles and 3D surface topographies. The error in the estimation of joint roughness heterogeneity using the partition sampling method was analyzed. It illustrates that the maximum relative error between the population parameter and the expected value first increased with the length of the joint sample and then decreased when the sample length was exceeding the half of the length of the original profile. The characterization of joint roughness heterogeneity is applied for the exploration of representative samples. Representative samples were obtained by finding the sample whose roughness was closest to the population parameter of all individual samples.