Using PFC2D to simulate the shear behaviour of joints in hard crystalline rock

Significant uncertainties remain regarding the assessment of the peak shear strength of rock joints. In recent years, Particle Flow Code (PFC) has been used to simulate shear tests of rock joints. Although previous studies showed PFC’s capability to simulate rock joint shear behaviour, it is uncertain how different parameters in PFC should be combined to realistically capture roughness and strength of asperities in contact of actual rock joints. Under low normal stresses, the shear behaviour of well-mated hard crystalline joints is governed by the interaction between asperities of some tenths of a millimetre. This paper investigates the capability of PFC2D to realistically simulate the peak shear strength of hard crystalline rock joints under different constant normal stress magnitudes. The simulated two-dimensional profiles were selected from the digitised joint surface obtained with optical scanning measurements. To realistically capture surface roughness and asperity strength in PFC2D, different values of joint segment length, particle resolution per segment, and bond strength between particles were studied and calibrated while taking into account the laboratory observations. The results of the numerical simulations in the PFC2D environment show that the simulated peak shear strength using the profile containing the steepest asperity is in good agreement with that measured in the laboratory. The joint profile needs to be represented by both a magnitude of segment length that captures the grain size, and at least two particles per segment. The bond strength calibration needs to account for both asperity size and the number of particles in contact during shearing.


Introduction
In rock mechanics, it is necessary to be able to assess the shear strength of rock joints in order to solve different rock mechanical problems; this includes, for example, slope stability, block and arching stability in tunnels and the sliding stability of dams. In recent decades, various attempts have been made to develop empirical and analytical criteria to increase the understanding of the peak shear strength of rock joints (e.g., Patton 1966;Ladanyi and Archambault 1969;Barton and Choubey 1977;Grasselli and Egger 2003;Johansson and Stille 2014;Yang et al. 2016;Casagrande et al. 2018;Ban et al. 2020;Ríos-Bayona et al. 2021, 2022. However, despite the proposed criteria, there is still uncertainty in the estimation of the shear strength of rock joints. In recent years, the emergence of advanced software has made it possible to use numerical modelling as a tool for studying the mechanical behaviour of jointed rock masses (Potyondy and Cundall 2004;Mas Ivars et al. 2011;Sinha and Walton 2020;Castro-Filgueira et al. 2020;Tang et al. 2020). Furthermore, various studies have used the software PFC developed by Itasca (2018) to numerically investigate the shear behaviour of rock joints. For instance, Cundall (2000), Park and Song (2009) and Asadi et al. (2012) used PFC to study the shear behaviour of rough rock joints with triangular, sinusoidal, and standard JRC profiles applying the bond removal method between particles on each side of the joint profile. In the bond removal method, the contact between particles that are within a specific distance on both sides of the joint profile is left unbounded. Although these studies demonstrated the capability of PFC to capture the principal shear behaviour of rough rock joints, the bond removal method suffered from unrealistic high micro-roughness due to the curvature of the particles used (Bahaaddini et al. 2013).
In an attempt to solve the shortcomings of the bond removal method and better reproduce the shear behaviour of rock joints numerically, the smooth-joint contact model (SJCM) proposed by Cundall in 2005(Mas Ivars et al. 2008 was later used and applied to the particles in contact on each side of the joint profile (Bahaaddini et al. 2013;Lambert and Coll 2014;Lazzari et al. 2014). In the SJCM, particles on different sides of the joint profile are allowed to pass through each other instead of moving around one another. Applying the properties of the SJCM, Lambert and Coll (2014) used PFC 3D to numerically study the shear behaviour of an actual granite rock joint whose surface roughness had been previously measured with optical scanning by Grasselli et al. (2002). However, the results from the numerical shear tests performed by Lambert and Coll (2014) were not compared to the actual shear test conducted in the laboratory. Additionally, the three-dimensional surface roughness utilised in the numerical calculations had a resolution of 1.4 by 1.4 mm. This resolution is not sufficient to accurately capture surface roughness at the grain scale (Grasselli 2001;Grasselli and Egger 2003;Tatone and Grasselli 2009). In a similar attempt, Lazzari et al. (2014) used the optical scanning measurement of the surface roughness of an actual rock joint to numerically study its shear behaviour in PFC 2D . They compared the results of the numerical simulation with the actual laboratory direct shear test. However, the comparison made between the numerical and the actual shear tests was not a clear success since interlocking of single particles presented a major problem in their shear test environment. According to Lazzari et al. (2014), it could happen that certain particles, initially located above the joint profile, lied below the joint profile after the application of the normal stress ( n ) used in the numerical simulations. Based on their results, interlocking due to these single particles occurred after a shear displacement of approximately one particle diameter, when the particles changing side were about to create new contacts with neighbour particles. This is the same phenomenon as the one observed by Bahaaddini et al. (2013). Furthermore, the two-dimensional profile that they used in the numerical simulation was based on a simplification of the optical scanning of the actual surface roughness. To overcome the particle interlocking problem, Bahaaddini et al. (2013) proposed the shear box genesis (SBG) approach. This methodology includes the detection of particles situated in the upper and lower parts of the numerical rock specimen and allows for the generation of new SJCM along the joint plane during the numerical shear test. Bahaaddini et al. (2013) used PFC 2D and applied the SBG approach in order to numerically study the shear behaviour of saw-tooth triangular joints and standard JRC profiles on Hawkesbury sandstone. Most recently, similar approaches have been used to numerically study the scale effect of rock joints (Bahaaddini et al. 2014), improve the application of the SJCM (Mehranpour and Kulatilake 2017), study the shear behaviour of rock-concrete surfaces in dam foundations (Gutiérrez-Ch et al. 2018), and study the mechanical behaviour of synthetic rock joint profiles generated with self-affine fractal theory (Ríos-Bayona et al. 2018).
The above review shows that a methodology in PFC to realistically account for the surface roughness of actual rock joints measured with optical scanning and successfully compare numerically simulated and actual laboratory shear tests is still lacking. The only exception is Li et al. (2020), who used PFC 2D to compare a numerically simulated direct shear test with an actual laboratory test conducted on a perfectly mated rock joint sample of granite tested under a constant n of 30 MPa. The two-dimensional joint profile tested by Li et al. (2020) had been selected through the steepest asperity. Although their results showed good agreement between the numerical and actual shear tests, it is uncertain if they are applicable for modelling the shear behaviour of hard crystalline rocks under lower n to joint compressive strength ratios, which are typical in civil engineering applications. Furthermore, the governing asperities in contact under low n , especially in well-mated joints in hard crystalline rocks, are on the order of some tenths of a millimetre (Johansson 2016). Therefore, modelling the strength and size of the active asperities in contact becomes more relevant in PFC to capture the actual shear behaviour of rock joints in these types of rocks. However, it is not clear how the resolution of the optical scanning measurements of hard crystalline rock joints should be combined with the calibration of the different parameters available in PFC to realistically account for their surface roughness and the strength of the asperities in contact in the numerical simulations.
Taking up this challenge, this paper investigates the capability of PFC 2D to realistically simulate the peak shear strength of hard crystalline rock joints under different applied constant n magnitudes. To achieve this aim, this paper presents a methodology in PFC 2D that uses the surface measurements of a rock joint obtained with optical scanning 1 3 at a resolution that captures grain scale detail. Furthermore, to realistically account for the surface roughness and the strength of the active asperities during the shearing process in PFC 2D , the scanning measurements are combined with the calibration of joint segment length, particle resolution per segment, and bond strength between particles.

Numerical shear test environment in PFC 2D
The numerical shear test environment in PFC 2D presented in this paper follows a procedure similar to the specimen generation procedure for the SBG approach suggested by Bahaaddini et al. (2013), which aims to avoid interlocking between particles. However, there are three main differences between the SBG approach described in Bahaaddini et al. (2013) and the shear test environment developed in PFC 2D Version 5 (Itasca 2018) presented in this study. First, only one wall representing the joint profile is inserted to divide the assembly of particles into two parts in the specimen generation step. The SBG approach by Bahaaddini et al. (2013), on the other hand, divides the numerical rock sample into two parts by inserting two walls that follow the joint profile. Second, the elimination of floating particles is not carried out in the numerical shear test environment presented in this paper. Instead, a predefined installation gap ensures that parallel bonds are installed at each particle-particle contact. Furthermore, the elimination of floaters is not available in the PFC version used in this study, since it was replaced by this new approach (i.e., installation gap). Third, the stress measurements during the installation of the isotropic stress are carried out on the whole upper and lower parts of the specimen separately. The mean stress in the methodology described by Bahaaddini et al. (2013) was measured by setting five measurement circles in each part of the rock sample.
The modelling procedure in the shear test environment in PFC 2D presented in this study involves the following four steps: 1. Specimen generation: In the first step, the generation of the rock assembly is carried out in six different substeps, as illustrated in Fig. 1a-f: Fig. 1 Specimen generation procedure in the numerical shear test environment: a vessel generation; b grain assembly; c joint wall generation; d isotropic stress installation; e Parallel bond installation; f specimen removed from vessel 1 3 (a) The material vessel consisting of four walls defining a rectangle in 2D with a specific width and height is generated. The total width of the vessel is determined by adding an extra length of twenty times the maximum particle radius ( R max ) to the actual width of the rock joint sample (see Fig. 1a). (b) In the second sub-step, the vessel is filled with an assembly of particles located in random positions (see Fig. 1b). The size of the particles varies between the range of the defined minimum particle radius ( R min ) and R max satisfying a uniform size distribution. The number of particles depends on the volume of the material vessel, the grain size distribution, and the overall grain porosity ( n c ), respectively. The n c value used in the numerical simulations in PFC 2D aims at obtaining a well packed assembly of particles, and it is not the porosity of the actual rock material. The materialmodelling support in PFC recommends a n c value of 0.08 for simulations in 2D (Potyondy 2018a). (c) The assembly of particles is then divided into two parts, upper and lower, by inserting a wall that follows the joint profile. This is illustrated in Fig. 1c. The particles are allowed to freely rearrange separately in the upper and lower parts under the condition of zero friction, interacting also with the wall that represents the joint profile. The kinetic energy in the numerical simulations is dissipated via the local damping factor available in PFC 2D . This local damping factor is equal to 0.7 to approximate quasi-static conditions (Potyondy 2018a). (d) The specified isotropic stress is applied in the upper and lower parts of the numerical specimen as if they were independent samples using the grain-scaling procedure described in Potyondy (2018a) (see Fig. 1d). The particle sizes are scaled iteratively until the target stress (0.1 MPa) was within tolerance (0.01 Pa) in both parts. Furthermore, the particles in the upper and lower parts, respectively, are identified by assigning a different zone ID. The zone ID is set to 1 for particles in the upper part and to 2 for particles in the lower part. (e) Parallel bonds are installed at each particle-particle contact with an installation gap less than or equal to 30% of the R min value, as illustrated in Fig. 1e. The aim of the installation gap is to ensure that no particles are left without a parallel bond at the end of the specimen generation procedure (i.e., floating particles). (f) Finally, the specimen is removed from the material vessel. The walls defining the vessel and the joint profile are deleted, as illustrated in Fig. 1f.

Smooth-joint contact installation:
In the second step, the joint is inserted in the numerical rock specimen by assigning SJCM to all particle-particle contacts intersected by the joint profile, as illustrated in Fig. 2a. The SJCM installation performed in the shear test environment in this study uses the discrete fracture network (DFN) logic embedded in PFC 2D . Thus, each segment of the rock joint profile is represented by a DFN joint, which is identified by an ID number different to zero. This ID number is used to assign the properties of a specific DFN joint to the SJCM generated by that joint segment. The normal vector of the SJCM assigned to each particle-particle contacts intersected by the joint profile is always directed normal to the joint fragment. For the case of particle-particle contact, the normal vector is directed along the line between their respective centres (Itasca 2018). In the case that the contact between two particles is intersected by two DFN joints, the installed SJCM corresponds the joint segment to the left of the contact. This is, however, a rare case that can be prevented with an appropriate particle resolution ( P res ). 3. Seating phase: In the third step, the specified n is applied to the upper part of the numerical specimen. Additionally, the left and right sides of the upper part of the rock specimen are carved, exposing two wings in the lower part (see Fig. 2a). The value assigned to the width of each carved wing is ten times the R max . The main purpose of the wings is to provide continuity to the joint and to prevent the upper part from bending after large horizontal displacements during the numerical direct shear tests. Furthermore, four walls are installed to provide the boundary conditions to the numerical specimen. The lower sample is fitted in a box defined by three fixed walls and a top wall is installed over the upper part. The n is applied by means of a servo-control mechanism installed on the top wall, making it constant during the numerical shear test. The servomechanism used in the numerical simulations in this study is a procedure available in PFC 2D that provides the ability to control the velocity of the top wall to apply or maintain a certain level of n (Itasca 2018; Potyondy 2018a). During the seating phase, the upper part of the specimen is compressed downwards until the applied n is within the specified tolerance. This process may lead to the generation of new contacts between both parts of the sample at the joint profile location. In this phase, the SJCM's installation is managed by the joint profile location defined by the DFN joint, as explained in step 2. Therefore, since the SJCM's installation is still active in this step, the SJCM is assigned to the new contacts created.

Direct shear test:
In the final step, the specimen is ready to be tested in the numerical shear test environment. Prior to the numerical shear testing, the DFN ID assigned in the SJCM installation is stored in the contact variable 1, thus pointing out the joint segment where the SJCM is located. Furthermore, the ball extra variable 1 is defined to store the ID of the closest DFN-joint to a certain particle, provided that the distance between the particle's perimeter and the DFN-joint is less than a specified tolerance set by default to R max . This is illustrated in Fig. 2a. The ball extra variable 1 defines two layers of particles above and below the joint profiles' boundary that is used to assign SJCM to new created contacts during the numerical shear tests. A new SJCM is installed at a new generated contact during a numerical test if the ball extra variable 1 of any of the particles defining the contact is set to a value different to zero; otherwise, a parallel bond will be assigned to the contact. The assigned properties of the SJCM correspond with the DFN ID stored in the ball extra variable 1 of the particle further away from the driving grip (i.e., where the velocity is applied). To simulate a laboratory direct shear test, two grips are defined at the left and right sides of the upper part of the specimen. The grip thickness is set to 1.5 times the R max . This is illustrated in Fig. 2b. The velocity is gradually applied to the right, or left, Fig. 2 Numerical rock joint specimen with saw-tooth asperities after SJCM installations: a prior the numerical shear test with carved wings and installed walls providing boundary conditions. The zoom in around the saw-tooth joint profile illustrates the two layers of particles defined by the ball extra variable 1 (magenta colour), the SJCM between particles on both sides of the joint profile (green colour), and the particles located further away from the joint profile with installed parallel bonds (blue colour); b during the numerical shear test with defined left and right grips simulating an actual laboratory test grip until a maximum value of 0.05 m/s is reached in the numerical direct shear test. This velocity in the shear test environment in PFC 2D is low enough to ensure that the sample remains in quasi-static equilibrium. Furthermore, the gradual application of shear velocity in the numerical simulations reduces the possibility of having a dynamic compression wave in the upper part of the numerical specimen.
There are, however, still situations, such as when applying high n or after the upper part has undergone a considerable shear displacement (more than one particle diameter), in which the generation method alone is not sufficient with respect to interlocking. This issue has also been discussed by Mehranpour and Kulatilake (2017) and Li et al. (2020). During the initial trial simulations with high applied n , it was observed that particles on both sides of the joint profile with SJCM were pushed against each other, thus generating large stress concentrations. The observed particle interlocking due to overlapping was remedied in this study by increasing the SJCM's normal stiffness ( k n ) to a constant value based on the uniaxial compressive strength of the intact rock ( ci ), given by where R avg is the average particle radius in the numerical specimen.

Intact rock parameters
The intact rock material utilised in the numerical shear tests presented in this study was generated with the materialmodelling package functions provided in PFC 2D using parallel bonded material (PBM) (Potyondy and Cundall 2004; Potyondy 2018a). The micro-properties of the PBM were calibrated against the intact rock properties of a granite rock joint sample previously tested in the laboratory by Johansson (2016). The rock material came from the Flivik quarry in Sweden and consisted of grey coarse-grained granite. The measured macro-properties of the intact granite rock are provided in Table 1.
In PFC 2D , a series of numerical uniaxial compressive strength (UCS) tests was conducted to match the ci , Young's modulus ( E i ) and Poisson's ratio ( i ) of the intact material measured in the laboratory. The numerical specimens generated in the UCS tests had a width of 50 mm and a height of 100 mm (see Fig. 3a). The specimen generation procedure in the numerical UCS tests is similar to the generation procedure of the shear test environment described in section "Numerical shear test environment in PFC 2D ". However, the shear test environment generates the upper and lower parts of the numerical rock specimen as two independent materials. To create the rock assembly, the rectangular vessels (2D) were filled with particles with a R min of 2.88 mm and a n c value of 0.08. The particle size distribution was assigned as uniform with a R max = 1.66 ⋅R min . This relationship gives a R avg of 3.83 mm and a P res of 6 particles per numerical specimen width. This adopted P res value aimed at finding a trade-off between modelling the surface roughness, the asperity strength of the actual rock joint sample, and computational effort. Parallel bonds were also installed at each particle-particle contact with an installation gap less than or equal to 30% of the R min value. To capture the strength variability of the PBM due to different particle packing arrangements, a series of material realisations was created and tested by varying the seed value of random-number generator ( S RN ). Figure 3b illustrates the failure of one of the generated specimens during the UCS tests performed. The results of the deviator stress and radial strain versus axial strain obtained in the numerical UCS tests with different values of S RN are presented in Fig. 4. The average value, standard deviation (SD), and coefficient of variation (COV) of the simulated macro-properties in PFC 2D are presented in Table 1. The calibrated micro-properties of the PBM are provided in Table 2. The parallel bonds with properties in Table 2 installed at each particle-particle contact provide the behaviour of two interfaces (Potyondy 2018a). The first interface is equivalent to the linear elastic model, and the second interface is the parallel bond. When the second interface is bonded, it resists relative rotation and its behaviour is linear elastic until the strength limit is exceeded.

Smooth-joint contact model parameters
During the numerical direct shear tests performed in PFC 2D , the SJCM was assigned to all particle-particle contacts intersected by the modelled rock joint profile. The SJCM parameters utilised in this study are based on results from three laboratory direct shear tests on sawn, planar rock joints performed by Lazzari et al. (2014) corresponding to the basic friction angle. The three laboratory samples had a length of 60 mm, and they were of the same type of granite from the same quarry as utilised by Johansson (2016) in the laboratory direct shear tests. The numerical simulations performed by Lazzari et al. (2014) to calibrate the properties of the SJCM had a R min of 0.38 mm and a R max of 0.62 mm, respectively. The SJCM parameters used in this study correspond to the macro-properties of the joint profile segments, which are independent of particle size. The exception in this study is the value of k n derived with Eq. (1) to prevent interlocking between particles. The values of the calibrated parameters of the SJCM utilised in the numerical simulations in PFC 2D are presented in Table 3. The friction coefficient of the SJCM is assumed to be constant in the numerical direct shear tests performed in this study.

Verification of the shear test environment in PFC 2D with idealised, saw-tooth triangular joint profiles
To check the applicability of the shear test environment in PFC 2D presented in this study to simulate the shear strength of rock joints, a series of numerical direct shear tests on rock joints with idealised, saw-tooth triangular asperities was conducted under constant normal load (CNL) conditions with different applied n . Furthermore, the peak shear strength of the idealised joint profiles obtained in the numerical calculations in PFC 2D was compared with the peak shear strength obtained analytically using the criteria developed by Patton (1966) and Ladanyi and Archambault (1969). These two criteria were originally developed by studying the mechanical behaviour of rock joints with idealised, sawtooth asperities.
The numerical specimens created in PFC 2D had a width of 100 mm and a height of 40 mm. Three different saw-tooth triangular joint profiles with a dip angle ( ) of 15°, 25° and 35° were tested under different applied n (1, 2, 3, 5 and 10 MPa). These values of are similar to the ones utilised in the study by Bahaaddini et al. (2013) and were chosen to capture both sliding over and shearing through the asperities. The base length of the saw-tooth asperities in the numerical specimens was 20 mm. Figure 2 shows one of the tested numerical specimens with saw-tooth asperities and a of 25°. The rock assembly was created by filling the rectangular vessels (2D) with particles with a R avg of 0.8 mm and a n c of 0.08. The particle packing arrangement of all generated numerical specimens was created with the same S RN = 10,001. As in the UCS tests performed, the relationship between R min and R max was 1.66. On average, each asperity of the saw-tooth Fig. 3 Numerical specimen generated in PFC 2D to conduct the UCS tests: a before the test starts; b after reaching the point of failure in the specimen. The red and black colours indicate tensile and shear failure, respectively joint profile had a P res of 6 particles. This is similar to the P res of the generated numerical specimens in the UCS tests for the calibration of the micro-properties of the PBM. In the calibration process shown in section "Calibration of the intact rock and smooth-joint contact model parameters", it is the number of particles (i.e., P res ) perpendicular to the application of the major principal stress that governs the simulated macro-properties of the numerical UCS tests, rather than the particle radius. The micro-parameters of the PBM and installed SJCM are provided in Tables 2 and 3. The results of shear stress ( ) versus shear displacement ( s ) for all the direct shear tests performed in PFC 2D are illustrated in Fig. 5a-c.
The comparison made between and s shows that for all the tested saw-tooth profiles the peak shear stress ( p ) increases with higher values of n . The values of shear displacement at peak ( p ) also increase with higher values of n . Furthermore, the results of the direct shear tests performed show that sliding failure controls the shearing process. The exception is the saw-tooth joint profile with = 35°. In this joint profile, a higher number of tensile and shear cracks could be observed on the surface of the triangular asperities in the simulated shear tests with n higher than 1 MPa (see Fig. 6a-c). This observed mechanical behaviour in the numerical direct shear tests on rock joints with idealised, saw-tooth asperities is consistent with other results reported in the literature by Bahaaddini et al. (2013) and Li et al. (2020).
The values of p obtained in the numerical simulations under the different applied n were then compared to the p Fig. 4 Results of the numerical UCS tests on the numerical specimens generated with different values of S RN : a deviator stress vs. axial strain; b radial strain vs. axial strain 1 3 obtained analytically with the shear strength criteria developed by Patton (1966), and Ladanyi and Archambault (1969). Based on the observations made in Figs. 5a-c and 6a-c on the predominant failure mode controlling the shearing process, only the part of Patton's failure criterion accounting for sliding failure was used. The shear strength criteria of Patton (1966) and Ladanyi and Archambault (1969) can be expressed as: and where b is the basic friction angle, i is the inclination of the saw-tooth asperities with respect to the shear direction, a s is the shear area ratio and ̇v is the rate of dilation at failure. The value of b of the tested granite rock material is 31.4°. The comparison between the p simulated in PFC 2D and that calculated with Eqs. (2) and (3) is illustrated in Fig. 5df. The values of p simulated in PFC 2D were in good agreement with the p calculated with Patton (1966) and Ladanyi and Archambault (1969). This shows that the shear test environment in PFC 2D presented in this study can simulate the mechanical behaviour of rock joints with idealised, sawtooth asperities with applied values of n up to 10 MPa. The reason for the observed discrepancies between predicted p with the criteria of Patton (1966) and Ladanyi and Archambault (1969) in Fig. 5d and e, respectively, is that sliding failure controls the shearing process for the tested rock material at low values. The criterion by Ladanyi and Archambault (1969) accounts for the contribution from both sliding and shearing based on energy principles. However, at higher inclinations (i.e., = 35° in Fig. 5f), shearing through the asperities starts to occur. The contribution from shearing through the asperities is taken into account in Ladanyi and Archambault's peak shear strength criterion, and therefore, the simulated p has a better agreement with the predicted p using Eq. (3) for this value of .

Rock joint sample and joint profile geometry
The shear test environment in PFC 2D was utilised to simulate the peak shear strength of an actual rock joint obtained in laboratory testing. A tensile-induced rock joint sample in grey coarse-grained granite with dimensions of 60 by 60 mm, previously tested in the laboratory under CNL conditions by Johansson (2016), was utilised in this study. A comparison was made between the results of the laboratory and numerical direct shear tests performed with an applied n of 1 MPa. The numerical direct shear tests in PFC 2D were performed on three different two-dimensional profiles selected along the shear direction in the actual rock joint. These joint profiles were selected after studying the high-resolution optical scanning measurements of the surface roughness. The scanning measurements were taken with an ATOS III system and are presented in Johansson (2016). The measurements had a resolution of approximately 120-140 µm. The scanned surface was regenerated with a resolution of 0.3 by 0.3 mm. This selected resolution was based on the observations made in the laboratory by Johansson (2016) on the mechanical behaviour of the perfectly mated rock joint samples tested. He observed that, on average, the measured value of p in these samples was approximately 0.3 mm. He concluded that this gives an indication of the size of the effective asperities contributing to the shear strength at the peak. Furthermore, the resolution of 0.3 by 0.3 mm is assumed to be appropriate to capture the grain size according to previous recommendations by Grasselli and Egger (2003), and Tatone and Grasselli (2009). The first joint profile (2D) was selected in the centre of the analysed sample. The second and third profiles were selected at approximately the centre of the left and right remaining halves, respectively. In addition, profile 3 was taken along the asperity with maximum value of apparent dip angle facing the shear direction ( * ). The digitised rock joint surface in three dimensions and the three two-dimensional profiles tested in the numerical shear test environment in PFC 2D are illustrated in Fig. 7. To be consistent with the laboratory observations from Johansson (2016) and to achieve a good representation of the surface roughness in PFC 2D , the joint profiles in the numerical direct shear tests had a segment length ( SL ) of 0.3 mm. The shear direction indicated in Fig. 7 refers to the shear direction of the lower parts of the tested numerical rock specimens.

Shear test simulations
The numerical specimens created in PFC 2D to test the three different joint profiles had a width of 60 mm and a height of 15 mm. The relationship between height and width (1:4) was based on a preliminary sensitivity analysis in the shear test environment intended to find the optimal sample size that numerically simulates stable shear behaviour without cracking the samples and at the same time be calculated within a reasonable computational time.
The roughness resolution in the simulated rock joint in PFC 2D depends on the magnitude of SL and P res in each segment. In the numerical shear tests performed, the rectangular vessels were filled with particles with a R avg of 0.08 mm. The relationship between R min and R max was 1.66. The particle packing arrangement of the three numerical specimens was created with the same S RN = 10001. On average, each segment of the tested joint profiles had a P res of approximately 2 particles. This P res differs from the one utilised in both the UCS tests in section "Calibration of the intact rock and smooth-joint contact model parameters" and the direct shear tests with the saw-tooth profiles in section "Verification of the shear test environment in PFC 2D with idealised, saw-tooth triangular joint profiles". This was mainly due to limitations in the computational time. To further reduce the numerical instabilities when using small particle radius, a scaling algorithm was used where the particle radius and the dimensions of the numerical shear test were scaled with a factor of 20. The generation and direct shear test procedures in PFC 2D of the numerical specimens follow the principles explained in section "Numerical shear test environment in PFC 2D ". Furthermore, the upper parts of the numerical specimens were sheared in the same direction as the laboratory shear tests (from left to right) by gradually applying velocity on the left grip (see Fig. 8). The n applied during the numerical direct shear tests was equal to 1 MPa, as in the laboratory test. The micro-parameters of the PBM and properties of the SJCM installed in the numerical shear tests are provided in Tables 2 and 3. A comparison between the results from the laboratory direct shear test and the numerical simulations in PFC 2D on the selected joint profiles in the form of mobilised friction angle ( ) versus s is presented in Fig. 9. The values of are calculated from the recorded values of and n during the numerical simulations in PFC 2D , which is given by The measured mobilised peak friction angle ( p ) of the rock joint sample tested in the laboratory by Johansson (2016) was 65°. The results of the numerical direct shear tests with profiles 1, 2 and 3 showed a p of 52.8°, 58.2° and 65.8°, respectively. The p of profile 3 simulated in PFC 2D was 13° and 7.6° higher than the p of profiles 1 and 2, respectively. The difference between the p measured in the laboratory and that simulated in PFC 2D , as expressed in absolute values, varied between 0.8° and 12.2°. The measured post-peak in the laboratory test after a s of 1.5 mm was 52°. The postpeak of profiles 1, 2, and 3 simulated in PFC 2D after a s of 1.5 mm was 24°, 38°, and 33°, respectively. The difference between measured post-peak in the laboratory and simulated in PFC 2D varied between 14° and 28°. The measured p (4) = arctan n . Fig. 5 Results of the numerical shear tests on the rock joint specimens with idealised, saw-tooth asperities with different applied values of n , and comparison with the analytical criteria by Patton (1966), and Ladanyi and Archambault (1969): shear stress, vs. shear displacement, s with a = 15°; b = 25°; c = 35°; and peak shear stress, p vs. normal stress, n with d = 15°; e = 25°; f = 35°◂ Fig. 6 Asperity degradation during the numerical direct shear tests in PFC 2D on the rock joint samples with saw-tooth asperities with = 35° after a s of 5 mm: a n = 1 MPa; b n = 3 MPa; c n = 10 MPa. The red and black colours indicate tensile and shear failure, respectively of the laboratory test was 0.35 mm. The values of p in the numerical simulations with profiles 1, 2, and 3 were 0.41, 0.18, and 0.15 mm, respectively.
The numerical direct shear tests in PFC 2D performed on profiles 1, 2, and 3 after a s of 1.5 mm is illustrated in Fig. 8. The results of the numerical shear tests showed that sliding failure mostly controlled the shearing process in PFC 2D . There were, however, some areas along the tested joint profiles where tensile and shear cracking could be observed. These cracks developed superficially and were located in those areas where the asperities are steeper, i.e., higher values of * . These results from the numerical simulations agree with the observations made by Johansson (2016) on the actual rock joint tested in the laboratory.

Numerical shear tests in PFC 2D with actual rock joint profiles
The obtained results of the numerical direct shear tests performed in PFC 2D with profiles 1, 2, and 3 show that the shear test environment, together with the methodology presented in this study, has the capability of realistically capturing the principal shear behaviour of actual hard crystalline rock joints under the applied value of n .
From a qualitative point of view, the three joint profiles captured well the peak and post-peak behaviour compared with the actual analysed rock joint. Furthermore, the results obtained in the numerical shear tests indicate that the effect of increasing the k n of the SJCM as in Eq. (1) has a negligible influence on the simulated peak shear strength of the tested rock joint profiles. However, the elastic part of the stress-shear displacement curve could be slightly influenced by this assumption, due to the contribution of the local normal stiffness of the asperities to the overall shear stiffness of the sample (see Figs. 5 and 9). From a quantitative point of view, only profile 3, which was taken through the steepest asperity, was in reasonably good agreement when comparing its p simulated in PFC 2D with that measured in the laboratory. However, it is still uncertain if the number of contact points generated in a joint profile in 2D containing the steepest asperity can capture the three-dimensional characteristics of surface roughness of actual rock joints. A three-dimensional approach may therefore capture the actual contact area more accurately. However, the computational time needed to conduct a numerical direct shear test in 3D would increase significantly. A discrepancy was observed when comparing the postpeak simulated in PFC 2D to the residual behaviour measured in the actual rock joint tested in the laboratory. The simulated post-peak in the numerical shear tests on profiles 1, 2, and 3 after a s of 1.5 mm was, on average, 31.7°. After the same s , the measured post-peak in the rock joint tested in the laboratory was 52°. A possible reason for this discrepancy between the simulated and measured post-peak may be due to the fact that the SJCM installed in the numerical simulations performed in PFC 2D was set to small strain mode. This means that one particle in PFC 2D may have more than one active SJCM after the upper part of the numerical specimen has been sheared more than one particle diameter. This implies that the shear behaviour after large displacements (more than one particle diameter) may not be totally accurate. This may contribute to the observed discrepancy between the post-peak simulated in PFC 2D and that measured in the actual rock joint tested in the laboratory. Although the large strain mode can be assigned to the installed SJCM, the validity of the results obtained by assigning this mode has not been proven in this study. Another possible reason for this discrepancy may be due to fewer contact points generated along the tested joint profiles in PFC 2D in comparison with the actual rock joint in 3D. Additionally, the upper part of the numerical specimen in PFC 2D was always sheared horizontally, and it was not able to tilt and accommodate during the numerical shear tests. This further contributes to fewer contact points compared with the actual shear test. This approach where the upper part is not able to tilt and accommodate has also been applied in previous research that uses PFC 2D to conduct numerical direct shear tests (Bahaaddini et al. 2013(Bahaaddini et al. , 2014. A three-dimensional approach where the upper part of the sample can tilt and accommodate may therefore capture rock shear behaviour more accurately. A second discrepancy was observed when comparing the p in the numerical direct shear tests of the two-dimensional joint profiles in PFC 2D to that measured in the laboratory on the actual rock joint. The value of p simulated in PFC 2D on the three joint profiles was, on average, 0.25 mm. The measured p of the actual rock joint in the laboratory was 0.35 mm. The possible reason for this discrepancy between the simulated and measured p may be due to the upper and lower parts of the actual rock joint having a small mismatch prior the laboratory test. This contributes to a less steep strength-displacement curve in the beginning of the shear test and a larger measured p since this depends on ig. 9 Results of mobilised friction angle, vs. shear displacement, s of the laboratory direct shear test on the 60 by 60 mm rock joint sample tested by Johansson (2016) and the numerical direct shear tests in PFC 2D with profiles 1, 2 and 3 under a n of 1 MPa the shear stiffness of the asperities in contact during the shearing process (see Fig. 9). On the contrary, the upper and lower parts of the numerical specimens in PFC 2D had particle-particle contact along the joint profile before the direct shear tests were simulated. Furthermore, the slope of the strength-displacement curve in the numerical model in PFC 2D depends on both the shear and normal stiffness of the installed SJCM, and the stiffness of the PBM installed between the particles modelling the asperities. This discrepancy in the slope of the strength-displacement curve between shear tests performed in PFC 2D and in the laboratory has also been observed in previous studies (Li et al. 2020). This discrepancy may be reduced using a lower value of the shear stiffness of the SJCM assigned to the particles intersected by the simulated joint profile, to compensate for the high k n value used to prevent particle interlocking derived using Eq. (1).

Influence of surface roughness resolution in the simulated shear strength in PFC 2D
To study the influence of the surface roughness resolution used in the joint profiles in PFC 2D on their simulated shear behaviour, a sensitivity analysis was performed on the SL and the P res . Due to the relatively long computational time, this sensitivity analysis only included profile 3, which showed better agreement between the simulated p in PFC 2D and that measured in the laboratory (see Fig. 9). This sensitivity analysis included four additional numerical direct shear tests under a constant n of 1 MPa. The micro-parameters of the PBM and the properties of the SJCM used in the numerical simulations are provided in Tables 2 and 3. Two tests were first performed in PFC 2D using values of SL of 0.6 and 1.2 mm, respectively. The value of R avg was 0.08 mm, as in section "Numerical direct shear test in PFC 2D of an actual rock joint in granite". In these two numerical direct shear tests, the value of P res increased from 2 to 3.4 and 6.9 particles, respectively. The results of versus s for these two numerical shear tests are illustrated in Fig. 10a. The other two additional tests were performed with value of SL of 0.3 mm (i.e., as in section "Numerical direct shear test in PFC 2D of an actual rock joint in granite") and values of R avg of 0.18 and 0.27 mm, respectively. The value of P res in these numerical tests decreased from 2 to 0.75 and 0.5 particles, respectively. The results of versus s for these two numerical shear tests are illustrated in Fig. 10b.
The results illustrated in Fig. 10a show that the simulated p in PFC 2D with a R avg of 0.08 mm and values of SL of 0.6 and 1.2 mm were 61.5° and 61.1°, respectively. The values of p in the numerical simulations were 0.33 and 0.25 mm, respectively. These values of p are lower than both the simulated p and the measured p of the actual rock joint (see Fig. 9). Furthermore, the numerical simulations in Fig. 10a had a larger p and a smoother transition between peak and post-peak behaviour than the numerical simulation in Fig. 9. These results are not surprising, since part of the joint surface information captured in the measurements with optical scanning is missed when SL values larger than the grain size are used in the numerical simulations. Therefore, the results presented in Fig. 10a show the importance of using a value of SL in PFC 2D that realistically captures the surface roughness of the actual rock joint at grain scale.
The results illustrated in Fig. 10b show that the simulated p with a SL of 0.3 mm and values of R avg of 0.18 and 0.27 mm were 70.7° and 61.6°, respectively. The values of p in the numerical simulations were 0.32 and 0.4 mm, respectively. These values of p are both higher and lower than the simulated p shown in Fig. 9 using the same joint profile. The values of p in the two numerical simulations in Fig. 10b were also larger than the p obtained in Fig. 9 using the same joint profile. Furthermore, the match between the p obtained in these two additional simulations in PFC 2D and the p measured in the actual rock joint was not good. A possible reason for this mechanical behaviour observed in PFC 2D is that using a value of SL that captures surface roughness at grain scale but with an overly large particle size is incorrect. This means that filling the simulated SL in PFC 2D with sufficient P res is necessary to achieve a good representation of the surface roughness of an actual rock joint at grain scale. The results obtained in the numerical simulation with profile 3 using a SL of 0.3 mm and a R avg of 0.08 mm show that a P res of at least 2 particles per joint segment captures the mechanical behaviour of the actual rock joint. Thus, using SL of 0.3 mm and a P res of at least 2 particles per joint segment yields results in line with the observations made in the laboratory by Johansson (2016) on the actual granite samples. Therefore, prior to numerical shear testing of perfectly mated rock joints using PFC 2D , it is recommended that a SL of 0.3 mm and a P res of at least 2 particles per joint segment should be used.

Influence of the calibrated micro-properties and particle packing arrangement on the simulated shear strength in PFC 2D
To study the influence of the calibrated micro-properties of the PBM and the particle packing arrangement on the simulated shear strength in PFC 2D , a sensitivity analysis using profile 3 was performed. This sensitivity analysis included two additional series of numerical direct shear tests performed under a n of 1 MPa using two different calibrated micro-properties of the PBM and five different values of S RN . The values of SL and R avg utilised were 0.3 and 0.08 mm, respectively (as in section "Numerical direct shear test in PFC 2D of an actual rock joint in granite"). The first series of numerical direct shear tests was performed using the micro-parameters of the PBM and properties of the installed SJCM provided in Tables 2 and 3. The results of versus s for this series of numerical simulations using five different values of S RN are illustrated in Fig. 11a. The numerical shear test with S RN = 10,001 is the same as the one illustrated in Fig. 9.   Fig. 11 Results of mobilised friction angle, vs. shear displacement, s of profile 3 under a n of 1 MPa and using five different values of S RN in the specimen generation in PFC 2D : a with calibrated micro-properties of the PBM in the numerical UCS tests using a P res of 6 particles; b with calibrated micro-properties of the PBM in the numerical UCS tests using a P res of 20 particles The second series of numerical direct shear tests was performed using different micro-properties of the PBM. A new intact rock calibration process following the steps described in section "Calibration of the intact rock and smooth-joint contact model parameters" was conducted. The main difference between the new calibration and that presented in section "Calibration of the intact rock and smoothjoint contact model parameters" is that the 50 by 100 mm numerical specimens were filled with particles having values of R avg = 1.15 mm, n c = 0.08, and P res = 20 particles. Parallel bonds were installed at each particle-particle contact with an installation gap less than or equal to 30% of the R min value. The average values of ci , E i , and i simulated in the new calibration were 195.9 MPa, 84.4 GPa, and 0.26, respectively. These simulated values are in good agreement with the macro-properties measured for the intact granite rock in the laboratory (see Table 1). The simulated values of ci , E i and i had a SD of 15.9 MPa, 1.42 GPa, and 0.006, respectively. This gave a COV of 8.1% for ci , 1.7% for E i and 2.3% for i . The calibrated micro-properties of the PBM obtained with the new calibration are provided in Table 4. The results of versus s using the micro-parameters of the PBM in Table 4, and the properties of the installed SJCM provided in Table 3 with five different values of S RN are illustrated in Fig. 11b.
The results illustrated in Fig. 11a show that the series of numerical shear tests performed with the micro-properties in Table 2 had, on average, a simulated p in PFC 2D of 65.5°. The difference between the average value of the simulated p and that measured in the actual rock joint tested in the laboratory was 0.5°. The numerical shear tests with different S RN had a SD of 0.2°. On the other hand, the results illustrated in Fig. 11b show that the series of numerical shear tests performed with the micro-properties in Table 4, which were lower than the micro-properties in Table 2, had, on average, a simulated p in PFC 2D of 63.7°. Compared to the p measured in the laboratory test, the average value of simulated p in Fig. 11b was 1.3° lower. Furthermore, the SD in this second series of numerical shear tests was 1.9°.
The results illustrated Fig. 11a and b show that the average values of p simulated with the micro-properties presented in Tables 2 and 4, respectively, and different values of S RN are in good agreement with the measured p of the actual rock joint. A possible reason for this observation is that under low applied n , the predominant failure mode governing the shearing process is sliding along the active asperities. This can be observed in Fig. 8, where only some parts of the tested profiles developed tensile cracks superficially. For this reason, under a n of 1 MPa, the difference between p simulated in PFC 2D using different micro-parameters of the PBM was not significantly large.
To increase the completeness of this sensitivity analysis, a third series of direct shear tests was performed in which the applied n was increased to 5 and 10 MPa, respectively. The numerical specimens were generated using both the simulated micro-properties of the PBM in Tables 2 and 4, respectively. For simplicity, only two different values of S RN were used. The values of p simulated in PFC 2D were compared with the p extrapolated for higher applied n using the backcalculated value of JRC of the actual rock joint tested in the laboratory under a n of 1 MPa (Barton and Choubey 1977). The value of JRC back-calculated based on the laboratory test of the actual rock joint was 14.9. The results of p versus n are illustrated in Fig. 12. The results show that under an applied n of 5 and 10 MPa, the values of p simulated in the numerical specimens using the micro-properties of the PBM in Table 4 were, on average, 2.1° and 3.8° lower than the extrapolated values of p with a JRC of 14.9, respectively. Compared with the extrapolated values of p using a JRC of 14.9, the values of p simulated in the numerical specimens with micro-properties of the PBM in Table 2 were, on average, 1.0° and 1.2° lower, respectively. The values of p simulated in the numerical specimens with micro-properties of the PBM in Table 2 were in better agreement with the extrapolated values p for the actual rock joint. The reason for this is that the P res during calibration of PBM microproperties was 6 particles for properties in Table 2, and 20 particles for properties in Table 4, respectively. Thus, the micro-properties in Table 4 obtained with a P res of 20 particles overestimate the number of particles that simulates the asperities in contact during the shearing process with profile 3, leading to an overall softer behaviour of the intact rock.
The results presented in this sensitivity analysis show the importance of understanding well how the number and size of the contact points generated in the joint profile during a numerical direct shear test in PFC 2D contribute to simulating the peak shear strength. Additionally, the results presented in Fig. 12 show that the shear test environment in  [-] 5.0 Normal to shear stiffness ratio [-] 5.0 PFC 2D , which has been validated against a laboratory shear test under a n of 1 MPa, may be able to simulate the peak shear strength of actual rock joint profiles (2D) up to a n of 10 MPa. However, it is possible that above this magnitude of n , asperity breakage increases, and it may be necessary to use flat-jointed bonded-particle material to simulate the mechanical behaviour of the intact rock under such a high level of confinement (Potyondy 2018b). The main advantage of the flat-joint model is that it allows for the calibration of the bonded material against the compressive to tensile strength ration. This is an ability that parallel-bonded material lacks. However, the used of flat-jointed material requires also more computational power.

Conclusions
In this paper, we investigate the capability of PFC 2D to realistically simulate the shear strength of hard crystalline rock joints under different applied magnitudes of n with surface roughness obtained from optical scanning measurements. Based on the simulations performed with two-dimensional profiles selected from the digitised rock joint surface, it can be concluded that the numerical shear test environment in PFC 2D presented in this paper has the capability of capturing the principal shear behaviour of actual rock joints under the conditions tested in this study.
The results of the simulations performed in the shear test environment in PFC 2D show that it is possible to match the peak shear strength of the actual rock joint of granite tested in the laboratory by selecting the two-dimensional profile in the shear direction that contains the steepest asperity. The results of the sensitivity analysis performed in PFC 2D show that to accurately capture the surface roughness of actual rock joints, a magnitude of joint segment length that captures grain size and at least two particles per segment are both needed. Additionally, the sensitivity analysis performed with different properties of the PBM shows that weaker bonds have a greater influence on the simulated peak shear strength with applied normal stresses higher than 10 MPa. This highlights the importance of adapting the calibration of the bond properties between particles in PFC 2D to simulate the size and failure mechanism of the active asperities in contact during the shearing process.
The post-peak shear strength obtained in the numerical simulations on the joint profiles was lower than that measured in the laboratory. The use of a small strain mode during the numerical simulations in the shear test environment in PFC 2D may have an influence on the simulated residual behaviour. Furthermore, the 2D approach together with the inability of the upper part of the samples to freely rotate in the PFC 2D environment generates fewer contact points compared to the 3D approach. In order to be able to realistically model the complete shear behaviour of rock joints in the future (i.e., including residual shear strength), a threedimensional approach is recommended. In addition, the shear box should have the ability to let the upper part of the sample freely rotate in the same way as in actual shear tests. By doing so, a more accurate number and size for the contact points will be obtained.
The methodology in PFC 2D presented in this study has been tested to simulate the shear behaviour of an actual rock joint sample of granite to realistically account for its surface roughness and asperity strength. However, the segment length value required to accurately capture grain size may be different for other rock types. Both experimental work and numerical simulations are recommended to investigate the combination of scanning measurements and calibration of PFC parameters in order to simulate the shear behaviour of rock joints in other type of rocks. The practical implications of this methodology are that surface scanning measurements may be used to simulate the shear behaviour in 2D using the joint profile containing the steepest asperity. However, the presented methodology has only been used to simulate the shear behaviour of an actual rock joint sample with perfect match between its contact surfaces. Further studies are required to study the applicability of numerical techniques to simulate the shear behaviour of actual rock joints without perfect match between the contact surfaces. Extrapolation with JRC=14.9 PBM in Table 2, S RN =10001 PBM in Table 2, S RN =10005 PBM in Table 4, S RN =10001 PBM in Table 4, S RN =10005

Fig. 12
Results of p under different applied n simulated in PFC 2D with profile 3 using two different values of S RN , calibrated microproperties of the PBM with a P res of 6 and 20 particles, respectively, and its comparison with the p extrapolated by back-calculating the value of JRC of the actual rock joint, as proposed by Barton and Choubey (1977), based on the laboratory test with a n of 1 MPa Funding Open access funding provided by KTH Royal Institute of Technology. The research presented in this paper was supported by the Rock Engineering Research Foundation (BeFo) (Grant No.364) and the Swedish Hydropower Centre (SVC) (Grant No.VK10798). SVC was established by the Swedish Energy Agency, Elforsk and Svenska Kraftnät together with Luleå University, KTH Royal Institute of Technology, Chalmers University of Technology and Uppsala University. http:// www. svc. nu. The authors acknowledge the support of the Swedish Nuclear Fuel and Waste Management Company (SKB) (order number 22483). The authors also acknowledge the support of Itasca Consultants AB for providing with a PFC 2D license for the development of the numerical shear test environment.

Declarations
Competing interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.