Study on Effect of Scaling and Anisotropy on Roughness of Natural Fractured Rock Surfaces

Abstract

In the modified formula for the cubic law, the roughness correction coefficient C serves as a parameter that indicates the characteristics of the roughness of a rock’s surface. In this study, surface data for natural rock samples were acquired through high-precision 3D scanning and combined with publicly accessible CT scan data on rough rock fractures to generate spatial coordinates. The roughness correction coefficient C was calculated and analyzed using both formulaic and numerical methods. The analysis revealed significant effects of scale on the roughness correction coefficient for rock fractures within a size range of 5 to 10 cm (determined based on the actual sample size), with tensile fractures demonstrating greater variations compared to shear fractures. When calculating the roughness correction coefficient on the same fracture surface in different directions, significant directional effects were observed. Furthermore, elliptical fitting demonstrated favorable results. The conclusion drawn was that the roughness correction coefficient for fracture surfaces can be represented effectively using a tensor form, thereby simplifying the expression of directionality. Calculating and analyzing the fractal dimension of a rough surface further confirmed the existence of effects of size on roughness.

1. Introduction

Roughness affects the rock geometry, mechanics, and transport characteristics within fault zones [1]. In studying fractured seepage, the roughness of a fractured rock surface represents the degree of fracture in the surface structure and explains the mechanics and seepage characteristics. To date, three main types of methods have been used to describe fractured roughness: the raised height characterization method, the joint roughness coefficient (JRC) characterization method, and the fractional dimension characterization method [2]. In recent research, Miao et al. [3] studied the impact of fracture roughness on solute transport, examining the flow rate and velocity against the fractal dimension. Meanwhile, Gan et al. [4] leveraged the JRC to quantify the effect of roughness on seepage and derived expressions for this relationship.

A roughness correction coefficient C is added to formulas for the cubic method to reflect the influence of roughness. From the “Barton Ten Rules” [5] put forward by N. Barton and V. Choubey to the fraction dimension representation method developed later [6,7], many researchers have formulated their own revised cubic laws and found different roughness correction coefficient values. The basic elements of a fractured rock mass network are a single fracture and its various parameters, such as its roughness, filling, and degree of fractured surface contact, which directly affect the fractured network and thus the permeability of the whole fractured rock mass [8,9]. Therefore, as a coefficient that describes the fractured roughness, the roughness correction coefficient is indispensable within a corrected formula for the cubic law. Many scholars have devised formulas for calculating roughness [2,10]. The roughness correction coefficient C varies within these but tends to be constant.

Lang et al. [11] used a tensor form to represent the equivalent permeability coefficient for a fractured rock mass, obtaining two-dimensional and three-dimensional tensor representations of permeability. Equally, by studying a seepage model, Liu Weiqun et al. [12] found that the seepage rates differ in fractured shale in three directions and are about three times higher in the vertical or horizontal direction. These findings led us to wonder whether the roughness correction coefficient should also be directional in a corrected formula for the cubic law and whether it can be expressed in a corresponding tensor form if so.

To this end, this study used a high-precision 3D scanner to scan the surfaces of natural rock samples to obtain accurate rock surface data. In combining these data with high-precision data from CT scans of rough fractured rocks published online, a rough fractured model was generated, and after processing, a single fractured roughness correction coefficient was calculated. The rough fractured rock was divided into top and bottom surfaces. After processing the data, the spatial position coordinates for each point were read, and the gap width was calculated. The proportions of the largest square that could be taken from a single slit surface were recorded at a scale of 100% and then continuously reduced (from 100% to 10% in 10 groups), and the gap width and the spatial coordinate data were intercepted and read to study the effect of the sample size on the roughness correction coefficient. Rotating samples of the same size at a given angle, we intercepted the spatial coordinate data for different angle gap widths and fractured surfaces and calculated the change in the roughness correction coefficient in different directions to study the directional characteristics. Then, we put the results of calculating the correction coefficient at each angle into the polar coordinate system and applied elliptic fitting to observe the overall characteristics of directionality. According to the laws of changes in direction and scale, a relevant law for fractured roughness was obtained.

2. Data

Compared with the artificial samples used in certain studies, in natural rock fractures, the degree of roughness is more realistic due to the influence of structural forces and weathering. Therefore, the fractured rock surfaces analyzed in this study were all naturally formed, some of which were sampled from the Gejiu Gaosong mineral field in Yunnan Province, while some of were derived from data on natural fractured rocks published online [13,14].

The study area is characterized by the development of folds and faults, with the fault permeability significantly impacting the seepage of groundwater. The surface in this region is mainly covered by sedimentary rocks. Based on its exposure and the sedimentary characteristics of the strata, the geological formations in this area can be roughly divided into two types: bedrock mountainous areas and basin areas. According to the structural laws and characteristics, the area can be primarily divided into two major tectonic belts, namely east–west and north–south orientations, with the north–south tectonic belt being particularly prominent.

2.1. Data Sources

In this study, we employed the high-precision 3D scanner TY-3Dscan to scan the rock fractures, as illustrated in Figure 1a. We collected the scanned samples from natural fractured rock surfaces during a geological survey of the Gaosong minefield in Gejiu, Yunnan Province. The rock fractures observed in this study were shear-type fractures in structures predominantly composed of limestone. We selected three samples for this study, namely S1, S2, and S3. Figure 1b depicts photographs captured during the on-site survey, while Figure 2(a₁–c₃) display images of the samples collected.

To enhance the reliability of the data and reduce the impact of randomness on the sample selection, we selected two certified and highly credible sets of fracture data made publicly available online in order to perform identical computational analyses. In the original representations of the computed tomography (CT) data presented in Figure 2(d₁–e₃), two types of materials exhibiting tensile fractures are shown. Material “Figure 2d₁” corresponds to tuff and is denoted as S4, whereas material “Figure 2e₁” represents granite and is designated as S5. The fracture divides the superior and inferior surfaces, and the gap widths are computed through subsequent subtraction of the data, with the corresponding spatial coordinates derived. A rectangular Data Extraction Frame that captured the largest quantity of original data from both surfaces was implemented to simplify the calculation and analysis. Figure 2(d₂, d₃, e₂ and e₃) shows three-dimensional views of the gap width and the bottom surface.

Normal and tangential displacement methods [15] were applied to the scanned fractured surfaces to obtain data on the top and bottom surfaces, calculate the gap width, and export the data on the bottom surfaces and the spatial coordinates of the data on the gap width. In order to facilitate the calculations and analyses and retain as much original data as possible, the rectangular interception frame captured the gap width and the bottom surface so that the maximum amount of original data on the top and bottom surfaces were accommodated.

2.2. Data Processing

Rectangular intercept boxes can be either rectangular or square. Given that a difference in side length has the potential to interfere with the experimental results, using a square intercept box is recommended, as this is considered more conducive to analysis of the results. The procedure for intercepting the data is outlined below.

To explore the anisotropy of the roughness correction coefficient, we used a square intercept box to intercept the fracture surface data and the gap width data, respectively, from the spatial coordinate data for each sample described above.

Based on the rectangular frame defined in Section 2.1, we established the centroid of this frame as the reference point and conducted a series of interception operations using the maximum size as the basis, with the interception ratios set to 10%, 20%, 30%, 50%, 60%, 70%, 80%, 90%, and 100% of the frame’s dimensions. The aim was to investigate how variations in size impact the roughness correction coefficient. Subsequently, the interception frame was kept at a constant size while being rotated counterclockwise. During the rotation, a set of data samples containing information on the fracture surface and gap width was intercepted every 5 degrees until the rectangular frame had rotated a full 360 degrees, signifying successful completion of the data collection phase for the interception model. Figure 3 is a clear schematic diagram illustrating this operational procedure.

During the data organization phase, we found that, using the method mentioned, each sample generated a total of 72 specific data points. These points could be further consolidated into 36 sets of sample data, with each set comprising one top surface and one bottom surface, resulting in a total of 720 independent data points in this study. However, given the symmetry between the data from 0° to 180° and from 180° to 360°, we limited our calculations to the 0° to 180° range. Subsequently, leveraging the principle of symmetry, we expanded this range of data and plotted them onto a comprehensive chart spanning the full 360° spectrum. Ultimately, we utilized a total of 1080 sets of data (taking into account the efficient utilization of symmetry) for the subsequent analysis.

3. Methodology

In discussing the seepage laws for a single fracture, there are generally three ways to derive them: theoretical research [16], experimental research [17], and numerical research [18]. Using data on natural fractured surfaces and CT scan data on rough rocks as the samples [13,14], we analyzed spatial coordinate data on the fracture surfaces on the basis of theoretical research (the formula method) and numerical study (the numerical simulation method) to determine whether the corrected coefficients reflect effects of directionality and scaling on roughness.

3.1. Calculating the Roughness Correction Coefficient—The Formula Method

The cubic law is generally adopted in studying seepage in fractured rock. A plethora of researchers have derived their own modified cubic law or local cubic law and added a roughness correction coefficient C, such as that given by Zhang Youtian [19]:

$$\mathrm{q}={\displaystyle \frac{\mathrm{g}{a}^3}{12\mathsf{\nu }\mathrm{C}}}\mathrm{J}$$

(1)

where q is the fractured flow, g is the gravity acceleration, a is the fractured width (mechanical gap width), ν is the motion viscosity coefficient of water, and J is the hydraulic gradient.

We can write the cubic law above in the form of Darcy’s law:

$$\mathrm{q}=-k_{f}a\mathrm{J}$$

(2)

where ${ k}_f$ is the hydraulic conduction coefficient of the fractured rock, $k_f={\displaystyle \frac{g{a}^2}{12vC}}$,

and

$$k_f={\left(k_{in}\right)}_f{\displaystyle \frac{\rho g}{\mu }}$$

(3)

with

$${\left(k_{in}\right)}_f={\displaystyle \frac{a^2}{12C}}$$

(4)

As observed in the equation above, the intrinsic hydraulic conduction coefficient of the fracture, ${\left(k_{in}\right)}_f$, is affected by the gap width and by the roughness correction coefficient. Upon incorporating the hydraulic gap width e, and with

$${\left(k_{in}\right)}_f={\displaystyle \frac{e^2}{12}}$$

(5)

the relationship between the mechanical gap width a, the hydraulic gap width e, and the roughness correction coefficient C can be derived according to Equations (4) and (5):

$${\displaystyle \frac{a^2}{e^2}}=C$$

(6)

In research on the effect of the gap width on seepage, different corrected formulas for the cubic law and for the roughness correction coefficient have been proposed. R. C. Liu [20] summarized these formulas, and the collated results are presented in

Table 1. Representative empirical relations that mathematically link the hydraulic aperture and the mechanical aperture.

Author/Year	Formula	Description
Lomize, 1951 [20]	${\mathrm{e}}^2={\displaystyle \frac{{\mathrm{a}}^2}{\mathrm{C}}}$, $\mathrm{C}=1+\mathrm{m}{\left({\displaystyle \frac{∆}{2\mathrm{a}}}\right)}^{1.5}$, $\mathrm{m}=17$	$\mathrm{a}$ is the magnitude of the discontinuity surface roughness, and m is a coefficient.
Louis, 1969 [21]	${\mathrm{e}}^2={\displaystyle \frac{{\mathrm{a}}^2}{\mathrm{C}}}$, $\mathrm{C}=1+\mathrm{m}{\left({\displaystyle \frac{∆}{2\mathrm{a}}}\right)}^{1.5}$, $\mathrm{m}=8.8$
Quadros, 1982 [20]	${\mathrm{e}}^2={\displaystyle \frac{{\mathrm{a}}^2}{\mathrm{C}}}$, $\mathrm{C}=1+\mathrm{m}{\left({\displaystyle \frac{∆}{2\mathrm{a}}}\right)}^{1.5}$, $\mathrm{m}=20.5$
Patir and Cheng, 1978 [22]	$\mathrm{e}=\mathrm{a}{\left[1-0.9\mathrm{e}\mathrm{x}\mathrm{p}\left(-0.56{\displaystyle \frac{\mathrm{a}}{{\mathsf{\sigma }}_{\mathrm{E}}}}\right)\right]}^{{\displaystyle \frac{1}{3}}}$	${\mathsf{\sigma }}_{\mathrm{E}}$ is the standard deviation of the varying aperture over the mechanical aperture.
Walsh, 1981 [23]	${\mathrm{e}}^3={\mathrm{a}}^3{\displaystyle \frac{1-\mathrm{C}}{1+\mathrm{C}}}$	C = 0.25 in the study by Zimmerman et al. (1996) [24].
Barton, 1978 [5]	$\mathrm{e}={\displaystyle \frac{{\mathrm{a}}^2}{{\mathrm{J}\mathrm{R}\mathrm{C}}^{2.5}}}\left({\mathrm{e}}^2={\displaystyle \frac{{\mathrm{a}}^2}{\mathrm{C}}}\mathrm{a}\mathrm{n}\mathrm{d} \mathrm{C}={\displaystyle \frac{{\mathrm{J}\mathrm{R}\mathrm{C}}^5}{{\mathrm{e}}^2}}\right)$
Hakami, 1995 [25]	${\mathrm{e}}^2={\displaystyle \frac{{\mathrm{a}}^2}{\mathrm{C}}}$, $\mathrm{C}=1.1~1.7$	a = 100–500 µm.
Renshaw, 1995 [26]	$\mathrm{e}=\mathrm{a}{\left[1+{\displaystyle \frac{{\mathsf{\sigma }}_{\mathrm{E}}^2}{{\mathrm{a}}^2}}\right]}^{-{\displaystyle \frac{1}{2}}}$
Zimmerman and Bodvarsson, 1996 [24]	${\mathrm{e}}^3\approx {\mathrm{a}}^3\left[1-1.5{\displaystyle \frac{{\mathsf{\sigma }}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}^2}{{\mathrm{a}}^2}}+\dots \right]\left(1-2\mathrm{c}\right)$	C is the contact ratio.
Waite et al., 1999 [27]	$\mathrm{e}=\mathrm{a}〈\mathrm{a}〉/{\mathsf{\tau }}^{1/3}$	〈a〉 is the harmonic mean of the true aperture, and τ is tortuosity.
Matsuki et al., 1999 [28]	${\mathrm{e}}^3={\mathrm{a}}^3\left[1-{\displaystyle \frac{1.13}{{1+0.191\left(2\mathrm{a}/{\mathsf{\sigma }}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}\right)}^{1.93}}}\right]$	${\mathsf{\sigma }}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}$ is the standard deviation of the mean mechanical aperture.
Olsson and Barton, 2001 [29]	$\mathrm{e}={\displaystyle \frac{{\mathrm{a}}^2}{{\mathrm{J}\mathrm{R}\mathrm{C}}^{2.5}}}\left({\mathrm{u}}_{\mathrm{s}}\le 0.75{\mathrm{u}}_{\mathrm{s}\mathrm{p}}\right)$
Xiong et al., 2011 [30]	$$\begin{gathered} \mathrm{e}={\mathrm{a}}^{1/2}{\mathrm{J}\mathrm{R}\mathrm{C}}_{\mathrm{m}\mathrm{o}\mathrm{b}}\left({\mathrm{u}}_{\mathrm{s}}\ge {\mathrm{u}}_{\mathrm{s}\mathrm{p}}\right) \\ {\mathrm{e}}^3={\mathrm{a}}^3\left[1-1.0{\displaystyle \frac{{\mathsf{\sigma }}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}}{\mathrm{a}}}\right]\left(1-{\displaystyle \frac{{\mathsf{\sigma }}_{\mathrm{a}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}}}{\mathrm{a}}}{\displaystyle \frac{\sqrt{{\mathsf{\sigma }}_{\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}}}{10}}\sqrt{\mathrm{R}\mathrm{e}}\right) \end{gathered}$$	${\mathrm{J}\mathrm{R}\mathrm{C}}_{\mathrm{m}\mathrm{o}\mathrm{b}}$ is the mobilized value of the JRC, and ${\mathrm{u}}_{\mathrm{s}}$ is the shear displacement, which does not exceed 75% of the peak shear displacement ${\mathrm{u}}_{\mathrm{s}\mathrm{p}}$. ${\sigma }_{slope}$ is the standard deviation of the local slope of the fracture surface.
Rasouli and Hosseinian, 2011 [31]	$\mathrm{e}=\mathrm{a}{\left[{\left(1-0.03{\mathrm{d}}_{\mathrm{m}\mathrm{c}}^{-0.565}\right)}^{{\mathrm{J}\mathrm{R}\mathrm{C}}_{\mathrm{a}}}\right]}^{{\displaystyle \frac{1}{3}}}$ or $\mathrm{e}=\mathrm{a}{\left(1-2.25{\displaystyle \frac{{\mathsf{\sigma }}_{\mathrm{E}}}{\mathrm{a}}}\right)}^{{\displaystyle \frac{1}{3}}}$	$d_{mc}$ is the minimum closure distance, and ${JRC}_a$ is the average of the joint roughness coefficients for the upper and lower rock fracture profiles.
Li and Jiang, 2013 [15]	$\mathrm{e}={\displaystyle \frac{\mathrm{a}}{1+{\mathrm{Z}}_2^{2.25}}}$, Re < 1 $\mathrm{e}={\displaystyle \frac{\mathrm{a}}{1+{\mathrm{Z}}_2^{2.25}+\left(0.0006+0.004{\mathrm{Z}}_2^{2.25}\right)\left(\mathrm{R}\mathrm{e}-1\right)}}$, Re ≧ 1	$Z_2$ is the root mean square of the first deviation of the profile.
Xie et al., 2015 [32]	${\mathrm{e}}^3={\mathrm{a}}^3\left[0.94-5.0{\displaystyle \frac{{\mathsf{\sigma }}_{\mathrm{s}}^2}{{\mathrm{E}}^2}}\right]$	${\sigma }_s$ is the standard deviation of the mechanical aperture during shear.
Zhang et al., 2019 [33]	${\mathrm{e}}^2={\displaystyle \frac{{\mathrm{a}}^2}{\mathrm{C}}}$, $\mathrm{C}=1+\mathrm{m}{\left({\displaystyle \frac{∆}{2\mathrm{a}}}\right)}^{0.61912}$, $\mathrm{m}=0.6138$

:

Since calculating the JRC (joint roughness coefficient) using a two-dimensional rough curve solution is not applicable in this study, the formula for calculating the JRC will not be discussed further in this paper. In comparing the formulas presented in Table 1, it is evident that Louis’ formula is older within the related research and has thus been cited more frequently. Many subsequent corrected formulas have been derived from Louis’ formula according to discussion and modification. Notably, Zhang proposed an m value of less than 1, which is uncommon among similar corrected formulas. Therefore, we have chosen to discuss the formulas proposed by Louis and Zhang in detail.

Louis’ corrected formula for the cubic law recommends a roughness correction coefficient (C).

$$\mathrm{C}=1+8.8{(∆/2\mathrm{a})}^{1.5}$$

(7)

where $∆/a$ is the relative variation in the terrain. However, Louis did not explore $∆/a$, while Zhang Youtian summarized this aspect in detail. The formula $∆={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left|h_i^{∆}-h_{i-1}^{∆}\right|$ is such that as n increases, the value of $∆$ becomes smaller, which clearly does not hold. Therefore, the statistical value ($h_i^{∆}-h_{i-1}^{∆}$) is used to calculate the variance ${\left({\sigma }_h^{∆}\right)}^2$, where greater variance represents a rougher surface [19]. Since the value ${\sigma }_h^{∆}$ is the statistical value within the range for the crack, it is more representative (hereinafter, ${ \sigma }_h^{∆}$ is referred to as ${\sigma }_s$).

Zimmerman [24], Xiong [30], and Matsuki [28] also proposed corrected formulas for the cubic law, and all three of them discussed the relationship between the hydraulic gap width (e, in mm units) and the mechanical gap width (a, in mm units) by analyzing and studying the local cubic law, with each suggesting a formula suited to their respective studies. Of these, the final relationships determined by Zimmerman and Xiong only differ in the coefficient of the (${\sigma }_s^2/a^2$) term. To observe the difference, both are calculated together. The formula proposed by Matsuki et al. is represented in Table 1 as a non-integer power formula. Ultimately, various formulas express the relationship between the mechanical gap width and the hydraulic gap width as a ratio, as shown in

Table 2. Formulas for the relationship between fractured width and roughness.

Author/Year	Fieldwork Symbols	Formula	Description
Louis, 1969 [21]	Eq-L	${\displaystyle \frac{\mathit{e}}{\mathbf{a}}}=\mathbf{C}$	$\mathrm{C}=1+\mathrm{m}{\left({\displaystyle \frac{∆}{2\mathrm{a}}}\right)}^{1.5}$, $and \mathrm{m}=8.8$.
Zimmerman R.W. et al., 1996 [24]	Eq-Zi	${\displaystyle \frac{\mathit{e}}{\mathbf{a}}}=\sqrt[\mathbf{3}]{\mathbf{1}-\mathbf{1.5} {\mathsf{\sigma }}_{\mathbf{s}}^{\mathbf{2}}/{\mathit{a}}^{\mathbf{2}}}$	${\mathsf{\sigma }}_{\mathrm{s}}$ equals 0, and the ratio of e/a approaches 1.
Matsuki K. et al., 1999 [28]	Eq-M	${\displaystyle \frac{\mathit{e}}{\mathbf{a}}}=\sqrt[\mathbf{3}]{\mathbf{1}-{\displaystyle \frac{\mathbf{1.13}}{\mathbf{1}+\mathbf{0.191}{\left(\mathbf{2}\left(\mathbf{a}/{\mathsf{\sigma }}_{\mathbf{s}\mathbf{0}}\right)\right)}^{\mathbf{1.93}}}}}$	${\mathsf{\sigma }}_{\mathrm{s}0}$ is the standard deviation of the initial aperture value.
Xiong X. et al., 2011 [30]	Eq-X	${\displaystyle \frac{\mathit{e}}{\mathbf{a}}}=\sqrt[\mathbf{3}]{\mathbf{1}-\mathbf{1.0}{\mathsf{\sigma }}_{\mathbf{s}}^{\mathbf{2}}/{\mathit{a}}^{\mathbf{2}}}$
Zhang et al., 2019 [19]	Eq-Zh	${\displaystyle \frac{\mathit{e}}{\mathbf{a}}}=\mathbf{C}$	$\mathrm{C}=1+\mathrm{m}{\left({\displaystyle \frac{∆}{2\mathrm{a}}}\right)}^{0.61912}$, and $\mathrm{m}=0.6138$.

.

Let A represent the right part of each formula in Table 2. Combining (6) with each formula yields different expressions for C, as follows:

$$\mathrm{C}={\displaystyle \frac{1}{{\mathrm{A}}^2}}$$

(8)

In incorporating the respective formula A, different expressions for C can be obtained. Each expression of C can then be used to operate the numerical model. This paper substituted the formulas derived from Table 2 (a total of 5) into Equation (8) to determine the roughness correction coefficient.

3.2. The Roughness Correction Coefficient—The Numerical Method

To calculate the roughness correction coefficient for each set of roughness data, a fracture flow model was developed using a flow field simulation. Normal temperature and pressure conditions were set in the model, the fluid material was water, the permeability model applied involved the cubic law and a steady state, and the formula for the roughness correction coefficient was consistent with Equation (8). The pore size was calculated using the corresponding imported slit width data. To simplify the numerical calculations, the intercepted data on the fractured surfaces and their gap widths were rotated by the same angle in the opposite direction according to the angle of rotation of the intercept box to ensure that the limit of the data range was parallel with the axis of the coordinates. The water flow was determined in the x direction, and the initial head was set to 0. The data for each group, including data on the fractured surfaces and the corresponding gap widths, were imported into the established model to calculate the roughness correction coefficient. Figure 4 is a schematic diagram of the model flow and the flow field.

4. Results

4.1. Results of the Scale-Effect Calculation

The roughness correction coefficient was determined for each formula and numerical method listed in Table 2 (represented as “Num”). The coefficient was calculated for three rock fracture samples consisting of granite, tuff, and griotte, respectively. The results are presented in

Table 3. Results of calculating roughness correction coefficient (S1).

Formula	Truncated Frame Size (Percent)
	10%	20%	30%	40%	50%	60%	70%	80%	90%	100%
Eq-L	1.0243	1.0186	1.0167	1.0159	1.0153	1.0147	1.0142	1.0142	1.0143	1.0145
Eq-Zh	1.7910	2.5223	3.1271	3.7963	4.5254	5.4118	6.5592	8.0026	9.7472	11.586
Eq-X	1.0824	1.0886	1.1422	1.1354	1.1310	1.1284	1.1262	1.1267	1.1247	1.1241
Eq-Zi	1.1304	1.1408	1.2352	1.2227	1.2146	1.2099	1.2060	1.2069	1.2034	1.2023
Eq-M	1.1241	1.1321	1.1961	1.1886	1.1836	1.1806	1.1781	1.1787	1.1764	1.1757
Num	1.2589	1.3928	1.8186	1.6210	1.6111	1.6000	1.5193	1.5424	1.5100	1.2376

,

Table 4. Results of calculating roughness correction coefficient (S2).

Formula	Truncated Frame Size (Percent)
	10%	20%	30%	40%	50%	60%	70%	80%	90%	100%
Eq-L	1.0158	1.0089	1.0090	1.0076	1.0065	1.0061	1.0059	1.0059	1.0060	1.0060
Eq-Zh	3.4128	4.8593	6.6603	8.4522	10.614	13.780	17.176	20.749	24.481	28.357
Eq-X	1.0366	1.0424	1.0565	1.0515	1.0489	1.0488	1.0506	1.0500	1.0489	1.0477
Eq-Zi	1.0563	1.0653	1.0879	1.0799	1.0758	1.0756	1.0784	1.0775	1.0758	1.0737
Eq-M	1.0600	1.0686	1.0889	1.0819	1.0782	1.0780	1.0805	1.0798	1.0782	1.0763
Num	1.1498	1.1579	1.2712	1.1850	1.1850	1.1373	1.1688	1.2037	1.1960	1.1958

,

Table 5. Results of calculating roughness correction coefficient (S3).

Formula	Truncated Frame Size (Percent)
	10%	20%	30%	40%	50%	60%	70%	80%	90%	100%
Eq-L	1.0611	1.0602	1.0354	1.0244	1.0213	1.0170	1.0161	1.0155	1.0148	1.0136
Eq-Zh	2.2180	4.3623	5.0022	5.3133	6.2334	6.5405	7.6451	8.9055	10.049	10.661
Eq-X	1.1434	1.1758	1.1470	1.1599	1.1581	1.1631	1.1642	1.1576	1.1446	1.1382
Eq-Zi	1.2372	1.2978	1.2439	1.2675	1.2644	1.2736	1.2758	1.2632	1.2393	1.2275
Eq-M	1.1977	1.2324	1.2014	1.2156	1.2136	1.2192	1.2203	1.2132	1.1991	1.1918
Num	1.4868	2.0036	1.4851	1.4733	1.4509	1.5515	1.5969	1.6171	1.5186	1.6022

,

Table 6. Results of calculating roughness correction coefficient (S4).

Formula	Truncated Frame Size (Percent)
	10%	20%	30%	40%	50%	60%	70%	80%	90%	100%
Eq-L	1.0816	1.0582	1.0118	1.0099	1.0089	1.0078	1.0073	1.0069	1.0067	1.0066
Eq-Zh	2.5210	2.9724	1.6342	1.7491	1.8106	1.7956	1.7666	1.7706	1.8230	1.9185
Eq-X	1.0112	1.0255	1.3840	1.5634	1.4650	1.4364	1.4604	1.4696	1.4296	1.3916
Eq-Zi	1.0169	1.0392	1.8956	2.4361	2.0307	1.9503	2.0221	2.0467	1.9156	1.8015
Eq-M	1.0198	1.0418	1.3862	1.5077	1.4563	1.4375	1.4527	1.4589	1.4351	1.4108
Num	1.0916	1.1460	1.4409	2.2475	2.0752	1.9182	2.0088	2.1018	2.0138	1.7916

and

Table 7. Results of calculating roughness correction coefficient (S5).

Formula	Truncated Frame Size (Percent)
	10%	20%	30%	40%	50%	60%	70%	80%	90%	100%
Eq-L	1.0598	1.0376	1.0360	1.0291	1.0231	1.0189	1.0141	1.0132	1.0129	1.0126
Eq-Zh	1.7806	2.3314	2.9868	3.4119	3.8425	4.2180	4.0435	4.3534	4.6486	4.8899
Eq-X	2.3916	1.3681	1.3673	1.3575	1.3238	1.3018	1.2979	1.3043	1.3203	1.3562
Eq-Zi	9.8563	1.7369	1.7371	1.7080	1.6216	1.5691	1.5587	1.5741	1.6139	1.7067
Eq-M	1.7474	1.3946	1.3936	1.3873	1.3626	1.3453	1.3425	1.3475	1.3597	1.3859
Num	5.0726	1.8825	2.4237	2.1008	2.0180	1.8159	1.6998	1.7720	1.8208	1.9584

. Figure 5 depicts the corresponding imagery.

From the results of the calculations in Table 3 and Figure 5, we determine the following:

(1)

The roughness correction coefficient obtained using the Eq-Zh calculations exhibits a significant overestimation, deviating greatly from the actual values, particularly for S2, attaining a maximum value of 30. On the other hand, the roughness correction coefficient determined using Eq-L is excessively uniform, limiting its ability to capture the variations in the roughness levels across fissure surfaces. Consequently, we will not discuss or analyze the results obtained using Eq-Zh and Eq-L further.

(2)

The calculations using the different types of formulas indicate that the roughness correction factor changes minimally when varying the fracture size for shear fractures (S1, S2, S3), whereas the roughness correction factor varies with fluctuations in the fracture size for tensile fractures (S4, S5), but the extent of the variation differs. For instance, the variation when using Eq-Zi is relatively significant, aligning with our actual understanding of the fracture roughness.

(3)

Irrespective of the type of rock, the roughness correction factors vary across varying sizes, predominantly ranging from 1.0 to 3.0, with the exception of a minimal number of outlier data points. Across these differences, S5 exhibits the highest degree of variation in the roughness correction factor compared to the others.

(4)

The results of the calculations using the Eq-X and Eq-M formulas exhibit the greatest stability across various types of rock.

(5)

The roughness correction coefficient calculated using these numerical methods varies significantly with sample size, with fewer fluctuations observed in the compressive-shear rock samples compared to the tensile rock samples. For the different calculation methods used, all of the samples except for S5 exhibit a noticeable trend of the roughness correction coefficient gradually stabilizing as the size increases, indicating that scale has an effect on roughness. In comparison, compressive-shear fractures tend to stabilize more readily in terms of their roughness, while larger sizes are required for tensile fractures to stabilize.

4.2. Results of Calculating Directionality

Because of the square frame applied to the data when using the formulas, we cannot judge directionality. Zhang Youtian [19] proposed using an elliptic diagram for permeability to derive a tensor of the roughness correction coefficient. Analyzing the elliptic diagram provides a reliable estimate of the tensor. In addition, this method has been combined with other previous approaches to determine the directionality of the permeability coefficient in studies on seepage in fracture networks. To determine whether fractured rock masses can be represented by equivalent continuous media, Long et al. [25] drew ellipses to assess the permeability coefficient. Similarly, Liu Richeng et al. [26] utilized equivalent elliptical curves for the permeability coefficient, drawn based on measurements of fractured openings and their intercept lengths, to assess whether the requisite conditions for using an equivalent continuum were satisfied. After taking osmotic measurements, Marcus (1962) and other researchers [24] drew ellipsoids to determine the direction of penetration, represented by the maximum direction of the primary axis of the spheroid. Liu Richeng [26] also noted that the difference between the long and short axes of the elliptic curve and its smoothness are critical to determining directionality.

This study was inspired by this previous research. Here, the roughness correction coefficient value was obtained using the numerical method, and elliptical graphs were analyzed for research and evaluation purposes. We surmised that being able to obtain a well-defined elliptical curve would suggest that the correction coefficient could be represented as a tensor. The results obtained by calculating S1–S5 using numerical methods were fitted using ellipses, as shown in Figure 6.

The long and short axes of the ellipse indicate the maximum and minimum values of the roughness correction coefficient, respectively, and directionality is judged according to the size difference between the angle between the two axes and the horizontal axis [25]. The statistics measured show the angles between the short axis and the horizontal axis in Figure 6, and the results are shown in

Table 8. Oval rotation angles in Figure 6.

Graphic Number	S1	S2	S3	S4	S5
Oval rotation angle (°)	53.60	−31.93	−3.21	43.76	14.68
Long axis (mm)	3.11	1.27	1.84	1.908	1.32
Short axis (mm)	1.11	1.08	1.26	1.700	0.98

. In the table, we can see the following:

(1)

The angle between the long axis, the short axis, and the horizontal line varies in each figure, resulting in an anisotropic roughness correction coefficient.

(2)

The roughness correction coefficient varies with direction for different types of rocks.

(3)

The variation in the roughness correction coefficient is indicative of the variation in the roughness of the rock samples. Additionally, the roughness will differ in different directions within the same sample, leading to variations in the flow velocity of fluid in different directions.

The rotation angle can be expressed as a matrix for each angle considering the angle of rotation in the horizontal direction. It should be noted that θ represents the rotation angle, and its corresponding matrix R is [9]:

$$\mathrm{R}=\left(\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }& \mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }\\ -\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\theta }& \mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }\end{array}\right).$$

The rotation matrix from the ellipse method for permeability can be applied to solving the roughness correction coefficient [10,12]. In the relation formula,

C_an = R C_digR^T

C_an represents the roughness correction coefficient tensor matrix for the off-diagonal array, while the array of diagonal roughness correction coefficients is referred to as C_dig. Additionally, R denotes the rotation matrix in this equation. Since the ellipses studied are mainly two-dimensional with polar coordinates, we express them in the two-dimensional form $C_{\mathrm{d}\mathrm{i}\mathrm{g}}$=$\left(\begin{array}{cc}{\mathrm{c}}_{\mathrm{a}}& 0\\ 0& {\mathrm{c}}_{\mathrm{b}}\end{array}\right)$. Here, the longer and shorter axes of the ellipse are represented by a and b, respectively. Furthermore, C_an can be used, with $C_{an}$=$\left(\begin{array}{cc}{\mathrm{c}}_{11}& {\mathrm{c}}_{12}\\ {\mathrm{c}}_{21}& {\mathrm{c}}_{22}\end{array}\right)$, and we use C_ij as a simple notation for the coefficient tensor matrix C_an, denoted as C_an = C_ij. Analyzing Figure 5 and Figure 6 indicates that increased ellipse sizes improve the alignment. To apply the formula, the angles of the fitted ellipses in Table 8 should be put in order based on their maximum size and substituted into the formula above.

We observe the following data in

Table 9. Results of calculating roughness correction coefficient tensors.

Graphic Number	S1		S2		S3		S4		S5
Can	Ci1	Ci2	Ci1	Ci2	Ci1	Ci2	Ci1	Ci2	Ci1	Ci2
C1j	1.8126	−0.9559	1.2188	0.0881	1.8145	0.1293	1.8080	−0.1038	0.1504	1.2293
C2j	−0.9559	2.4045	0.0881	1.1324	0.1293	1.2937	−0.1038	1.7990	0.1504	1.2293

:

(1)

The absolute value of C_ij falls between 0 and 3, which mostly correlates with the calculated results from Table 8 (1.0~3.0); that is, the results calculated for the roughness coefficient matrix are similar to the original results.

(2)

The roughness correction coefficient matrix and the roughness coefficient matrix for the granite sample both exhibit symmetric patterns on the secondary diagonal, which suggests stable surface roughness in all directions.

In summary, we can conclude that the roughness correction coefficient is both anisotropic and symmetrical and can be expressed in tensor format. As a parameter reflecting fractured roughness, it is apparent that roughness is also subject to orientation and symmetry. Additionally, it has been observed that the surface roughness of rocks depends on both the dimensions of the sample and the scale of observation [27]. However, these results require further interpretation.

5. Discussion

According to fractal theory, the fractal dimension is the main indicator of the complexity and irregularity of rock pores. It is also an important parameter in that it represents the degree of self-similarity of the pore distribution, the proportion of pits and protrusions on a rough surface, and the rough, uneven characteristics of rock joint surfaces fully [28]. The fractal dimension of a fractured surface, in reflecting the roughness and gap width, is indicative of the permeability of a single slit [9,29]. In his study, Ju Yang [9] confirmed that roughness and the fractal dimension influence fractured permeability and proposed an equivalent fractal permeability coefficient for fractures, as well as a formula for calculating it, analyzing the fractal dimension of fractures and incorporating the pore structure and the direction and rate of seepage.

We calculated the fractal dimensions for our data on the fractured surfaces and their gap widths, and our results are illustrated in

Table 10. Results of calculating the fractal dimension based on the roughness correction coefficient.

Sample	Cut Frame Size (Percentage)/%
	10	20	30	40	50	60	70	80	90	100
S1-AP	2.1552	2.2319	2.3103	2.3594	2.3945	2.4137	2.4375	2.4416	2.4592	2.4759
S1-Sur	2.0949	2.1315	2.1057	2.0985	2.1143	2.1114	2.1003	2.0985	2.0908	2.2086
S2-AP	2.1919	2.2632	2.287	2.3206	2.3223	2.3458	2.3647	2.3818	2.4056	2.4222
S2-Sur	2.0864	2.0748	2.0926	2.1026	2.0762	2.0632	2.0602	2.06	2.0532	2.0515
S3-AP	2.2115	2.3295	2.3651	2.3832	2.3988	2.4132	2.4175	2.4239	2.4327	2.4157
S3-Sur	2.1055	2.1005	2.1138	2.0697	2.0762	2.0664	2.0683	2.0627	2.0681	2.0862
S4-AP	2.2344	2.3521	2.23	2.2741	2.2747	2.2343	2.2412	2.2043	2.2032	2.201
S4-Sur	2.1524	2.1358	2.1112	2.1299	2.1693	2.1643	2.1714	2.1899	2.184	2.2084
S5-AP	2.2739	2.322	2.3264	2.3644	2.3517	2.3373	2.3075	2.3200	2.3256	2.3333
S5-Sur	2.1074	2.0665	2.0978	2.1478	2.1307	2.1429	2.1524	2.1477	2.163	2.1537

Note: S1-AP represents the fracture width of sample S1, S1-Sur represents the surface of sample S1, and the same principles apply to the notation for the other samples and their fracture widths and surfaces.

. When these are combined with Figure 7, the following is evident:

(1)

The fractal dimension of each material is different. Additionally, for a given material, the fractal dimension of the fractured width is higher than that of its fractured surface. These observations hold true for the roughness of granite, as investigated in this study.

(2)

The fractal dimension of a given face changes with the scale of the observation. However, the fractal dimension curve stabilizes as the scale of observation increases, as previously found by Fardin [30]. Some of the curves, such as S1-Sur, were not flattened, indicating that the size of our samples may have been insufficient. Indeed, in their investigation of the heterogeneity of roughness, Shi-Gui, D. et al. [31] emphasized that their results were significantly influenced by sample size.

(3)

The fractal dimension of natural rough joints tends to be low and was found to be a non-integer varying between 2 and 2.5 for each study area here. This variation is suggestive of the anisotropy of the surfaces studied. In conclusion, the roughness of the fractured surfaces is subject to the effects of scale and anisotropy.

6. Conclusions

In conclusion, by combining the detailed results of our calculations and our analyses of data from four samples, we derive the following findings:

(1) The roughness correction coefficients of the samples consisting of shear-type rock exhibit slight fluctuations as the scale varies, while the same values for the tensile rock samples experience significant fluctuations. Nonetheless, the roughness correction coefficients of all the rock samples are affected by the scale of the fracture samples, thereby indicating that scale has an effect on roughness.

(2) Different formulas for calculating the roughness correction coefficient produce distinct results and exhibit varying sensitivity to size. However, in general, these coefficients undergo some degree of change as the scale of the fracture sampling changes, providing further evidence that scale affects the roughness of fracture surfaces. The roughness of rock fracture surfaces also demonstrates anisotropy, which signifies its directional characteristics.

(3) The roughness correction coefficient in a polar coordinate system can be fit to an elliptical shape, suggesting that the anisotropy of the coefficient can be represented using a tensor form.

(4) One reason why the effect of scale is more pronounced in the tensile rock samples compared to the shear-type rock samples is the substantial disparity in the fissure range and the distance between the local top and bottom surfaces, with less variation in the fissures in the shear-type rock samples.