#### Abstract

The joint roughness coefficient (JRC) is an important factor affecting the shear properties of rock joints, and its accurate estimation is a challenging task in rock engineering. Existing JRC evaluation approaches such as the empirical comparison method and the statistical parameter method have some unresolved defects. In this study, a new method is proposed for JRC estimation to overcome the deficiencies of existing approaches based on back calculation of shear strength. First, the 10 standard roughness joints are established in numerical rock samples generated by the bonded particle method (BPM). Secondly, the microscopic parameters of the intact rock and joints are calibrated, and a series of direct shear tests of joint samples are carried out under different normal stresses. Finally, the empirical relationships between shear strength and JRC are proposed under high correlation conditions. The results show that the modified smooth joint model (MSJM) is proved to better simulate the mechanical properties of rough joints than the smooth joint model (SJM). When the shear strength of target joint is substituted in the corresponding relationship, the JRC of joint along the shear direction can be conveniently obtained. In addition, the JRC values of 10 standard roughness joint profiles under shear direction of from right to left (FRTL) are obtained. By estimating the JRC of 9 target joints in the literature, it can be seen that the new method proposed in this paper can well reflect the directionality of roughness and it is convenient to apply.

#### 1. Introduction

The mechanical behaviour of rock mass is much more complex than intact rock because of its intricate joints and structural planes. The safety and stability of rock mass are mainly controlled by shear deformation along the joints, which is particularly prominent in slope and tunnel engineering [1–3]. For unfilled rock joint, the irregular morphology of joint surface is one important factor to control the shear strength of joints; therefore, how to describe it quickly and accurately is of crucial importance [4–7]. Barton [6] proposed a joint shear strength formula through the direct shear test of 136 natural joint specimens. The joint roughness coefficient (JRC) was put forward to quantitatively describe the irregular morphology of joint at the first time.

Based on the results of direct shear tests, Barton and Choubey [7] proposed 10 standard roughness profiles. The JRC value of the standard profile is obtained by back calculation from equation (1) with a range of 0 to 20, as shown in Table 1. By comparing the surface morphology of natural joints with the standard profiles, the most similar profile and corresponding JRC value are selected (for example, JRC = 8–10). This process is called the empirical comparison method. The International Society for Rock Mechanics [8] recommend these 10 profiles as standard method for evaluating the joint roughness coefficient, because the empirical comparison method gives the way to determine JRC value for the first time.

The empirical comparison method mentioned above is very convenient and fast. However, the process is greatly affected by experience of engineer and the JRC value fluctuates greatly [9–12]. Therefore, in order to obtain the roughness coefficient of joints more accurately, statistical parameters such as the straight edge method [13], the fractal dimension (D) [14, 15], root mean square roughness index (RMS) and mean square value roughness index (MSV) [16], root mean square of the first deviation of profiles (*Z*_{2}) [16–21], structure function (SF) [16–20, 22], standard deviation of the angle (SDi) [17], and the roughness profile index (R_{P} − 1) [18–22] have been presented. These parameters can be calculated from geometric coordinates of joint profile as shown in Figure 1, where *x*_{i} and *y*_{i} are the abscissa and ordinate of joint profile, Δx is the distance between the abscissa of adjacent coordinate point and becomes a fixed value at the same sampling interval, and *L* is the horizontal length of joint profile [23].

At present, many researchers believe that the rock joint surface has anisotropy. The joint usually shows different strength and deformation characteristics when shearing from different directions [24–26]. In the JRC-JCS model proposed by Barton, parameters other than JRC cannot reflect the directionality. Therefore, it can be confirmed that the directivity of the strength and the deformation is caused by the directivity of JRC. The surface morphology of natural joints is random, and it is impossible to have the same feature in different directions. The default shear direction of 10 standard profiles proposed by Barton is from left to right (FLTR), which means that the JRC values in parentheses (Table 1) represent the roughness under the condition of FLTR. JRC values of 10 standard profiles under the condition of from right to left (FRTL) is unknown. When the above parameters are used to predict the JRC, the estimated shear strengths of rough joints are the same at the shear condition of FLTR and FRTL. Obviously, this situation does not meet the test rules.

In order to solve this problem, new parameters such as maximum inclination () [19, 21], modified root mean square () [27, 28], and average negative angle (*β*_{100%}) [20] were proposed to reflect the directionality of the joint morphology. When defining these parameters, only one side of the joint morphology is considered and the other side is completely ignored. Although these parameters were verified to be highly correlated with JRC, there are still some problems in the above assumptions. When considering the joint 1 and the joint 2 in Figure 2, their JRC under the condition of FLTR is the same, which is obviously unreasonable.

In addition, it has been found that most of the existing parameters cannot distinguish between the joint 1 and joint 2 in Figure 2, because the inclination and height of the joint profile cannot be described by one parameter at the same time [27]. Although a comprehensive parameter that takes into account inclination and height is proposed by Zhang et al. [27] to compensate for this deficiency, its correlation with the JRC is not sufficiently good. The contribution of joint morphology to strength and deformation characteristics during shearing is unknown at the moment, and the application of the statistical parameter method is greatly limited.

In the past ten years, numerical simulation is becoming more and more widely used in geotechnical engineering [29–35]. With the continuous innovation of computational mechanics, the discrete element method (DEM) in particulate and blocky systems has a rapid development [30–39]. The potential advantage of this approach is that the response of jointed rock mass can be derived based on simple particle contact logic in jointed rock [40].

In order to solve the shortcomings of the current roughness estimation methods, this paper presents a back calculation idea for JRC estimation based on numerical simulation methods. Numerical rock samples containing rough joints were established using Particle Flow Software (PFC^{2D}). The joints are modelled using smooth joint model (SJM) and the modified smooth joint model (MSJM) [33–35, 39, 40]. A series of direct shear tests of rock joints are carried out under different normal stresses. Based on the test results, the fitting relationship between JRC and shear strength was established with high correlation. The JRC of target joint profile can be obtained by back calculation of shear strength. This process mentioned above avoids the difficulty of describing the morphology of joint surface and the directionality of roughness can be well considered. This method is very simple and convenient and can be well applied in practical engineering.

#### 2. Proposed New Method for Estimating JRC of Joint Profile

Based on the 10 standard roughness joint profiles, a new JRC determination method is proposed. In the first step, the geometric coordinates of the 10 standard roughness joint profiles are acquired using image processing technology. In the second step, these joints are set into a numerical model of rock sample generated by the bonded particle method (BPM) in Particle Flow Code (PFC^{2D}) and simulated by the modified smooth joint model (MSJM) [30–40]. Secondly, a servo mechanism is added to make the normal stress on the upper and lower surfaces of joint sample maintaining constant. Next, a series of direct shear tests of joint samples are carried out under different normal stresses. Finally, the fitting relationships between shear strength and JRC are established under different normal stresses. Under the same conditions of numerical test, the direct shear test of target joint can be carried out in a certain shear direction. When the shear strength of target joint is substituted in the corresponding relationship, a JRC value is obtained for each normal stress condition and the mean of JRC values obtained by back calculation is recommended to characterize the roughness of target joint profile along the shear direction. The procedure of this method is shown in Figure 3.

In this method, the direction of the joint roughness is consistent with the direct shear test and can be easily considered. Once the fitting relationships between shear strength and JRC are established, the written command text and model can be reused. This provides great convenience for the application of the method proposed in this paper.

#### 3. Preparation for Research

##### 3.1. Digitization of Ten Standard Roughness Profiles

The original data of the 10 standard roughness profiles proposed by Barton are measured by a profile comb (Figure 4). The sampling interval of the obtained data is 0.5 mm while the spacing between the comb teeth is 1 mm. It can be seen that the measured data have a slight error from the actual morphology of joint surface. Digitization of joint profile is equivalent to reacquisition of the original data. It is not reliable to obtain the discrete data of the standard roughness profiles with sampling interval less than 0.5 mm because this sampling interval has exceeded the accuracy of the original points. Therefore, the minimum sampling interval of joint profile is set to 0.5 mm.

We firstly saved these profiles shown in Table 1 as images in jpg format by screenshot software and imported these images into GETDATA software. After setting the axes, the geometric coordinates of the profile could be acquired by the autocapture mode. Finally, the point capture mode was used to fill the vacant part and modified the singularity point. This process was similar to Jang et al. [19] as shown in Figure 5. The horizontal lengths of 10 standard roughness joint profiles were close to 100 mm and the average horizontal length was about 99 mm. Among them, the seventh profile was the shortest (96 mm) and the eighth profile was the longest (101 mm) [7, 19]. For determining the geometric coordinates, the horizontal length of the 10 joint profiles was assumed to be 100 mm and the coordinate of starting point was set to (0, 0). The horizontal length of seventh joint profile had a maximum variation and an increase by 4.167%, according to the processing idea of Zheng and Qi [20]. The JRC values of 10 joint profiles were assumed to be unchanged. Figure 6 shows the modified 10 standard roughness profiles.

##### 3.2. Verification of Profile Data

In order to verify the accuracy of joint morphology acquired previously, statistical parameters *Z*_{2} and R_{P}-1 which are widely used in engineering are recommend as representative for verification. Their definitions are as follows [16–22]:where *N* is the number of discrete points obtained from joint profile and it is 200 at a 0.5 mm sampling interval. The values of *Z*_{2} and *R*_{P} − 1 were calculated from the standard roughness profiles as shown in Figure 7. Except for the singularity obtained from the fourth profile, the values of *Z*_{2} and *R*_{P} − 1 tended to increase with the increasing of JRC values. The fitting relationship of statistical parameters and JRC were established by using equation (6) proposed by Jang et al. [19]:where *P* represents the statistical parameters and *a*, *b*, and *c* represent the regression coefficients. Under the condition of 0.5 mm sampling interval, the correlation coefficients of *Z*_{2} and *R*_{P} − 1 with JRC calculated in our paper were larger than that of Zheng and Qi [20] and Jang et al. [19] using the same function (Table 2). This was further confirmed that data of the Barton standard profiles digitized in our paper was highly accurate.

**(a)**

**(b)**

#### 4. Realization of Direct Shear Tests of Joint Sample

##### 4.1. Simulation of Joint Sample

The numerical models in this study were established by the Particle Flow Code software (PFC^{2D}). The mechanical behaviour of materials can be simulated by displacement and force interaction in particles. The basic relationship between forces and motion of particles is Newton’s second law of motion [40]. The joints were simulated using the smooth joint model and the modified smooth joint model, respectively.

###### 4.1.1. Simulation of Joint with Smooth Joint Model

The bond removal method is the preliminary method of simulating joints. But, this method has a congenital deficiency that the joint roughness is easily affected by the distribution of particles. In order to overcome the shortcomings of the bond removal method, the smooth joint model (SJM) was introduced into PFC and has been widely used in numerical studies of jointed rock masses [33–35, 39, 40]. The smooth joint model is placed at contacts between the particles whose centre is lying on the opposite sides of the joint surface. At these contacts, the existing bonds have been removed, and this process is the same as the idea of the bond removal method. In addition, the smooth joint model allows particles to overlap or pass through each other rather than to spin around each other. This function has been verified to better simulate the mechanical properties of rough joints [33–35].

For building the numerical joint model, a box made up of 6 frictionless walls was generated firstly, which could provide the boundary for particle generation or loading board to applying velocity boundary conditions. The model size was 100 mm in width and 60 mm in height. Randomly placed particles were generated in the box. The particle diameters satisfied uniform size distribution. The minimum radius was 0.15 mm, and the maximum radius was 0.24 mm. Each numerical sample contained 41 883 particles at a porosity of 0.15. In the beginning, there was a large amount of overlap between the particles. After giving an initial stiffness, the particles generated previously were allowed to rearrange and reached the static equilibrium under zero friction condition. Before installing bond, particles with less than one contact were eliminated to produce a series of highly interlocked particles [33, 34]. Then, the bonded particle model (BPM) was installed at a contact point between all particles to simulate intact rock [30–40]. Discrete fractures network (DFN) was applied to input local geometric information including midpoint, dip angle, and length of joints [33–35, 39, 40]. Since the sampling interval of the geometric coordinates was 0.5 mm, the 10 standard roughness joint profiles were sequentially divided into 200 small segments and set into the numerical model. Then, the SJM was installed in DFN to simulate joint [33–35, 39, 40]. The joint sample established by the smooth joint model is shown in Figure 8.

###### 4.1.2. Simulation of Joint with Modified Smooth Joint Model

Bahaaddini et al. [34] found that when the shear displacement exceeded the minimum diameter of the particles, the smooth joint model might lead to an unreliable increase in the shear strength and normal displacement of joint sample. This was due to the interlocking phenomenon of some particles near the joint. Therefore, the modified smooth joint model (MSJM) has been proposed, and there are some changes to the modelling steps of the joint sample [34].

In this method, two contact boxes were generated first, and each box has four frictionless walls. The bottom wall of the upper box and the upper wall of the box below were overlapping. In order to build a rough joint model, they were generated from the geometric coordinates of the rough joints. Since only one side of the wall is active, the input of the joint node followed the right-hand rule. The particle size and porosity were the same as Section 4.1.1. Then, the upper and lower parts of the joint sample were generated separately rather than once as Section 4.1.1. The particles were then rearranged and brought to a static equilibrium state. After eliminating the floating particles, the BPM was installed on the upper and lower parts of the particles, respectively, and the two irregular walls in the middle of the sample were deleted. Next, a suitable normal pressure was applied to the sample, and new contacts were generated between the particles of the upper and lower portions of the sample. Finally, the DFN was applied to determine the location of rough joint, and the smooth joint model is applied at DFN [33–35, 40]. The joint sample established by the modified smooth joint model is shown in Figure 9.

##### 4.2. Calibration of Microscopic Parameters

In PFC, the macroscopic behaviour of the material that we are expecting to simulate is derived by the interaction of microscopic components such as particle and bond. However, the microscopic parameters of these constituents are usually unknown. Therefore, these parameters must be determined by the calibration process by trial and error.

###### 4.2.1. Calibration of Microscopic Parameters of BPM

Based on the procedure explained in Section 4.1, numerical specimens with length of 50 mm and height of 100 mm were generated for calibration of microscopic parameters of BPM applied in intact rock [30–40]. The particle size and porosity are the same as Section 4.1.1. The specimen consisted of 34901 particles. The uniaxial compression test results showed the uniaxial compressive strength *σ*_{c} = 128.91 MPa and the deformation modulus *E* = 33.05 GPa for the intact rock. The macroscopic mechanical properties and failure mode of the numerical model reproduced the behaviour of the physical experiments (Figure 10). Table 3 summarized the calibrated microscopic parameters of BPM.

###### 4.2.2. Calibration of Microscopic Parameters of SJM

Microscopic parameters were calibrated by direct shear tests on planar joints. Numerical rock samples with dimensions of 100 mm × 60 mm were generated, and a planar joint was inserted in the middle of the samples by the SJM. In the direct shear tests, the walls below the joint position were kept stationary and a constant normal stress was applied to the top wall by using a servo mechanism. The sliding behaviour of the joint could be observed for different normal stresses and the peak shear strength obeyed Coulomb’s sliding model. The results showed the shear stiffness was 3.12 GPa/m and the friction angle was 32.3° for the joint. The calibrated microscopic parameters of SJM are presented in Table 4.

##### 4.3. Servo Mode

The direct shear test under constant stress conditions needs to maintain the normal stress applied to the upper and lower surfaces of sample to a constant state. However, it is not possible to apply loads directly to the wall in the PFC software. The constant load state can be maintained by changing the normal moving speed of the walls at each calculation step. This process is controlled by a servo function. The normal velocity of the upper and lower walls of sample can be set as follows:where *G* is the servo parameter, *σ*^{measure} is the normal stress actually applied to the upper and lower surfaces of the joint specimen, *σ*^{require} is the target normal stress set by us, and Δ*σ* is the stress difference between *σ*^{measure} and *σ*^{require}. The maximum value of the stress difference is defined as

In equation (2), is the average stiffness of the particles in contact with the upper and lower walls, *N*_{c} is the number of particles, and *A* is the area of the wall. In order to keep the difference between *σ*^{measure} and *σ*^{require} in a small range, it is assumed that there is a release factor *α* as shown in equation (6):

The following relationship can be obtained by substituting equations (4) and (5) in equation (6):

Thus, the servo coefficient *G* can be obtained as follows:

When the software calculates a time step, the servo parameter *G* of the next time step is calculated by the servo function. The normal velocity of the upper and lower walls is updated by (4), and the servo parameter *G* required for the next time step is also calculated simultaneously.

In order to verify the rationality of this method, the direct shear tests of planar joint under normal stress conditions of 5 MPa, 6 MPa, and 7 MPa were carried out. The normal stresses applied to the upper and lower walls were tracked as shown in Figure 11. It could be seen that the difference between *σ*^{measure} and *σ*^{require} is very small, so it could be considered that the servo process was reliable.

##### 4.4. Shearing Process

The scheme of the direct shear test is shown in Figure 12. The upper shear box consisted of 3#, 4#, and 5# walls, and the lower shear box consisted of 1#, 2#, and 6# walls. During the direct shear test, a constant horizontal velocity of 4 mm/s was applied to the upper shear box and the lower shear box was fixed. Shear stress and shear strain were calculated by recording the force and displacement of the corresponding wall. The calculations in PFC were carried out by a time-stepping algorithm. The time between each calculation step had a small value (Δ*t* ≈ 2.904 × 10^{−7} s), and therefore, the upper and lower walls were moved at the rate of 11.616 × 10^{−7} mm per time step. This movement rate was slow enough to ensure that the test process was in quasistatic equilibrium.

The direct shear tests of 10 standard roughness joint profiles were carried out, and the shearing direction was FLTR. The joints were simulated by the SJM and the MSJM, respectively [33–35, 39, 40]. The normal stresses applied to the joint specimens were 5 MPa, 6 MPa, 7 MPa, and 8 MPa respectively.

#### 5. Test Results

##### 5.1. Numerical Joint Samples Generate by the SJM

The shear direction of the direct shear test carried out previously was FLTR, which was consistent with Barton’s research work. For the convenience of describing shear direction, JRC^{LR} and JRC^{RL} were defined as joint roughness coefficient of joint under the shear direction of FLTR and FRTL, respectively. With the same idea, *τ*^{LR} and *τ*^{RL} were defined as shear strength of joint under the shear direction of FLTR and FRTL, respectively. When the shear displacement reached 5 mm, the shear stress of all joint generated by the SJM showed a peak point. The *τ*^{LR} values of 10 standard roughness joints under normal stresses of 5–8 MPa are summarized in Table 5.

It could be seen that the *τ*^{LR} value of the third, seventh, and ninth joints are less than those of the preceding joint. As the joint number increased, the *τ*^{LR} values appeared to increase as a whole. With the increase of normal stress, the *τ*^{LR} value showed an increasing trend. This was consistent with the laws of laboratory tests, and the correctness of servo mechanism in this paper was well verified. Joint surfaces usually have different friction properties, so even if the joint profile is a straight line, the joint still has a certain shear strength. Therefore, a function identical to equation (3) was used to establish the fitting relationship between JRC^{LR} and *τ*^{LR} as shown in equation (9):where *a*, *b* and *c* are regression parameters.

The JRC^{LR} values of each joint had been obtained by back calculation using equation (1) in Barton’s paper (listed in Table 1). The fitting relationships between JRC^{LR} and *τ*^{LR} under different normal stresses were established by equation (9), as shown in Figure 13. The correlation coefficient *R*^{2} of fitting relationships had a maximum value of 0.9418 and a minimum value of 0.9168. It could be seen that the mechanical properties of joints were not well simulated by the SJM.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.2. Numerical Joint Samples Generated by the MSJM

The *τ*^{LR} values of 10 standard roughness joints simulated by the MSJM are summarized in Table 6. It was worth noting that the *τ*^{LR} values of joints under different normal stress conditions increased with the increase of joint number, and there was no singularity to appear. As the normal stress applied to the specimen increased, the *τ*^{LR} values increased linearly.

When the normal stress is 5 MPa, 6 MPa, 7 MPa, and 8 MPa, the fitting relationship between JRC^{LR} and *τ*^{LR} was established, respectively, by equation (9) using the ideas in Section 5.1. It could be seen from the results in Figure 14, the correlation coefficient *R*^{2} of fitting relationships under different normal stresses was greater than 0.97. When the normal stress was 8 MPa, the correlation coefficient had a maximum value of 0.9940. When the normal stress was 5 MPa, the correlation coefficient had a minimum value of 0.9799. The mechanical properties of joints can be well represented by MSJM. The MSJM proposed by Bahaaddini et al. [34] was more capable of simulating rough joints than the SJM.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.3. Establishment of the Fitting Relationships

According to the results of Section 5.2 and 5.3, the joint sample was recommended to be simulated using the MSJM. When the shear direction is FRTL, JRC^{RL}, and *τ*^{RL} will also necessarily satisfy the fitting relationships between JRC^{LR} and *τ*^{LR}. Therefore, the fitting relationships showed in Figure 12 were transformed and integrated as equation (10) to equation (13). These equations could be used as a basis for roughness estimation and the directionality of JRC and shear strength must be consistent. It is worth noting that the simulation conditions of the target joint and the 10 standard roughness joints must be consistent; otherwise, it will lead to unknown errors.

#### 6. Verification of Proposed Method

##### 6.1. JRC^{RL} Values of 10 Standard Joint Profiles

Based on the numerical joint samples established in Section 5.2, a series of direct shear tests were carried out under the shear direction of FRTL. The *τ*^{RL} values under different normal stresses are summarized in Table 7.

In order to further verify that the directionality of joint roughness, the *τ*^{RL} values in Table 7 was used to establish fitting relationship with JRC^{LR} values. It could be seen from Figure 15 that the correlation between JRC^{LR} and *τ*^{RL} was not very good, and all correlation coefficients *R*^{2} were lower than 0.9. This was due to the fact that they are derived from different shear directions and should not be mixed.

**(a)**

**(b)**

**(c)**

**(d)**

In the physical test, the JRC values of natural joint in different directions cannot be back-calculated at the same time, because the joint surface usually has damage after the shearing process. Therefore, the JRC^{RL} values of 10 standard joint profiles are unknown. To solve this problem, the *τ*^{RL} values in Table 7 were substituted into equation (10) to equation (13) according to the method described in Section 5.3. Then, the JRC^{RL} values of 10 standard joint profiles were obtained, as showed in Table 8.

It can be found that for the first six profiles, the JRC^{RL} values did not always increased. When the shear direction was FRTL, the order of the joint profiles was rearranged according to the value of JRC^{RL} as shown in the rightmost column of Table 8. The JRC^{RL} values had a little large difference with the JRC^{LR} values. In most cases, the front was smaller than the latter except for the first, the second, and the fourth joint profile. The difference between JRC^{RL} and JRC^{LR} verified the important impact of shear direction, and they were usually affected by the irregular morphology of the joint surface.

##### 6.2. JRC^{LR} and JRC^{RL} Values of Target Joints

The 9 joint profiles with different morphologies in the literature [41] were set as the target joint and digitized by the image processing technique mentioned in Section 3.1. Their direct shear tests were carried out when the shear directions were FLTR and FRTL, respectively. The test steps and conditions for the target joints were identical to the 10 standard roughness joints. In the command text required for software calculation, only the geometric coordinates of the joint need to be changed, which saves a lot of time. The shear strengths of target joints were traced under normal stresses of 5 MPa, 6 MPa, 7 MPa, and 8 MPa and then substituted in equation (10) to equation (13). The estimated JRC^{LR} and JRC^{RL} are summarized in Table 9.

The results showed that the values of JRC^{LR} and JRC^{RL} obtained from 9 target joints were obviously different. In the third joint profile, the difference between JRC^{LR} and JRC^{RL} was as high as 12.98. In the seventh joint profile, the difference between JRC^{LR} and JRC^{RL} was only 0.1. Compared with the JRC range given in the literature, the proposed method in this paper can consider the directionality of roughness well and further reasonably predict the shear strength of joints.

The JRC estimation method proposed in this paper is based on the numerical shear test of joint profile. The values of shear strength used in equation (10) to equation (13) are only the simulation results of the numerical tests, and they are different from the meaning of actual shear strength of joint profiles in the actual engineering. When the simulation conditions of target joint profiles are the same as the 10 standard roughness joint profiles, we can get the JRC of target joint profile by equation (10) to equation (13). Subsequently, if we know the true values of other parameters in the JRC-JCS model proposed by Barton, the true shear strength of target joint profile can be easily predicted [6, 7].

#### 7. Conclusion

This paper proposes a new method for estimating the joint roughness coefficient (JRC) of joint profile. In this method, the description of the joint morphology is avoided and the shear direction can be easily considered. First, 10 standard roughness joint profiles proposed by Barton were digitized using data processing software. The calculation results of existing parameters such as *Z*_{2} and *R*_{P} − 1 were similar to previous studies, which indicated the data acquisition process was reliable.

The bonded particle model (BPM) embedded into the Particle Flow Software PFC was applied to establish the numerical model of joint sample. Then, 10 standard roughness joints were set in numerical samples through DFN and simulated by the SJM and the MSJM, respectively. After calibrating the microscopic parameters of the intact rock and joints, a series of direct shear tests of joint samples were carried out to determine shear strength. The test results indicated that the MSJM can better simulate the mechanical properties of joint than the SJM.

The fitting relationships between shear strength(*τ*) and JRC were established under different normal stress conditions. The shear strength of the target joint can be obtained by a direct shear test with same modelling conditions as the 10 standard roughness joints. When the shear strengths of the target joint under different normal stresses are substituted into the corresponding fitting relationships, the JRC value of target joint along shear direction will be obtained after an averaging process. The direct shear tests of 10 standard roughness joints were carried out at the shear direction of FRTL and the shear strengths (*τ*^{RL}) and joint roughness coefficients (JRC^{RL}) in this direction were obtained. The JRC^{RL} values were significantly different with the JRC^{LR}, which was a good proof of the importance of shear direction.

In order to verify the accuracy of proposed method, the joint roughness coefficients of 9 joint profiles in the literature were obtained. The values of the JRC^{LR} and the JRC^{RL} were significantly different. It could be well confirmed that the method proposed in this paper can not only accurately estimate the roughness of the joint, but also can easily and conveniently consider the influence of the shear direction.

In the following work, we will apply this method to the roughness evaluation of three-dimensional joint surface, and the anisotropy of joint roughness will be focused.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was funded by the National Natural Science Foundation of China (11572246, 51179153), Scientific Research Program of Shaanxi Provincial Education Department (17JS091), and Natural Science Basic Research Plan in Shaanxi Province of China (2019JQ395).