Numerical Simulation of Shear Behavior and Permeability Evolution of Rock Joints with Variable Roughness and Infilling Thickness

State key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China Key Laboratory of Deep Coal Resource Mining (CUMT), Ministry of Education of China, School of Mines, China University of Mining and Technology, Xuzhou, Jiangsu 221116, China Department of Energy and Mineral Engineering, EMS Energy Institute, and G3 Center, The Pennsylvania State University, University Park, PA 16802, USA


Introduction
The mechanical and hydraulic properties of rock joints are of great significance in rock engineering, including near-field modeling of activated faults [1][2][3], petroleum, shale, and geothermal reservoirs [4][5][6][7][8][9].Shear stimulation of existing fracture networks is considered a method for permeability enhancement in geothermal systems [10][11][12].Similar to fault zones (Figures 1(a) and 1(b)) with continuous infilled gouge layers [13,14], weathered joints (Figure 1(c)) are generally filled with undeformed minerals, which can significantly influence the mechanical properties and permeability of rock masses [15].The shear behavior and permeability of filled discontinuities in rock are determined by several parameters, including the infill thickness and joint surface roughness [16][17][18].Although the influence of joint roughness and infill thickness on rock mechanics has been studied previously, the coupled effect of these parameters on the mechanical behavior and permeability evolution of filled discontinuities is rarely reported.
Shear failure occurs when a frictional threshold value is reached.A sand-infilled rock joint can dramatically reduce the shear failure strength of a rock mass, and the infilled materials can dominate the shear behavior of intact rock due to its relatively low frictional properties [19][20][21][22].The infill thickness plays a vital role in estimating the mechanical behavior of a rock mass.Specifically, a thicker infill implies a smaller peak shear strength [23,24].Several models have been proposed to predict the peak shear strength of infilled joints by considering the ratio of infill thickness (t) to the height of an idealized (i.e., regular or planar) joint wall asperity (a) [25][26][27][28][29].In these models, shear strength decreases with t prior to reaching the critical value t/a.However, natural rock joints commonly have undulating surfaces, and the resistance of the roughness may pose a significant difference to that of a hypothetical smooth surface.Hence, it is difficult to estimate the shear strength of natural joints because a reasonable t/a ratio is not easily determined.
In order to characterize natural joint roughness, Barton [30] introduced the joint roughness coefficient (JRC), which varies from 0 to 20, corresponding to the smoothest to the roughest joint surfaces.Based on the JRC, Barton and Choubey [31] established an empirical law of rock joint friction to estimate the shear strength of rock joints.Numerous subsequent studies have been conducted to investigate shear strength criteria [32][33][34][35], mechanical behavior [36][37][38][39], and porosity and/or permeability evolution of infilled joints [40][41][42][43].Coupled shear-flow experiments of a single rock joint have also been performed to study the effect of joint shear deformation on the permeability evolution of rock joints [44][45][46][47][48][49][50][51][52].Specifically, Ye et al. conducted a series of novel shear stimulation experiments to probe the fundamental principle of permeability enhancement in rock joints.The experimental results demonstrate that dilatant shear deformation and cracks propagation are two crucial and integral mechanisms for permeability creation in rock joints.
Fluid flow through a single rough fracture varies with shear displacement.More specifically, the permeability of a rock mass decreases with loading stress and increases with shear-induced dilation.However, research on the coupled effects of JRCs and infill thickness on permeability evolution is underreported in the literature.In this study, we numerically simulate the shearing process of infilled rock joints using particle flow code (PFC) software.The paper is organized as follows: the distinct element numerical method and material parameters are introduced in Section 2, the simulation procedure is illustrated in Section 3, and results and discussion are presented in Section 4. Finally, several important conclusions are drawn in Section 5.

Numerical Method and Material Parameters
2.1.Description of DEM and PFC2D.The distinct element method (DEM) is based on solving Newton's second law of motion for an assembly of particles with a predefined constitutive model applied to the grain-grain contacts [53].Particle flow code 2D (PFC2D) is a two dimensional program based on DEM theory [54,55].In PFC2D, particles are idealized as rigid particles that interact with one another at their contacts.
The calculation cycle in PFC2D is a time-stepping algorithm, and the contact force is updated by the forcedisplacement law.When particles come into contact, the contact forces are calculated as a function of relative displacements and specified stiffness.The linear relationships can be expressed as follows: where F n i represents the total normal force, K n is the normal stiffness, U n is the total normal displacement, n i is the unit normal vector to the contact plane, ΔF i s is the relative shear force, and K s and ΔU i s are the shear stiffness and relative shear displacement, respectively.
For sand-infilled discontinuities, we only consider the frictional parts of particles.Thus, in the direct shear numerical simulation, shear slip occurs when the iteratively checked contact force exceeds the maximum allowable shear contact force.
where μ represents the friction coefficient at the contact.  2 Geofluids mechanical behavior of a cement-like material substance between two contacting particles [56].The parallel bond component acts in parallel with the linear component and establishes an elastic interaction between the pieces (Figure 2(c)).Parallel bonds, which can transmit both force and moment between pieces, have been successfully employed in estimating the mechanical behavior of rocks [57][58][59].

Permeability Enhancement Model.
A sigmoidal logarithmic function is employed to illustrate the permeability versus shear displacement curve [60].
where k max is the maximum increment in permeability, k is the permeability change, d is the shear displacement, n is a constant related to the intrinsic property of the materials, and d 5 and d 95 are shear displacements at which 5% and 95% of the total permeability enhancement has occurred, respectively.According to equation ( 3), the relative difference in permeability is Combining equations ( 3) and ( 4), we obtain Equation ( 4) can be rearranged as where k 0 is the initial permeability, ε is the shear strain, ε 5 and ε 95 are shear strains at which 5% and 95% of the total permeability has occurred, respectively, and n equals 11.0.

Material Parameters.
As shown in Figure 2(a), the numerically simulated granite specimen consists of upper and the lower halves and the infilled material (assumed as nonbonded to which a zero bond strength is assigned).The elastic behavior of the granular assembly is represented by assigning an effective modulus and normal-to-shear stiffness ratio to the grain-to-grain contacts.The average, maximum, and minimum ball radius were 3 9 × 10 −3 mm, 3 0 × 10 −3 mm, and 4 8 × 10 −3 mm, respectively.Each ball had at least three contacts with other neighboring balls, and balls without contacts were algorithmically identified as floaters.Different material parameters, including the normal-toshear stiffness ratio, calibrated Young's modulus, and uniaxial compressive strength (UCS) are listed in Table 1.

Numerical Simulation Design
The direct shear configuration consists of two sets of densely packed particles of nonuniform size and a gouge sandwiched between them (Figure 3(a)).The two sets represent the upper and lower rocks, respectively.In the DEM modelling, 4811 particles were packed into a rectangular box with a length of 60 mm and height of 40 mm.The profiles of the infilled assembly were attributed by fitting five typical joint roughness curves [31].In the numerical simulation of the shear tests, five groups of specimens (i.e., T1 to T5) were 3 Geofluids compacted, and five typical JRCs (i.e., 0-2, 4-6, 8-10, 14-16, and 18-20) and five infill thickness ratios ranging from 0.02 to 0.2 were considered (Figure 3(b)).The simulation parameters are listed in Table 2. Figure 4 shows several assembled specimens with variable roughness and infill thickness.
The simulated specimens were compacted and tested as follows: (1) generation of rock specimens, (2) identification of the roughed infill with a certain thickness and JRC, (3) removal of floating particles, (4) application of a linear parallel contact bond to the intact rock and fill, (5) application of a constant normal stress of 20 MPa on the wall using the servo control algorithm in all numerical runs, (6) shearing of the upper half using a shear rate of 0.1 mm/s, and (7) termination of the shearing process when a shear strain of 0.10 is reached.During the shearing process, mechanical and physical parameters (e.g., shear stress and porosity) were recorded in real time.

Results and Discussion
Our numerical simulation results indicate that evolution of the shear strength and permeability of rock joints is strongly related to the infill thickness ratio (t/h) and JRC.Section 4.1 describes the effect of the thickness ratio on the shear behavior of infilled joints, and Section 4.2 demonstrates   As can been seen in Figures 5(a-2)-5(e-2), a greater infill thickness ratio implies a higher peak shear strain.Similarly, the residual shear strength increases with the infill thickness ratio.However, the residual shear strengths are similar after a relative high thickness ratio is exceeded (e.g., t/h = 0 10).The shear strength is therefore predominantly governed by the infill and trends toward a stable value [26].
Peak and residual shear strengths are plotted in Figures 6(a) and 6(b), respectively, and the simulated data are listed in Table 3.A higher thickness ratio clearly implies a smaller peak shear strength, which is consistent with previous studies [22,23,61].A critical thickness of 0.05 divides the peak shear strength versus thickness ratio curves into two stages.Prior to the critical value, the peak shear strength drops sharply with increasing thickness ratio, while a slight decrease can be identified afterwards.Interestingly, a higher residual shear strength, corresponding to a greater thickness ratio, is observed before reaching a critical thickness ratio of 0.10.5 Geofluids specifically, a thicker infill implies a smaller shear strength [23,24].The relationship between shear strength and thickness ratio can be expressed by the following function [26]: where τ/σ n is the normalized shear strength of the infilled or unfilled rock joints, △τ is the variation of the stress (MPa), σ n is the normal load (MPa), t is the thickness of the infill (m), and a refers to the height of the rock (m) and is a constant.α and β are constants that depend on the normal load and surface roughness, and their magnitudes can be obtained by a regression of the simulated data.
According to the hyperbolic function, Indraratna et al. [24] developed a conceptual normalized shear strength model for infilled joints that can be expressed as where φ p is the friction angle of intact rock, k cr is the ratio of infill thickness to the critical thickness (k cr = t/t cr ), φ f is the friction angle of the infilling material, and m, a, and b are constants.Curves fitted according to equation (8) and key regressive parameters are presented in Figure 6(a).The mechanical parameters are listed in Table 3.

Coupled Effects of t/h and JRC on the Permeability
Evolution of Infilled Joints 4.2.1.Permeability Evolution.Dilation or compaction usually occur during the direct shear test [31,62] and can have a substantial effect on the permeability evolution of an infilled rock mass, which can be directly measured in physical experiments by monitoring change in specimen layer JRC: 14-16, t/h = 0.02 Fitted curve t/h = 0.02 JRC: 14-16, t/h = 0.05 Fitted curve t/h = 0.05 JRC: 14-16, t/h = 0.10 Fitted curve t/h = 0.10 JRC: 14  Assuming that the rock matrix is assumed impermeable and fluid only flows through the porous infills within the void spaces of the fracture.Permeability is normally enhanced during shear slip and dilation.To illustrate the relationship between permeability and shear strain, we employed an existing model to predict the permeability of infilled joints by iterating the derived porosity.The model proposed by Segall and Rice [63] using the test data from Marone et al. [46] indicates that permeability evolution can be estimated by either porosity or specimens thickness [64,65].An estimation of the permeability can be expressed as follows: where k/k 0 is the relative difference in permeability, △H/H is the relative difference in specimen thickness, and △ϕ is the relative difference in porosity.
According to equations ( 9) and ( 10), the relative difference in permeability (for convenience, termed as permeability) can be obtained.The evolution of permeability is illustrated in Figure 7.It is clear that permeability increases approximately linearly with increasing shear strain before the yield stage is reached.However, permeability tends to maintain a stable value after the critical shear strain is attained.The thickness ratio has a key effect on the permeability of infilled rock joints.More specifically, permeability increases with the thickness ratio.The simulated data are included in Table 4.
As shown in Figure 7, permeability increases rapidly with shear strain before the yield stage is reached, during which the permeability trends to a constant value.The relation between the relative difference in permeability and shear displacement can be represented by a sigmoidal logarithmic function as shown in equation (6).The calibrated curves are plotted in Figure 7.The regressive curves show good agreement with the simulated results.In addition, the increasing stage, yield stage, and stable stage are also well illustrated.

Effects of t/h and JRC on Permeability Evolution.
Figure 8 and Table 4 show the effects of the JRC and thickness ratio on the relative difference in permeability.It is clear that permeability increases with both t/h and JRC.In the permeability versus thickness ratio curves, permeability increases rapidly at first; however, the growth ratio then decreases gradually after the peak permeability is reached.This result is consistent with the findings by Marone et al. [46].
Permeability increases with both the thickness ratio and JRC, which may be explained by the following two reasons.Firstly, loose media in joints can be assumed as tectonically crushed rock materials or the product of the decomposition or joint weathering, and the bond strength between joint surfaces is zero [23,24].The particles of infilled materials are unbonded but can still resist one another through surface friction, which means that these particles would be rearranged throughout the entire test.This rearrangement is affected by the applied normal stress and joint roughness.The rearrangement activity would be more intense with a thicker infill, which may result in a larger porosity variation as particles mechanically degrade one another.Secondly, as can be seen from Figure 8(b), a higher JRC corresponding to greater permeability due to greater joint dilation can be expected with a relatively higher JRC.

Conclusions
In this paper, the coupled effects of JRC and sand-infilled thickness on the shear behavior and permeability evolution were studied by conducting numerical simulations using PFC2D.Five groups of JRCs were adopted to represent joints with different roughness, and five groups of infill thickness ratios were used to study the impact of infill thickness on the shear behavior and permeability evolution.Several important conclusions can be drawn.
(1) Peak and residual shear strength were determined in the simulated direct shear tests, and the effects of the infill thickness are significant.The peak shear strength decreases with the thickness ratio in the form of a hyperbolic function.A critical thickness ratio of 0.05 divides the thickness ratio-peak shear strength curves into two stages, i.e., sharp drop stage and slight decrease stage.However, the residual shear strength increases with the thickness ratio before a critical thickness ratio is reached, after which point the residual shear strength trends toward a constant value (2) The permeability evolution during a direct shear test can be estimated using the relative porosity difference.Permeability is enhanced by increasing shear strain, and a sigmoidal logarithmic function can be employed to illustrate their relationship

Figure 1 :
Figure 1: (a) 3D schematic diagram of the infill in a rock mass during shearing.(b) Photograph of a fault with internal alignment of infilled materials (modified from Kirkland et al. [66]).(c) Photograph of a joint filled with clay minerals.

Figure 2 :
Figure 2: Infilled joint and contact models.(a) Numerically simulated specimen with an infilled joint.The yellow and red particles represent the intact rock and infill, respectively.(b) Magnified picture of Figure 2(a).(c) Components of the linear parallel contact bond model.

Figure 3 :
Figure 3: Simulation design and distribution of infill profiles with different joint rough coefficients (JRCs).(a) The height and width of the simulated specimens are 40 mm and 60 mm, respectively.The upper half of the specimen moved to the right at a constant shear rate, while the lower part was confined by the servo-controlled wall to maintain the constant normal stress.(b) Profiles of the employed JRCs.

Radius 3 × 10 − 3 ~4 8
of the thickness ratio and JRC on the permeability evolution.4.1.Effect of t/h on the Shear Behavior of Infilled Joints 4.1.1.Peak and Residual Shear Strength.Figure5illustrates shear stress-shear strain curves of infilled specimens with various thickness ratios and JRCs.Shear behavior near the peak shear strength is shown in the magnified plots, i.e., Figures5(a-2)-5(e-2). Prior to the peak shear strength, elastic and yield stages increase quickly, and the shear stress-shear strain curves undergo a gradual decline after reaching the peak strength.

4. 1 . 2 .
Peak Shear Strength Model.Infilled material in joints can reduce the peak shear strength of rock masses.More

Figure 6 :
Figure 6: Relation between the infill thickness ratio and normalized shear and residual strengths.(a) Normalized peak shear strength versus thickness ratio.(b) Normalized residual shear strength versus thickness ratio.

Figure 8 :
Figure 8: Effects of the (a) thickness ratio (t/h) and (b) joint roughness coefficient (JRC) on permeability evolution.

Table 1 :
Mechanical and physical parameters of granite and infill.

Table 3 :
Normalized peak and residual shear strengths.

Table 4 :
Effect of thickness ratio (t/h) and the joint roughness coefficient (JRC) on yield permeability.The thickness ratio and JRC have coupled effects on the permeability evolution of infilled rock joints.More specifically, permeability increases with both the thickness ratio (t/h) and JRC.The permeability versus thickness ratio curves show that permeability increases significantly at first; however, the growth ratio decreases gradually when the thickness ratio reaches a critical value