Influence of fracture roughness on shear strength, slip stability and permeability: A mechanistic analysis by three-dimensional digital rock modeling

Chaoyi Wang, Derek Elsworth, Yi Fang, Fengshou Zhang

a Department of Physics and Astronomy, Purdue University, West Lafayette, IN, 47907, USA

b Department of Energy and Mineral Engineering, The Pennsylvania State University, University Park, PA,16802, USA

c Department of Geosciences, The Pennsylvania State University, University Park, PA,16802, USA

d Institute for Geophysics, Jackson School of Geosciences, The University of Texas at Austin, Austin, TX, 78712, USA

e Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai, 200092, China

f Department of Geotechnical Engineering, Tongji University, Shanghai, 200092, China

ABSTRACT Subsurface fluid injections can disturb the effective stress regime by elevating pore pressure and potentially reactivate faults and fractures. Laboratory studies indicate that fracture rheology and permeability in such reactivation events are linked to the roughness of the fracture surfaces.In this study,we construct numerical models using discrete element method (DEM) to explore the influence of fracture surface roughness on the shear strength, slip stability, and permeability evolution during such slip events. For each simulation, a pair of analog rock coupons (three-dimensional bonded quartz particle analogs) representing a mated fracture is sheared under a velocity-stepping scheme. The roughness of the fracture is defined in terms of asperity height and asperity wavelength.Results show that(1)Samples with larger asperity heights (rougher), when sheared, exhibit a higher peak strength which quickly devolves to a residual strength after reaching a threshold shear displacement;(2)These rougher samples also exhibit greater slip stability due to a high degree of asperity wear and resultant production of wear products;(3)Long-term suppression of permeability is observed with rougher fractures, possibly due to the removal of asperities and redistribution of wear products, which locally reduces porosity in the dilating fracture; and (4) Increasing shear-parallel asperity wavelength reduces magnitudes of stress drops after peak strength and enhances fracture permeability, while increasing shear-perpendicular asperity wavelength results in sequential stress drops and a delay in permeability enhancement. This study provides insights into understanding of the mechanisms of frictional and rheological evolution of rough fractures anticipated during reactivation events.

Keywords:Fracture reactivation Fracture permeability evolution Fracture roughness Roughness anisotropy Slip stability

1. Introduction

Human interventions into the subsurface, such as hydraulic fracturing, enhanced geothermal stimulation, and carbon sequestration, involve injecting large volumes of fluid at high overpressures. Such interventions disturb the stress field by elevating pore pressure and altering far-field stress (Elsworth et al., 2016),potentially resulting in reactivation and seismic rupture of preexisting faults and fractures. Hydraulic fracturing, in particular,attempts to create engineered fracture networks by fluid injection to stimulate hydrocarbon production. These hydraulic fractures,while creating the possibility to extract hydrocarbon resources from tight shale, can be extremely vulnerable to seismic failure upon stress perturbation (Zoback and Gorelick, 2012; Ellsworth,2013; Walsh and Zoback, 2015), causing hazardous consequences.One key question in understanding the seismic cycle is in unraveling the evolution of shear strength,stability,and permeability of faults and fractures that may contribute to dynamic slip events.

Previous laboratory shear experiments on faults and fractures showed that permeability declines over shear slip, which may be caused by clay swelling or clogging of wear products (Fang et al.,2017a, b; Im et al., 2017; Ishibashi et al., 2018). These observations reflect the friction-stability-permeability relationships at low confining stresses for fractures with small roughness, i.e. asperity size of the order of micrometers.A recent study reports that in the injection-induced shear slip experiment,asperity degrades but the surface roughness would still contribute to increasing permeability(Ye and Ghassemi, 2018). Investigations on fabricated fractures with controlled roughness suggest that roughness patterns exert a strong control on permeability evolution via competitive effects of compaction and dilation during shearing(Zhang et al.,2017,2019;Fang et al., 2018). However, due to the lack of direct tracking in laboratory conditions, it is still not clear how mechanistically the asperities of fractures evolve during shear, which further controls the permeability evolution.

Classic experimental studies have suggested that the shear strength of fractures are closely linked to surface roughness(Barton,1973; Barton and Choubey,1977). Empirical indices have been developed to describe the roughness of rock surfaces in the order of millimeters using parameters such as joint roughness coefficient (JRC) and joint compressive strength (JCS). Correspondingly, rougher surfaces undergo greater dilation when sheared, resulting in an increase in aperture and enhancement of permeability. However, breakage and degradation of asperities can lead to impeded dilation, reduced aperture, and reduced permeability.Laboratory shear experiments have been conducted on rough samples which feature“saw-tooth”or sinusoidal shaped asperities (Asadi et al., 2013), and the results suggest that the breakage and degradation of asperities are linked to normal stress,bonding strength,and asperity geometry.These studies,however,lack proper reproduction of the stochastic characteristics of natural rough rock surfaces. Mathematical algorithms have been developed to describe the natural roughness of rock surfaces(Brown and Scholz, 1985). The theories suggest three key parameters to describe a rough rock surface: (1) The root mean square (RMS) roughness including variance of amplitude and distribution of the asperities; (2) The length scale for degree of mismatch; and (3) Fractal parameters (Brown,1995).Moreover, a recent study provides an optimized method in approaching roughness of natural rock joints using Fourier series (Yong et al.,2018).

Numerical approaches have been adopted to investigate the effects of roughness on the shear strength of rock joints. While continuum numerical models are implemented to simulate the onset of shear failure of rough joints and fractures, the discrete element method (DEM) (Cundall and Strack, 1979) provides numerical analogs to follow the progression of failure of rough fractures. In DEM models, the surface profiles can be reproduced, and the damage of asperities can be tracked during a simulated shear test (Cundall, 2000). DEM studies have examined the evolution of the shear strength of fractures described by JRC profiles and shown that tensile failure dominates the breakage of asperities and development of micro-cracks (Park and Song, 2009). DEM models have also been used to investigate the evolution of shear strength,slip stability, and permeability of gouge materials (Morgan,1999;Morgan and Boettcher, 1999; Guo and Morgan, 2004; Sun et al.,2016; Wang et al., 2017, 2019). However, no DEM studies to date have explored the ensemble of rough surface profiles featuring stochastic characteristics, and the linkage of shear strength, slip stability, and permeability evolution.

In this study, we report DEM simulation results regarding the influences of surface roughness on the evolution of shear strength,slip stability, and permeability of fractures by utilizing a series of rough surface profiles with variations in asperity roughness,wavelength, and degree of wavelength anisotropy.

2. Numerical method

We use a simplified stochastic algorithm to create rough fracture surfaces and utilize the DEM to construct our numerical model. A modified slip-weakening friction constitutive model is implemented on particle—particle contacts to represent slip evolution and the evolution of fracture porosity and permeability.

2.1. Fracture roughness

Three key parameters may be used to define fracture roughness:(1) The RMS roughness; (2) The length scale for degree of mismatch; and (3) A fractal parameter/dimension. Each of these parameters plays a role in influencing the dynamic response of fractures during reactivation events. In this study,we focus on the first-order effect of roughness on the shear strength, slip stability,and permeability of mated fractures prone to reactivation. Therefore, we simplify the characterization of the surface roughness by considering only the statistical size and distribution of asperities.We use the two key parameters: (1) The RMS height (Sq) of the asperities; and (2) The wavelength (λ) describing the distance between two statistically independent points.The RMS height of the asperities in a sample surface of areaAcan be expressed as

wherezis the individual asperity height and(x，y)is the location of the asperity. In this study, we independently vary the RMS height(Sq), and the wavelength (λ) in the two orthogonal directions (xandy-direction) within the mean fracture plane to characterize different roughness profiles (see Section 2.3).

2.2. Model construction

The model (Fig.1) in this study is developed via Particle Flow Code 3D (PFC3D) utilizing the principles of DEM (Cundall and Strack, 1979). The applicability of DEM to simulate the dynamic response of rocks and faults is summarized elsewhere(Antonellini and Pollard,1995; Morgan,1999, 2004; Abe et al., 2002; Burbidge and Braun, 2002; Guo and Morgan, 2004; Sun et al., 2016; Gao et al., 2018; Wang et al., 2017). The majority of DEM simulations are in two-dimensional (2D) configurations to reduce computational cost yet still produce representative results. However, the negligence of out-of-plane interactions causes the inability of such approaches to reproduce the interaction of fracture planes with rough textures. In this study, we use a three-dimensional (3D)configuration to reproduce fracture surfaces with predefined roughness features.

Virtual rock samples (10 cm × 1 cm × ～0.5 cm(length × width × thickness) each half) with a predefined rough fracture (mated) are sheared to simulate the dynamic response of intact fractures during reactivation events. Specifically, each rock coupon is generated by filling particles in a virtual container with one side replaced by a predefined rough surface (Fig. 1a). The infilled particles are equilibrated to dissipate the kinetic energy caused by initial infilling, i.e. cycling until the ratio of total unbalanced force to total body force is less than 0.001.Particles are then linked into ensemble “l(fā)ithified” samples by bonding (Fig.1b). The pair of virtual rough rock fracture coupons with mated fracture surfaces are slowly brought together and confined (Fig. 1c). Once the incremental confining stress reaches the prescribed magnitude(10 MPa), the upper coupon is displaced laterally to shear against the static lower coupon at a prescribed shear velocity(1 μm/s).The shear velocity is increased to 10 μm/s after 5 mm of shear displacement and cycled back to 1 μm/s after another 5 mm of shear displacement. The cycles continue until a total shear displacement of 25 mm is reached. We present typical fracture roughness profiles with asperity wavelength anisotropy, together with asperity wavelengths implemented in this study,in Fig.1d and e, respectively.

Fig.1. Model construction.(a)Two numerically generated rough surfaces are imported,ready to be brought together.(b)A virtual box(not shown)with one face replaced with one of the rough surfaces serves as a mold for one analog rock coupon.Particles are generated,equilibrated,and bonded inside the mold.A pair of analog rock coupons is generated at the same time.The rough surfaces are removed upon completion of bonding.The particles located closest to the fracture surfaces are marked in blue and red.(c)The two coupons are confined under a prescribed normal stress(10 MPa)and the upper coupon of the sample is loaded to initiate the shear test.(d)Typical roughness profile of the fracture lowerhalf(d-1),wavelength anisotropy is introduced by increasing wavelength in the shear-parallel direction(d-2),and shear-perpendicular direction(d-3);d-1,d-2,and d-3 correspond to simulations rss6, rss9, and rss12, respectively. (e) Schematic of asperity wavelengths in this study, i.e. 0.5 cm,1 cm, 3 cm, and 5 cm.

Fig.2. Contact model between bonded particles:(a)Schematic of the modified linear parallel bond model. A rotation resistance component is included to restrict any free rolling motion; and (b) Evolution of friction coefficient at contacts upon local shear slip.

2.3. Contact model

The interaction of particles in the DEM model is determined by contact models. We use a modified parallel bond contact model(Fig. 2) to describe the particle interaction in the assembly. The modified parallel bond model consists of three main elastic components:(1)Linear elastic spring in the normal and shear directions of the contact,where frictional sliding is achieved by a slider in the shear direction; (2) Linear elastic bonds in the normal and shear directions with tensile strength and cohesion; and (3) Rotation resistance(Iwashita and Oda,1998;Ai et al.,2011;Jiang et al.,2015).The particles are bonded prior to the application of confining stresses. The breakage of bonds is tracked throughout the simulation. The contact model is illustrated in Fig. 2a.

The force and moment within a bonded contact are

whereFcis the contact force andMcis the contact moment. The contact force is resolved as the sum of linear forces(Fl)and parallel bond forcesThe contact moment is provided by parallel bond momentand the contact will not provide moment resistance after the parallel bond is broken.A detailed discussion of the linear component forces is reported by Potyondy and Cundall(2004).The parallel bond force is resolved into normal and shear components and the parallel bond moment resolved into twisting and bending moments:

The tensile and shear stresses within the parallel bond are defined as

The frictional response of the contact may be accommodated by rate-and-state friction (RSF) law. RSF law (Dieterich, 1978; Ruina,1983) has been developed to describe the evolution of friction during slip of faults and fractures. The constitutive relation of RSF may be described as

where μ is the friction coefficient; μ0is the reference friction coefficient;V,V0, andVlpare the current, reference, and load point velocities of the system, respectively; θ,Dc, andkare the state variable, characteristic slip distance, and system stiffness, respectively; andaandbare the stability parameters.

Although RSF is able to match many of the first- and secondorder features in the evolution of friction, it is computationally intensive when implemented in models at contact-level, such as DEM (Abe et al., 2002). A quasi-rate-and-state friction (quasi-RSF)law that replicates key features of the RSF law(Wang et al.,2017)is used here to reduce the computational cost. The constitutive relation (Fig. 2b) may be represented as

where μp, μref, and μssare the peak, reference, and steady-state friction coefficients on the evolving contact, respectively;VlpandVrefare the current and previous global shear velocities from the last velocity step, respectively; andDaccis the accumulated shear displacement on the contact.

A slip event is initiated if the resultant shear stress exceeding the frictional strength of the contact after the contact bond is broken.The evolution of the contact friction depends on the local accumulative shear displacement and the difference between current and previous global shear velocities. During a local slip event (single contact shear),the friction begins to evolve following the red path in Fig. 2b if current global shear velocity differs from previous global shear velocity.If the contact remains active(the two particles are in contact),the contact will evolve along the red path(Dacc＜Dc)and transfer to the purple path at steady state (Dacc≥Dc). The purple path is shown for velocity weakening,but can also follow the path of velocity strengthening,depending on the magnitude of the stability parametersa-bassigned to the contact. Where the contact becomes inactivebeforereaching steady-state(the two particles areno longer in contact),the friction of either contact will remain as-it-is on the red path (Wang et al., 2017). The material properties and parameters used in this study are enumerated in Table A1 in the Appendix. Material properties for uniaxial compressive strength,tensile strength,and Young’s modulus recovered from uniaxial numerical loading experiments(4 cm in diameter,and 8 cm in height)are reported in Table A2 in the Appendix.

2.4. Experiment matrix

We explore the influence of fracture roughness on the evolution of shear strength, slip stability, and permeability of rough fractures during slip events.Specifically, we compare the impacts of (1) RMS height of the asperities(rss1-rss6);(2)Spatial distribution(asperity wavelength)of asperities(rss6-rss9 forx-direction;rss6,and rss10-rss12 fory-direction); and (3) Tensile strength and cohesion of the wall-rock represented by particle bond strength (rss6, rss13, and rss14).The spectrum of experimental variables is noted in Table A3 in the Appendix.

3. Results and analysis

We perform direct shear simulations on analog rock coupons with predefined roughness profiles.The simulations are conducted in velocity stepping mode with velocities up-and down-stepped between 1 μm/s and 10 μm/s over incremented shear offset of 5 mm. The evolutions of shear strength, slip stability, and permeability are evaluated as functions of asperity height,wavelength,and strength.

3.1. Evolution of shear strength

The virtual rock coupons are sheared to a total relative shear displacementof 25mminfivevelocitysteps.Thelowercoupon isheld in place and the upper coupon translates while restrained to deform parallel to the long axis of the fracture.Confining stress is maintained constant at 10 MPa.Shear stress evolution is monitored by the evolution of the friction coefficient of the assembly. The friction coefficient is defined as the ratio of shear stress to confining stress.Since confining stresses are maintained constant during the test,the shear strength scales as friction coefficient.Fig.3 shows the fracture surface profiles and the evolution of the friction coefficient for tests rss1-rss6.The analog coupons feature RMS asperity heights from 0.005 cm to 0.05 cm. The comparison shows an anticipated trend of increasing peak shear strength with increasing RMS asperity height.Shear stress buildsuntil failurewitha pronounced post-peak stress/strengthdrop.Test rss6 shows the highest peak friction(～0.65),and rss1 exhibits the lowest(～0.36).The magnitude of the stress drop increases with the RMS height as does the shear displacement required to mobilize peak strength.All the samples stabilize at a similar coefficient of residual friction(～0.32)after failure.

3.2. Fracture dilation and permeability evolution

Fracture permeability is controlled by the local contribution to ensemble aperture along fluid channels formed by the interconnected pore network. These effects can be estimated by monitoring the evolution of sample dilatancy(scaled to change in sample thickness during shear)and local porosity along the fracture.In this study,the sample thickness is calculated as the distance between the top and bottom of the analog coupons.Local porosity is measured by averaging five evenly distributed and equally sized sampling windows(spheres)placed along the fracture.The evolutions of sample thickness for tests rss1-rss6 are shown in Fig.4.The samples begin with a total height/thickness of ～7.85 mm and gradually dilate to a peak magnitude maintained as a plateau or slight compaction.Samples with low RMS asperity heights (rss1 and rss2) reach a steady sample layer thickness after ～8 and ～10 mm of shear displacement, respectively. Samples with larger asperities (rss3-rss6) dilate more significantly in terms of observed increases in sample layer thickness and also undergo larger shear displacement to reach peak dilation. For example, it takes ～22 mm of shear displacement for the roughest sample(rss6)to reach the maximum dilation — this peak dilation reaches simultaneously with peak strength. After reaching the peak strength, the samples no longer dilate(rss1-rss3)and in some cases compact slightly(rss4-rss6).

The evolution of fracture permeability is estimated from the local change in porosity sampled along the fracture (Ouyang and Elsworth,1993; Samuelson et al., 2011):

wherek/k0is the change in permeability and Δφ is the change in porosity. The change in porosity is calculated as the difference between the initial porosity and the porosity during the simulation.Estimates of fracture permeability evolution(k/k0)during shear for tests rss1-rss6 are shown in Fig. 5. Permeability of the fracture decreases slightly during the first ～10 mm of shear displacement,due to the dominant initial shear compaction.The effect of dilation exceeds shear compaction after ～10 mm of shear displacement,resulting in a net permeability increase. Permeability shares a common trend of increasing in the less rough fractures(rss1-rss3)after a threshold shear displacement (～8 mm, ～10 mm, and～13 mm, respectively) with the permeability then reaching a plateau and stabilizing.However,for fractures with moderate RMS roughness (rss4), the fracture permeability evolves unstably after the initial onset of permeability growth. Finally, for the roughest fractures (rss5 and rss6), permeability decreases following the attainment of peak permeability enhancement.

3.3. Evolution of slip stability

Fig. 3. (a) Lower fracture surfaces of rss1-rss6 before shear with colored contours illustrating the topography of the surfaces (asperities); and (b) Evolution of shear strength interpreted as friction (τ/σ) for rss1-rss6 (RMS asperity heights ranging from 0.005 cm to 0.05 cm). The shear strength of the samples generally builds up until reaching peak strength,followed by a stress drop post-peak,sometimes comprising several successive stress drops.Samples with rougher fractures exhibit a higher peak shear strength and larger threshold shear displacement corresponding to the peak strength. All samples show similar residual shear strength after failure.

The evolution of slip stability of the simulated rough fractures is important in understanding the characteristics of the potential for induced seismicity.Stability may be indexed via the parametera-b. Positive values suggest aseismic behavior and negative values indicate potential seismic behavior.In this study,thea-bvalues are extracted by fitting the RSF law to the detrended friction evolution from the velocity steps. Fig. 6 shows the summarizeda-bvalues for samples rss1-rss6 (a few extreme outliers that result from the interference of closely occurring stress drops are excluded). Thea-bvalues generally scatter around the neutral(zero)line,showing mostly velocity neutral behavior.Thea-bvalues increase with RMS asperity height(from 0.005 cm to 0.05 cm).

4. Discussion

RMS asperity height is shown to play an important role in controlling the shear strength,slip stability,and permeability of the fracture. We observe that larger RMS asperity heights result in higher peak shear strengths and a larger threshold displacement for failure with a larger stress drop. A higher RMS asperity height also promotes greater shear dilation, resulting in fracture permeability increase. However, extremely rough fractures exhibit net reduction in post-peak permeability.Natural fractures are complex systems with anisotropy in roughness and asperity strength. We provide a brief discussion of relative roles of asperity wavelength anisotropy,asperity strength,and finally a proposed mechanism for permeability evolution of rough fractures.

4.1. Influence of RMS height on peak frictional strength

Fig. 4. Evolution of sample thickness for samples with RMS asperity heights ranging from 0.005 cm to 0.05 cm (rss1-rss6). The samples generally dilate until reaching a plateau where the sample thickness either ceases to increase(rss1-rss3)for small RMS asperity heights or slightly compacts (rss4-rss6) for large RMS asperity heights.

Fig. 5. Evolution of fracture permeability (k/k0) for tests rss1-rss6. Permeability decreases slightly with compaction during early shearing. Permeability increases rapidly upon a threshold shear displacement and continues to increase until reaching a peak,after which plateau (rss1-rss3) is observed. Fracture in rss4 shows unstable permeability evolution after reaching the peak. Fractures with large RMS asperity heights(rss5 and rss6) exhibit permeability reduction after reaching peak values.

Fig. 6. Summarized a-b values for different RMS asperity heights (0.005—0.05 cm).The a-b values scatter around zero, indicating velocity neutral behavior. The a- b values increase with RMS asperity height, implying increasing influences of asperity comminution generated wear products.

Fig. 7. Peak shear strength and corresponding threshold shear displacement versus RMS asperity height (from 0.005 cm to 0.05 cm). These two properties are positively correlated to RMS asperity height up to a threshold RMS height, i.e. 0.04 cm.

We have shown that RMS asperity height is closely related to the peak shear strength.Fig.7 shows the relationship between the peak shear strength, correlated threshold shear displacement, and RMS asperity height (0.005—0.05 cm in this study). The peak shear strength increases near linearly with RMS asperity height. The threshold shear displacement for the peak shear strength also increases with RMS asperity height following a similar trend except for RMS heights larger than 0.04 cm. When RMS height reaches 0.05 cm, the threshold shear displacement only increases slightly.These observed trends suggest that the shear strength and shear displacement required to reach failure are linearly related to the RMS asperity height, up to a limiting asperity height (0.04 cm in this study).Thus,the shear strength of rough fractures is not solely determined by asperity heights and wavelengths, but also by the strength of the asperities.

4.2. Influence of roughness anisotropy

The anisotropy of the roughness within the plane of the fracture plays an important role in determining the mechanical and rheological properties. In this study, the anisotropy of roughness is interpreted by varying asperity wavelength parallel to the shear direction (x-direction, rss7-rss9), and perpendicular to the shear direction (y-direction, rss10-rss12) in the plane of the fracture.Fracture surfaces with various degrees of roughness anisotropy(wavelengths of 0.5 cm,1 cm,3 cm,and 5 cm in one direction)are tested. The fracture surface profiles are shown in Fig. A1 in the Appendix.

The evolution of shear strength, permeability, and stability parameters resulting from these profiles are shown in Fig. 8.Increasing asperity wavelength in the shear parallel direction reduces the peak strength and the amount of the stress drop at failure. Fractures with larger wavelengths parallel to the shear direction (rss9) show permeability enhancement at a smaller threshold shear displacement.Sample rss9 also shows a stabilized permeability at near-peak levels without any tendency to decrease,while the permeabilities of rss6-rss8 (fractures with smaller asperity wavelengths in the shear direction)decrease slightly after reaching peak permeability.In terms of stability parameters,Fig.8b shows a slight trend of increasinga-bvalues with increasing wavelengths in the shear direction.

Increasing asperity wavelength perpendicular to the shear direction (Fig. 8c) reduces the peak shear strength and the magnitude of stress drops at failure.Additionally,the stress drop evolves from a single or several abrupt drops at failure with large magnitude(rss6 and rss10)to a series of smaller stress drops with a cyclic form (rss11 and rss12). It takes longer threshold shear displacement for the samples with larger shear-perpendicular asperity wavelength to reach residual shear strength. Moreover,permeability enhancement is also delayed. In terms of stability parameters, the values are broadly scattered with no obvious relationship between the values ofa-band the asperity wavelength.

Fig.8. The evolution of shear strength,fracture permeability,and stability parameters related to the anisotropy of asperity wavelength,i.e.in the shear direction and perpendicular to the shear direction:(a)Shear strength and permeability evolution of fractures with asperity wavelengths in shear direction(clx,Table A3)of 0.5 cm,1 cm,3 cm,and 5 cm(rss6-rss9); (b) Stability parameters of samples rss6-rss9; (c) Shear strength and permeability evolution of fractures with asperity wavelengths perpendicular to shear direction (cly,Table A3) of 0.5 cm,1 cm, 3 cm, and 5 cm (rss6, rss10-rss12); and (d) Stability parameters of samples rss6, and rss10-rss12.

4.3. Influence of asperity strength

The strength of asperities plays an important role in determining the ensemble mechanical strength of fractures. In this study, we discuss the influence of asperity strength by varying the tensile strength and cohesion for the contacts while keeping other parameters constant. Fig. 9 shows the evolution of shear strength, permeability, and stability of samples featuring parallel bond strengths and cohesion of 20 MPa, 50 MPa, and 500 MPa(rss14, rss4, and rss13, respectively). The resulting evolution of shear strength (Fig. 9) suggests that increasing tensile strength and cohesion results in an increase in the peak shear strength.Specifically, rss13 shows ～25% higher peak shear strength than rss4 and ～50% higher than that of rss14, respectively. This trend can be explained by the rationale that lower bonding strength and cohesion result in weaker asperity strength. Therefore,samples with lower bonding strength and cohesion are subjected to increased localized failure and asperity breakage, producing more wear products during shearing. This effect is shown by the number of broken bonds inside the analog coupons, as shown in Fig. 9a (the lower fracture coupons are shown with 95% transparency). Noticeably, the bond breakage tends to localize on the contacting faces of the samples where stress concentrates due to loading.

In terms of permeability evolution, increasing tensile strength and cohesion causes permeability enhancement to initiate both earlier and to a higher degree, as shown in Fig. 9b. Interestingly,rss14 shows almost no enhancement of permeability throughout the simulation. This may be related to the localization of bond breakage on the two sides of the sample and clogging of wear products between the fracture faces.Shear-induced bond breakage produces a relatively large amount of wear products, which plausibly clogs the fluid passage in the fracture, reducing fracture permeability.

Additionally, as shown in Fig. 9c, the stability parameters are mostly negative and show an even broader range of variation with larger bonding strength and cohesion. However, the stability parameters do not show any significant correlation with increasing tensile strength and cohesion.

4.4. Influence of shear-generated wear products

Fig. 9. (a) Fracture profiles, sample geometries, and bonding breakages after 25 mm of shear displacement for samples rss4, rss13, and rss14; (b) Evolution of shear strength and permeability with bond strengths and cohesions of 20 MPa, 50 MPa, and 500 MPa, respectively; and (c) Stability parameters plotted against bond strengths.

The generation of wear products during shearing are commonly observed(e.g.Bakker et al.,2016;Fang et al.,2017c;Im et al.,2018)and also in this study. The generated wear products can be transported and redistributed by mechanical interactions between the two fracture surfaces,and by fluid,if presented(Candela et al.,2014). These wear products can influence the mechanical and rheological properties of the fracture through a variety of mechanisms. The wear products may both impact shear strength of the fracture by localizing shear and clog pores and major fluid pathways within the fracture, staunching permeability evolution. Over longer timescales, mechanical and fluid interactions and reaction between wear products and asperities may result in geochemical transformations that alter surface properties of the asperities and impact the evolution of rheological and transport properties.In this study, the generation of wear products for samples rss1-rss6 is shown as red highlighted particles in Fig. 10. The corresponding numbers of broken bonds are shown in the lower right plot in Fig.10.

It is typically observed that rough fractures with higher RMS asperity heights produce significantly more wear products(rss6 vs.rss1 in Fig. 10). Samples rss1 and rss2 do not show a significant difference in the amount of generated wear products.Samples rss3 and rss4 show a large increase in the amount of generated wear products,and this increase slows down in rss5 and rss6.Noticeably,the threshold shear displacement for bond breakage (left plot of Fig.10), where the slope of the evolution curve becomes abruptly smaller,corresponds to the major stress drop at failure(Fig.3),and the peak in permeability evolution(Fig.5).This behavior indicates a dominant influence of wear products on the evolving shear strength and permeability of the fracture, as indicated by similar residual shear strengths for samples rss1-rss6.Large initial asperity height results in increased dilation and peak permeability, but generates more wear products, reducing permeability by clogging major fluid conduits. Therefore, permeability evolution is potentially dependent on the competitive influence of asperity height(contributing to dilation) and the volume of generated wear products (contributing to clogging), with the dominant process defining the response. When the clogging effect exceeds that of dilation, the permeability of the fracture may be reduced after reaching its peak,even though the fracture begins with higher RMS asperity heights.This mechanism is suggested by the permeability evolution of samples rss5 and rss6 (Fig. 5) and the amount of generated wear products(Fig.10).It is worth noting that the wear products mainly concentrate on the principal contacting portions of the fractures, suggesting that stress concentration (in the laboratory and otherwise)may result in the clustering of wear products at kinks in fractures and fault asperities.

5. Conclusions

In this study,we investigate the influence of fracture roughness on the shear strength, slip stability, and permeability evolution of fractures by DEM modeling. The rough fracture surfaces are generated based on RMS asperity heights and asperity wavelengths. We have discussed the influence of fracture roughness in the following aspects:RMS asperity heights,anisotropy of asperity wavelength,and strength of asperities.Also,we have analyzed the stability parameters,the relationship between peak shear strength and RMS asperity heights, and proposed mechanisms for the permeability evolution for rough fractures. We summarize the conclusions of this study as follows:

(1) Larger RMS asperity height yields higher peak shear strength while requiring more shear displacement to reach the peak strength.The relationship between asperity height and peak shear strength is positively correlated, but in a nonlinear fashion, i.e. under a given asperity strength, there is a limit for the peak shear strength for rough fractures with increasing RMS asperity heights.

(2) Increasing the RMS asperity height can alter slip stability of rough fractures from mostly velocity weakening to velocity strengthening. This transformation of slip stability can be related to the generation of wear products.

(3) Anisotropy of asperity wavelength can influence the shear strength and permeability evolution of rough fractures.Larger asperity wavelength parallel to the shear direction reduces the peak shear strength of the fracture while increasing the fracture permeability. Larger asperity wavelength perpendicular to shear direction can slightly reduce the peak shear strength, induce more frequent stress drops during failure while delaying, or suppress the permeability enhancement.

(4) The strengths of the asperities (bonding strength and cohesion) are crucial to determine the shear strength and permeability evolution of the rough fractures.Lower asperity strength results in lower shear strength and less permeability enhancement.

(5) The amount, distribution, and transport of shear-generated wear products can dominate the evolution of shear strength,slip stability,and permeability of rough fractures by localization and clogging effects. Fractures with more wear products exhibit lower shear strength, enhanced slip stability, and lower permeability during dynamic shear.

Conclusions drawn above are specifically applicable to the parameters and situations in this study with the potential of upscaling to the field. Future study may consider the deformation and crushing of individual analog particles. Nonetheless, our study provides a straightforward way to study the influence of surface roughness on the mechanical and rheological properties of fractures.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

This work is a partial result of the support provided by United States Department of Energy Grant DE-FE0023354.This support is gratefully acknowledged. This work utilizes data from literature which are cited in the main reference list,and data from numerical modeling of this study are shown in the main text.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2019.12.010.