Simulation of direct shear tests using a forces on fracture surfaces (FFS) approach

The aim of this work is to provide a complete data set of direct shear tests and to propose a corresponding simulation approach. Tests have been conducted on crystalline rock samples applying constant normal load (CNL) and constant normal stiffness (CNS) boundary conditions. A physical consistent algorithm which explicitly calculates the forces acting on the fracture surface (FFS) has been developed. This FFS approach can explain the occurrence of surface degradation and shows the main shear characteristics. After all, shearing of rough rock joints remains a complex process and the differences between laboratory and simulation results are still significant in some cases. All data and input files are provided free for download and testing.


Motivation
The shear characteristics of rock joints are important in geotechnical engineering. Shear forces act on most discontinuities due to gravitation, tectonic forces or stress redistributions. Shear movements along discontinuities such as joints or fracture planes can cause substantial damage or even the collapse of engineering structures at the surface as well as underground. If shear stresses on structures exceed the local shear strength, shear movements will arise which typically are coupled to normal displacements. This is commonly observed as joint dilation or closure and will cause a change in permeability. In case of CNS boundary conditions, a normal stress feedback is generated. This complex behavior is highly influenced by the roughness and local strength of the joint surfaces. Research on this topic started several decades ago and has not been finished yet. The aim of the presented study is to deepen the understanding of shear characteristics of rough rock joints and to develop a new constitutive law to simulate laboratory tests.

State of the art: (semi-) analytical solutions
The general form of a shear process, which is based on friction laws, follows the equation: where τ is the shear stress, σ n the normal stress acting on the surface, ϕ b the basis friction angle of the smooth material and roughness an additional shear strength parameter caused by the unevenness of the material.
This relation was adapted my numerous researchers. Barton and Choubey (1977) developed an easy and intuitivebut not quantitative-criterion to describe the roughness: the joint roughness coefficient (JRC): where JCS is the joint compressive strength-a measure of normal deformation to normal forces. Later, researchers like Tse and Cruden (1979) developed quantitative calculations of the JRC.
(1) = n ⋅ tan b + roughness This article is a part of the Topical Collection in Environmental Earth Sciences on "Sustainable Utilization of Geosystems" guest edited by Ulf Hünken, Peter Dietrich and Olaf Kolditz.
Later, Barton et al. (1985) introduced a constitutive model providing shear stress-shear displacement-dilation as well as normal stress-joint closure coupling. In their model, joint deformation and strength development is based on the JRC mobilization concept. This approach formulates the joint roughness mobilized or destroyed pre-peak or post-peak, respectively: The shear strength which is mobilized at any given shear displacement δ is denoted as φ(mobilized). It depends on the corresponding magnitude of JRC(mobilized), compressive strength capacity of the joint surface as well as on the residual friction angle φ r . Typically, joint roughness parameters available are valid for peak shear strength conditions and, therefore, denoted as JRC(peak). Since the dimensionless numbers JRC(mobilized)/JRC(peak) and δ(mobilized)/δ(peak) are found to increase almost identically during shear, standard tables of values of these dimensionless terms can be utilized to predict the shear stress-displacement behavior (Barton et al. 1985). Grasselli and Egger (2003) used a set of quantitative roughness parameters * max , C which are shear direction dependent but quite academic: Xia et al. (2014) also used Grasselli's parameters in their equation: Other researchers using Grasselli's parameters are Grasselli (2012, 2013) and Casagrande et al. (2017). Despite its sophisticated approach these parameters are perfect to combine with surface scan data which might be an explanation for their popularity. Zhang et al. (2016) used the contact ratio and the first order deviation of the surface facing the shear direction: Yang et al. (2016) combines Grasselli's and Barton's parameters: (3) (mobilized) = JRC(mobilized) × log JCS n + r . (4) (5) (7) An overview of 25 criteria including the ones given above and all used symbols can be found in Singh and Basu (2018).
Nevertheless, closed form solutions are necessary to calculate bigger models containing numerous rock joints, where the explicit geometry is unknown. The model by Barton et al. (1985) works with a set of functions [where (2) is the basis] which offers a good relation between complexity/calculation time versus accuracy. Li et al. (2018) and Gui et al. (2017) are using sets of functions, too. Their approach allows considering more physical details of the shear process.

State of the art: numerical methods
Another possibility to calculate the shear behavior of rock joints is to use the actual surface morphology and to simulate the shear process directly using this information. Surface scan data are nowadays widely used as a 3D-representation of the joint morphology. In a first step profile lines are scanned and an artificial average profile is generated (Fan and Cao 2020). More advanced approaches use the entire surface. Different numerical continuum and discontinuum based methods can handle the surface geometry and allow to backanalyse shear experiments, see for instance Nguyen (2013), Nguyen et al. (2014), Saadat and Taheri (2020a, b), Varela Valdez et al. (2018), Nemcik et al. (2014), Dang et al. (2019) or Jiang et al. (2020).
It is also possible to just use the geometric information of the surface for the calculations. The scan data have to be transformed in a grid structure, for example based on triangles or squares. The regions of the two surface parts in contact are of interest. The advantage of this FFS (Forces on Fracture Surfaces) approaches is the compactness. It fills the gap between complex numerical codes and simple analytical solutions. The analysis is done on a numerical grid. This allows full control of the simulation and focusing on one specific part of the problem. Examples for this kind of approach can be found in Fathi et al. (2016a, b), Park and Song (2013) and Casagrande et al. (2017). The basis of the new code presented in this paper is based on the algorithm by Casagrande et al. (2017). This approach calculates the shear forces on a triangulated surface and compares the forces which are needed to slide over a surface element: with the force needed to shear-off a surface element: If the sliding force is bigger than local shear strength, the surface element will be sheared-off and the calculation is repeated with the new geometry. This procedure will be iteratively repeated, until no further shear-off is recognized. The physical basis of this approach is solid. The code can just calculate peak and residual shear stresses but no shear stress-shear displacement relations and no dilatation.

Required functionality
The code by Casagrande et al. (2017) was used as a basis. The framework of the proposed procedure is based on the following functionalities: • Keep it physical motivated, no fitting parameters • Geometry of surface, basic rock parameters and boundary conditions of laboratory test as input • Calculation of whole shear stress-displacements curves • Calculation of dilatation • Calculation of CNS tests (CNL tests are obtained by setting the normal stiffness of the loading system K = 0, refer to Fig. 3) • Graphical user interface • Implementation for standard PC • Moderate simulation time (about 1 h)

Calculation scheme
The scheme of the proposed code is shown in Fig. 1. Basic input parameters are the 3D-coordinates of the points forming the shear plane as determined by the surface scanner. The surface geometry is stored as a matrix. A quadratic grid is used which is created from the raw scan data by a linear interpolation to the neighboring points. The surface is then duplicated and the shearing is done by shifting the matrices against each other by one pixel. Thus, the grid constant is the increment of the shear displacement. Evolution of shear stress and dilation with shear displacement is calculated in an incremental, stepwise manner until the user-set limit for maximum shear displacement is reached. At each shear step forces and new geometry are calculated. The formulas for the shear forces and sliding forces are (8) and (9), respectively. The normal force acting on the surface has a counter force due to the surface deformation. An elastic deformation is assumed: where a is the grid constant, E the elastic modulus, Δh i the deformation of a surface element and h 0 the height of the surface element-which is assumed to be the sample height. The height of the used samples is about h 0 ≈ 15 cm which is

Direct shear tests: lab testing
The direct shear tests were conducted in the rock mechanical laboratory of the TU Bergakademie Freiberg, Germany. To conduct the direct shear tests a special designed shear box device (Konietzky et al. 2012) was used. Figure 2 shows the used apparatus and Table 1 gives the key parameters.
A schematic visualization of a direct shear test is shown in Fig. 3. Two series of shear tests have been conducted in the framework of the GeomInt project: CNL tests using a granite specimen and CNS tests using a basalt specimen. The basic rock parameters are given in Table 2. Main results of the two experimental series are illustrated in Fig. 4. The CNL test at the first normal stress level shows a clear peak in the stress-displacement curve. After reaching the final shear displacement, the opposite joint faces were reset to their initial (pre-test) position under zero normal stress to avoid surface wear. Repeated shearing of the same sample with increasing normal stress results in shear curves which immediately reach the residual strength plateaus without producing peak values. For the CNS test the first step was performed as CNL test. Similar to the test mentioned before, a peak value with subsequent shear softening to a residual plateau was observed. However, in the CNS-test series repeated shearing was performed with increasing normal stiffness values K ranging from 0.25 to 16 MPa/mm (see inset in Fig. 4). Dilation in the shear plane triggers an increase in normal stress counteracting dilation with K being the proportional constant causing shear stresses to increase as well.

Constant normal load (CNL) test
To test the newly developed FFS approach under CNL boundary conditions the laboratory shear experiment on   a granite sample with identical boundary conditions as presented in Fig. 4a was duplicated with the FFS model approach. The granite sample was consecutively sheared four times with increasing level of normal stress in each step. The geometry of the shear plane before loading was determined by means of a 3D-Scanner using a grid constant of 0.5 mm. The quadratic grid interpolated from the 3D-coordinates gained by the scanning device serves as initial shear surface topography. As the FFS simulation approach solely requires fundamental rock parameters which commonly are determined by standardized testing procedures, the values given in Table 2 can directly be used as model input defining material properties. The initial normal stress on the shear plane σ n and the normal stiffness of the loading system K are assigned as loading boundary conditions. For the case of a CNL test K = 0 is chosen providing purely static boundary conditions. After each shearing step the sample was moved back to its original position, normal stress was increased to next load level and shear displacement started. Therefore, in the simulation as well as during the lab tests, the resulting geometry at the end of one loading step was used as initial geometry in the following step. The results of the first and last of the four loading steps are illustrated in Fig. 5.
In the first shearing step with a normal stress of σ n = 1 MPa the peak shear strength of the laboratory experiment is not reached by the simulation. Neither the new code nor the classical approach by Barton-Bandis does predict the pronounced peak strength followed by a significant stress drop. This deficiency might be attributed to lacking efficiency of the workflow of geometrical characterization of the shear plane showing a natural roughness and true 3D-features. The 3D-scanning device has technical limitations in terms of resolution and accuracy. By extensive calibration of the device and superposition of several scans from various directions accuracy of the 3D-data points was enhanced up to 50 µm and spatial resolution up to 0.25 mm. However, as the resolution in the new code is limited at present to a grid constant of a = 0.5 mm, information on features on a smaller scale is lost and not available in the calculation model. Effects of micro-roughness on strength development and on the peak strength can, therefore, not be accounted for in the current simulations. Shear stress values after increased shear displacement at the end of the shearing step is fairly well predicted with newly proposed calculation model.
In the fourth shearing step with a normal stress of σ n = 7.5 MPa both approaches-lab experiment as well as model calculation-show a steep increase of shear stress at initial shear displacement up to approximately 0.5 mm followed by a plateau at residual strength level. It is obvious that in this case the simulation over-predicts the shear stress. Clearly the proposed calculation model can well predict residual shear strength values for initial loading, but over-estimates the final strength for consecutive loading steps. Underestimation of damage of the shear plane, such as abrasion, asperity damage and local breakage, in the new calculation approach can explain this nonconforming behaviour. Thus, the shear plane at the start of a subsequent loading step in the calculation is more intact than in nature yielding higher shear resistance for the model compared with the experiment.
To sum up for the CNL-case, the quality of the calculation results using the FFS approach is fair considering the fact that no fitting parameter was used. In particular the general trend from an initially brittle failure mode showing a clear peak followed by a significant stress drop to an almost perfectly plastic behavior with a residual strength plateau is reproduced. Concerning development of stress with displacement, i.e., shear displacement to activate shear strength, the new model is an improvement with respect to the classical Barton-Bandis approach.

Constant normal stiffness (CNS) test
CNS boundary conditions have been implemented to backcalculate a lab shear test on a basalt sample (Fig. 4b). Sample dimensions, topography of the shear plane as well as test boundary conditions of the lab experiment have been duplicated in the model calculation. Main part of the test sequence is shearing the sample in seven consecutive steps with increasing normal stiffness K step by step (Fig. 3). After reaching final shear displacement in one step, normal load on shear plane was removed completely and the sample was reset into its original position before shearing starting again. A normal load of σ n = 1 MPa was used as a starting point for each step. Damage of the shear plane accumulated from step to step and was characterized my means of 3D-scanning before and after the experiment. Dimensions of the basalt sample exceeded dimensions of the granite sample used in CNL calculation. Thus, grid constant was chosen as a = 1.0 mm to limit calculation time to 1 h. For the CNS model normal displacement during shearing, i.e., dilation, is a key feature for mechanical behavior as it is directly linked to normal stress via normal stiffness K. For a first model assumption the dilatation angle was set equal to the maximum inclination angle in shear direction after correction for local asperity damage [Eq. (9)] and elastic deformation [Eq. (10)]. This is a pure geometrical assumption always keeping in mind the main objective to keep the model physical motivated, as simple as possible and definitely without fitting parameters. This model approach assumes far field constraints forcing the sample halves to move as rigid blocks and parallel to each other. This is the realistic case when the sample forms a local feature in a larger sized shear zone with its halves being fixed in the hanging and lying wall, respectively.
In Fig. 6, shear displacement-shear stress curves for normal stiffness values of 0.25 MPa/mm and 2.0 MPa/mm are presented. For the low stiffness value, the simulation is in good agreement with the experimental results, whereas for higher stiffness the shear stress is overpredicted by the new code. This problem is caused by an overestimation of dilation by the new code as shown in Fig. 7. The normal displacements calculated on basis of the 3D-surface data differ substantially from the data measured in the lab experiment. Dilation data in the new model is governed by the maximum surface gradient in shear direction. In the simplified model approach shear stresses are directly linked to normal stresses which are directly proportional to normal displacement. Thus, local variations, unrealistic or erroneous shear plane topography will significantly violate shear stress development. Clearly this circumstance will have its main impact for shear tests comprising several steps with damage accumulation during the test. The curves shown in Figs. 6b and 7b are the 4th step in the series. Any abnormal or unrealistic damage calculated in any previous steps will, therefore, unrealistically alter the shear plane behavior in a later stage of the test.
A review of the detailed picture documentation of the sample between consecutive steps in the lab shear test depicts that the upper left corner of the shear plane is heavily fractured, i.e., completely destroyed and eventually broken out and even missing in the physical experiment. However, the newly developed approach is not able to model such extreme behavior. It does consider asperity damage at a local scale (Eq. (9)) but does not predict deep-seated failure as experimentally observed for the sample examined here. To account for this deficiency, the upper left corner was manually removed from the shear plane model. The calculation was started again and the corresponding shear displacement-shear stress curves are given in Fig. 8. Shear stress development for a low normal stiffness value of 0.25 MPa/mm is now slightly underestimated with the new model. However, for a high normal stiffness value of 2.0 MPa/mm major improvement of the calculation results giving a more realistic result is noticed. Without doubt, the calculated shear stress curve shows high deviations especially after small shear displacements of 1 to 5 mm. However, considering the simple model approach on a strict physical basis not using any kind of fitting parameters the results are fairly realistic.

Damage evolution on the shear plane
A view of the damaged shear plane in basalt at the end of the CNS experiment is presented in Fig. 9. Direction of shear movement was towards left hand. The red color overlay marks zones which are marked as damaged by the new code. Original color of the shear plane is dark gray/black, damaged areas appear as light gray after the experiment. Care has to be taken, because broken and grinded rock particles, i.e., gouge material, is displaced due to shear movement and will color adjacent regions with grayish fines even if there is no physical damage observed. It is clearly visible that damage of the shear plane is a very localized phenomenon depending  (Barton et al. 1985) are presented alongside with the new model results For CNS boundary conditions the damage pattern of the shear plane is determined by the systems normal stiffness value K. The result of a parameter study with four different stiffness levels (2, 4, 8 and 16 MPa/mm) are given in Fig. 10. It is obvious, that size and distribution of damaged areas is increasing with increasing normal stiffness. Damage starts to evolve at dominant, primary features. For the sample presented in Fig. 10 these are the upper left corner and a curved structure starting from mid of the upper edge proceeding to the lower left corner as isolated domains. With the increase of normal stiffness of the surrounding rock mass, these isolated domains increase in size and eventually connect to bigger ones. Starting from K = 16 MPa/mm a considerable number of new areas is damaged, which has not been addressed before. Due to higher normal stresses in this case secondary, less distinct features establish contact and are damaged thereafter. This behavior has major impact on subsequent hydromechanical coupled phenomena of fractured rock (Dang et al. 2019). As the new calculation scheme can realistically duplicate this behavior based on a physical motivated and fast modelling approach it can not only be used for mechanical stress-displacement calculations but also as pre-processor generating a set of physically sound models for subsequent flow and transport simulations.  (Barton et al. 1985) are presented alongside with the new model results

Conclusions
To simulate shear behavior of rock joints in crystalline rock a MATLAB-code using an FFS approach was developed. Its main objective is to simulate shear behavior by just using physical-founded input parameters and realizing minimal calculation times.
The newly developed code was tested against results of lab experiments gained in series of direct shear tests on crystalline rock samples using different boundary conditions. The quality of agreement between lab test and simulation results is fair. The code can realistically predict shear stress development and residual strength values for a wide range of boundary conditions and load levels. Deviations between measured and calculated results appear whenever geometric properties of the sample are either not exactly defined due to the test-specific procedure (e.g., damage accumulation in subsequent test stages) or show a high amount of local bias (e.g., non-uniform surface gradients in parts of the shear plane). As the new model approach strictly calculates shear behavior on single grid points and does not use any spatial averaging or smearing of features, one single erroneous node will affect the behavior of the entire model. This problem is especially apparent for CNS boundary conditions, where local plane topography will directly feedback on shear stress development. At the current stage of code development manual adaption of shear plane's geometry data in accordance with observations in lab experiments is necessary and significantly improves the quality of the results. Identification and localization of damaged zones and areas performs very well and delivers physically sound results which can be the starting point for fluid and transport modelling on the shear planes. In comparison to the classical function set of Barton et al. (1985) the newly developed code delivers results of similar quality with less computer power needed.
One main advantage of the FFS approach is that only a few standardized (and, therefore, readily available) parameters for the rock matrix and the joint surface are needed for evaluating the coupled shear stress-displacement-dilation performance of a rock joint. A second point is that the calculation time is very short especially when compared to numerical calculation schemes. Therefore, the newly proposed simulation approach could be applied in the early stages of design of geotechnical constructions when (a) typically detailed rock mechanical parameters are not available, yet and (b) numerous alternative versions need to be considered and calculated to determine one or several best-fit solutions.
Further development steps of the code will comprise the consideration of two joint planes with different topography, i.e., a non-matching joint, and enhancements with respect to damage, fracturing and plastic deformation of the rock. A more realistic behavior will be achieved if deep-seated damage and plastic stress-redistribution will be incorporated in addition to local shear damage of asperities.

Open data
The full data sets of the CNL and the CNS test, the code files and a documentation of the code can be downloaded.

Conflict of interest
The funding of the GeomInt project by the Federal Ministry for Education and Research (BMBF) under Grant 03G0866A is highly acknowledged.
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/.