Three-dimensional double-rough-walled modeling of fluid flow through self-affine shear fractures

Richeng Liu, Ming He, N Hung, Yujing Jing, Liyun Yu

a State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou, 221116, China

b School of Engineering, Nagasaki University, Nagasaki, 852-8521, Japan

c School of Petroleum Engineering, China University of Petroleum, Qingdao, 266580, China

Keywords:Double-rough-walled fracture Navier-Stokes (NS) equations Anisotropy Fractal dimension

A B S T R A C T This study proposes a double-rough-walled fracture model to represent the natural geometries of rough fractures.The rough surface is generated using a modified successive random additions (SRA)algorithm and the aperture distribution during shearing is calculated using a mechanistic model. The shear-flow simulations are performed by directly solving the Navier-Stokes (NS) equations. The results show that the double-rough-walled fracture model can improve the accuracy of fluid flow simulations by approximately 14.99%-19.77%, compared with the commonly used single-rough-walled fracture model.The ratio of flow rate to hydraulic gradient increases by one order of magnitude for fluids in a linear flow regime with increment of shear displacement from 2.2 mm to 2.6 mm. By solving the NS equations, the inertial effect is taken into account and the significant eddies are simulated and numerically visualized,which are not easy to be captured in conventional experiments.The anisotropy of fluid flow in the linear regime during shearing is robustly enhanced as the shearing advances;however,it is either increased or decreased for fluids in the nonlinear flow regime, depending on the geometry of shear-induced void spaces between the two rough walls of the fracture. The present study provides a method to represent the real geometry of fractures during shearing and to simulate fluid flow by directly solving the NS equations, which can be potentially utilized in many applications such as heat and mass transfer,contaminant transport, and coupled hydro-thermo-mechanical processes within rock fractures/fracture networks.

1. Introduction

Hydraulic properties of fractured rock masses play an important role in estimate of the performances of underground projects such as enhanced oil recovery,geothermal energy development,nuclear waste disposal (Pruess et al., 1990; Cvetkovic et al., 2004). To facilitate the modeling of fluid flow,an assumption is made that the fluid flow within rock fractures obeys the cubic law derived from the complex Navier-Stokes (NS) equations by neglecting the nonlinear terms. The flow rate is linearly proportional to the pressure drop(Bear,1972;Tse and Cruden,1979),using the models that have two smooth and parallel walls(Witherspoon et al.,1980).The flow velocity in vertical (normal to the fracture surface) direction is also assumed to be negligible.However,in the vicinity of wells during pumping tests and/or in the karst systems, the flow rate is nonlinearly correlated with the pressure drop and the cubic law neglecting the inertial effects is not suitable(Chen et al.,2015;Xia et al., 2017). The simplified smooth parallel plate model can smooth the streamlines of fluid flow,which thus overestimates the conductivity of a fracture (Yeo et al., 1998; Bauget and Fourar,2008). Therefore, it is critical to establish a rough-walled fracture model to represent the natural geometries of fractures and to simulate fluid flow by directly solving the NS equations without any simplifications.

Although the parallel plate model(constant aperture across the fracture plane) is commonly used in large-scale practical reservoir simulation studies(Long et al.,1982;Cvetkovic et al.,2004;Latham et al., 2013; Zhao, 2013), the natural surface of fractures is rough,which gives rise to channeling flows(Figueiredo et al.,2016)and/or significant eddies when the Reynolds number is larger than the critical value (Zimmerman et al., 2004; Javadi et al., 2010; Li et al.,2016). The channeling flow induces anisotropically distributed streamlines,forming the flow paths of the fluid/solute(Neretnieks,1987;Dershowitz and Fidelibus,1999;Black et al.,2007;Zhao et al.,2014).The significant eddies are due to the inertial effect caused by the nonlinear flow, which has been verified to decrease the permeability of fractured rock masses to a large extent(Ranjith and Darlington,2007;Liu et al.,2016a). Therefore,for some cases such as contaminant transports and high-pressure packer tests, the fracture surface roughness should be incorporated when estimating the hydraulic properties of fractured reservoirs(Chen et al.,2015; Wang et al., 2016a).

To link the mechanical aperture to the hydraulic aperture,many factors such as contact ratio(Walsh,1981;Hakami,1995),standard deviation of varying aperture over mechanical aperture (Patir and Cheng, 1978), fracture surface roughness (Olsson and Barton,2001), harmonic mean of true aperture (Waite et al.,1999), standard deviation of local slope of fracture surface(Xiong et al.,2011),root mean square of the first deviation of the profile (Li and Jiang,2013), and standard deviation of mechanical aperture during shearing(Xie et al.,2015)have been considered.However,there is no empirical function that has negligibly small deviations.The most simple and effective way to simulate the geometries of a natural fracture is to establish a model that is formed by rough walls. The void spaces between the two walls come from fracture deformations subjected to a shear displacement and/or a normal displacement (Olsson and Barton, 2001; Koyama et al., 2006), as well as rock-fluid interactions (Yasuhara et al., 2015). When the local aperture of fractures is known from the in situ measurement or from empirical formula (Lapcevic et al.,1999), a “single-roughwalled model” is used, in which the lower wall of the fracture is smooth and the upper wall is rough(Javadi et al., 2010). However,this is a pseudo three-dimensional (3D) model, because the local fracture apertures are assigned to the elements of two-dimensional(2D) fracture plane at a certain interval for simplifying the model(Brown, 1987; Koyama et al., 2006; Matsuki et al., 2006; Huang et al., 2017a, b). The geometrical characteristics of void spaces between two fracture walls are thus implicitly represented by the anisotropic aperture distribution(Yeo et al.,1998;Watanabe et al.,2009; Vilarrasa et al., 2011). This simplification does not have a significant effect on calculation time while it yields a modeling condition far from reality.

To describe the two-rough-walled fracture, a “double-roughwalled model”is defnied in the present study.Such complex rough fracture models have been established for 2D and/or 3D fracture models (Zimmerman et al., 2004; Crandall et al., 2010; Zou et al.,2015; He et al., 2016). The 2D fracture model is used to calculate fluid flow with two walls being separated without contacts,because the presence of any contact will block the fluid flow(Zhou et al., 2018), which significantly deviates from the natural cases.There are two kinds of 3D double-rough-walled models that have been proposed for modeling nonlinear flow by solving the NS equations. The first one establishes a rough fracture wall and then adds a normally distributed aperture to each element to form the other rough wall (Wang et al., 2016a). Although anisotropic aperture distribution is considered in this model, it is, to some extent,different from the shear-induced void spaces.The other is to create two well-matched rough fracture walls and then artificially moves one wall to the other wall, for example, by assigning a shear displacement of 1 mm and a normal displacement of 0.65 mm(Zou et al., 2017). This simplified normal displacement-shear displacement relationship is difficult to represent the dilation of a rough fracture induced by shearing under a constant normal stiffness(CNS)boundary condition.Besides,the shear-induced anisotropy of hydraulic properties such as permeability anisotropy has been investigated for fluid flow in the linear regime(Huang et al.,2018);however, the anisotropic hydraulic properties of fractures during shearing in the nonlinear flow regime have not been studied in previous works, if any.

In the present study,we aim to establish a double-rough-walled fracture model. The two rough walls of a fracture are numerically generated using a modified successive random additions (SRA) algorithm and the aperture distribution is generated during shearing using a mechanistic model. Then the FLUENT module of ANSYS software using a finite volume method(FVM)code is used to solve the fluid flow in the aperture-anisotropically-distributed void spaces of a fracture by directly solving the NS equations. The difference between the double-rough-walled and the single-roughwalled models for simulating fluid flow is quantitatively illustrated. Finally, the influences of hydraulic gradient on anomalous streamlines and nonlinear flow properties,as well as shear-induced anisotropic hydraulic properties, are investigated.

2. Generation of a 3D rough fracture surface

A key feature of natural rough fracture surfaces is their complex topography following a self-affine fractal distribution (Power and Tullis, 1991; Plouraboue et al., 2000; Mourzenko et al., 2001;Drazer and Koplik, 2002; Jin et al., 2017; Li et al., 2018). This implies that the surface roughness does not contain an intrinsic wavelength and exhibits a repeating pattern at different scales.The self-affine fracture surface can be simulated using the fractional Brownian motion (fBm) (Mandelbrot and Van Ness, 1968). The asperity height of a rough fracture wall in a 2D fBm is described using a continuous and single-valued function Z (x, y). The stationary increment, Z (x + lx, y + ly) - Z (x, y), has a Gaussian distribution with a zero mean over a distance l =The statistical self-affine property of fBm is expressed as

where 〈·〉 is the mathematical expectation; H denotes the Hurst exponent that varies from 0 to 1 and is linearly correlated with the fractal dimension(D)by D=3-H(Family and Vicsek,1991);r is a constant; and σ2is the variance, which is defined as

Several approaches have been proposed to generate the related fBm, including the Weierstrass-Mandelbrot function, successive random additions (SRA) algorithm, and Fourier transformation(Madadi and Sahimi,2003).The SRA is considered to be an efficient and straightforward algorithm. However, the traditional SRA algorithms suffer from poor correction in the stochastic fractal distribution (McGaughey and Aitken, 2000). Thus, a more rigorous and accurate SRA algorithm is developed by Liu et al. (2004), which is used in the present study. A self-affine rough fracture surface that has a length of 14 mm and a width of 10 mm is generated by assigning H=0.5,as shown in Fig.1a.For natural rock fractures,H varies generally between 0.45 and 0.85 (Babadagli et al., 2015). A smaller H value indicates a lower spatial correlation and a rougher fracture surface.

Fig.1. (a)A rough fracture surface;(b)The calculated and averaged JRC values of 51 cutting lines;(c,d)Aperture distributions of single-rough-walled model with us=2.2 mm and 2.6 mm, respectively; (e, f) Aperture distributions of double-rough-walled model with us = 2.2 mm and 2.6 mm, respectively; (g) Enlarged view of meshing; and (h) Aperturefrequency curves for fractures at different shear displacements.

The joint roughness coefficient (JRC) is typically used for characterizing the roughness of fracture surface (Barton and Choubey,1977). According to the suggested method by the International Society for Rock Mechanics and Rock Engineering(ISRM,1981),the JRC of a 3D single fracture should be the average value of 2D cutting lines along the shear direction(Indraratna et al.,2015;Wang et al.,2016b). As shown in Fig.1b, although the JRC values of 2D cutting lines are distributed with a slight fluctuation,the calculated JRC of the fracture surface is 16.77, which is close to the upper bound of JRC value equal to 20.

3. Generation of a 3D double-rough-walled fracture model

The stresses applied on the fractures can be changed during earthquakes and/or tectonic movement, resulting in movement/shearing of faults/fractures(Xue et al.,2013).The void spaces in the natural rough fractures are changed due to the shear-induced dilation, which can be represented by the anisotropically distributed apertures(Lapcevic et al.,1999).The effects of aperture distribution and its change due to normal/shear displacement on the hydraulic properties of single fractures have been quantified experimentally and numerically(Olsson and Barton,2001;Li et al.,2008; Huang et al., 2019).

The hydraulic responses of fractures are significantly influenced by the relative displacement of the two walls of a fracture(Lee and Cho, 2002). In this context, the two rough walls of a fracture are assumed to be well-matched with zero apertures before shearing.The shear process is simulated by moving the upper surface,while the lower surface is fixed. The height of the lower surface is represented using the function ZL(x,y)generated by the modified SRA algorithm. The location of the upper surface ZU(x, y) at each usis estimated as ZU(x, y) = ZL(x, y) + uv, where uvrepresents the normal displacement. Indraratna et al. (2015) proposed an analytical model, uv= ∫˙vdus, to calculate uvduring shearing under CNS boundary conditions, where ˙v is the dilation rate that can be calculated as

where us-peakis the peak shear displacement,c0is the ratio of usto us-peakat which dilation initiates,c1and c2are the decay constants,and ˙vpeakis the peak dilation rate. The variation in uvversus usaddresses the effect of asperity damage along with other governing parameters such as JRC, joint compressive strength (JCS), initial normal stress(σn0),and boundary normal stiffness(Kn).The model parameters used in this study are tabulated in Table 1(after Huang et al.,2018).The aperture distribution b(x,y)at shear displacement us, in which the upper wall of the fracture is sheared along the negative side of x-axis while the lower wall of the fracture is fixed,can be calculated as

Two different fracture models with the same aperture distribution are generated. One is the “single-rough-walled model”, in which one of the fracture walls is a smooth plane and the other one has a rough surface by adding the aperture distribution on the basis of the smooth plane.The other model is the“double-rough-walled model”, in which the same aperture distribution as that in the single-rough-walled model is added on the lower rough plane to form the upper rough plane.Thus,the two walls of the fracture are both rough. As the shearing continues, the effective length of the upper and lower surfaces facing to each other decreases, thus across-sectional area of 10 mm×10 mm(x=4-14 mm and y=0-10 mm) in the xy plane is extracted from the original model to maintain a constant analyzed area.By solving Eq.(5),the aperture distributions at us=2.2 mm and 2.6 mm are obtained to represent the models containing a large amount of isolatedly distributed contact and a small amount of concentrated contact, respectively.Fig. 1c and d exhibits the single-rough-walled models at shear displacements of 2.2 mm and 2.6 mm, respectively, while Fig. 1e and f displays the double-rough-walled models at shear displacements of 2.2 mm and 2.6 mm, respectively. Comparisons between the aperture-frequency curves show a marked change in the aperture heterogeneity (Fig.1h). The heterogeneous aperture field implies a variable local permeability, which in turn impacts the macroscopic flow through the fractures.Larger shear displacement results in a smaller frequency value at aperture of 0 mm (contact)and larger void spaces, and thereafter, the larger permeability.

Table 1 Input parameters used for calculation of normal displacement during shearing(after Huang et al., 2018).

The four models shown in Fig.1c-f are meshed using ICEM CFD(an ANSYS module,ANSYS 13.0),in which the largest length of each tetrahedral cell is 0.1 mm while the spacing of the elements is 0.2 mm,as shown in Fig.1g.The total numbers of tetrahedral cells of the four models are 369,803, 595,016, 389,226 and 619,273,respectively. Then, these meshed files are imported into FLUENT,another ANSYS module (ANSYS 13.0), which is then used to calculate fluid flow by solving the NS equations. Two boundary conditions are estimated: unidirectional flow in x- (Fig.1e) and ydirection (Fig. 1f). To investigate the nonlinear flow properties of fluids through rough-walled fractures,a large variation in pressure drop varying from 0.1 Pa to 105Pa is applied between the inlet and outlet boundaries, resulting in a dimensionless hydraulic gradient(J) in the range of 10-3-103. Here, J is defined as the ratio of the hydraulic head difference applied on the opposing sides of the model to the model length along the fluid flow direction.The other two boundaries are impermeable, if fluid flows from the opposing two boundaries. The fluid is assumed to be water with density of 998.2 kg/m3at the room temperature (20 C), and its viscosity is 0.001 Pa s.

4. Impact of hydraulic gradient on anomalous streamlines and nonlinear flow properties

The Reynolds number defined as the ratio of inertial force to the viscous force is typically utilized to characterize fluid flow. However, in the fractured rock masses, there are a great number of fractures and it is rather difficult to calculate the Reynolds number of each fracture.Instead,J has been increasingly adopted,since it is also a dimensionless parameter and commonly has a known value in the hydraulic tests of fractured rock masses (Liu et al., 2016b;Wang et al., 2016a). Fig. 2a and b presents the x-directional anomalous streamlines for double-rough-walled fracture model using a total of 456 particles at us=2.2 mm with J=10-3and 103,respectively. Fig. 2c and d depicts the x-directional anomalous streamlines for double-rough-walled fracture model using a total of 830 particles at us=2.6 mm with J=10-3and 103,respectively.The results show that when J is small (i.e.10-3), although the streamlines are tortuous due to the aperture variation in different locations, the fluid steadily flows from the left to the right boundary.There is no eddy formed.When J is large(i.e.103),many eddies are formed (see Fig. 2b), indicating that the inertial effects cannot be negligible and fluid flow is within the nonlinear flow regime.With increasing usfrom 2.2 mm to 2.6 mm, the number of streamlines increases significantly,since the shear-induced dilation gives rise to an increase in the void spaces to allow fluid to flow. Fig. 2e and f shows the variations in Q/J and r1for x-directional flow with J = 10-3-103. Here, Q is the flow rate that can be calculated using the solved flow velocities in the NS equations as shown in Eq. (6),and r1is the relative error that can be calculated by Eq. (7).

Fig. 2. (a, b) x-directional anomalous streamlines for fractures at us = 2.2 mm with J = 10-3 and 103, respectively; (c, d) x-directional anomalous streamlines for fractures at us = 2.6 mm with J = 10-3 and 103, respectively; (e, f) Variations in Q/J and relative error for x-directional flow with varying J from 10-3 to 103 at us = 2.2 mm and 2.6 mm,respectively. The different color lines in the figure represent the streamlines of different particles.

ρ[?u/?t+(u·?)u] = -?P+?·T +ρf (6)

Prince Asmund dearly loved all outdoor sports and an open-air life, and from his earliest childhood he had longed to live entirely4 in the forest close by

where u is the flow velocity tensor,P is the hydraulic pressure,ρ is the fluid density,f is the body force tensor,t is the time,and T is the shear stress tensor.

where(Q/J)drand (Q/J)srare the values of Q/J calculated by solving fluid flow through double-rough-walled and single-rough-walled fractures, respectively. Note that in Eq. (7), (Q/J)dris used as a benchmark, because the double-rough-walled fracture model has the same geometry and void spaces of fractures with the natural cases. Therefore, r1is used to quantify the relative difference between the double-rough-walled and single-rough-walled models when characterizing fluid flow through fractures.

The fluid flow in the linear regime obeys the cubic law as shown in Eq.(8),and Q/J is a constant.The cubic law is derived from the NS equations by assuming that the flow velocity is sufficiently small(Zimmerman and Bodvarsson, 1996). With the increase in J, the inertial effect cannot be negligible,resulting in a decrease in Q/J(Liu et al., 2018).

where μ is the viscosity,g is the gravitational acceleration,w is the width of the fracture,and e is the hydraulic aperture.

As shown in Fig. 2e, with increasing J from 10-3to 103, Q/J for both single-rough-walled and double-rough-walled models changes slightly (i.e. J ≤10-1) and then significantly (i.e. J ≥100).This indicates that the fluid flow transits from a linear flow regime to a nonlinear flow regime,in which the critical hydraulic gradient that quantifies the onset of nonlinear flow is in the range of 10-1-100.According to the variations of the curve slopes,three zones are detected: linear flow zone (J = 10-3-10-1), weak nonlinear flow zone(J = 10-1-101)and strong nonlinear flow zone(J = 101-103).The value of Q/J calculated using the single-rough-walled model is larger than that calculated using the double-rough-walled model.This is because in the single-rough-walled model,the lower surface is straight and smooth,which needs a shorter distance for particle to flow, compared with the double-rough-walled model in which the lower surface is rough. Their difference represented by r1as shown in Eq. (7) holds an approximately constant value when the fluid flow is in the linear regime(i.e.J ≤10-1)and increases when the fluid flow enters the nonlinear flow regime (i.e. J ≥100). The single-rough-walled model overestimates the value of Q/J by approximately 14.99%-19.77% at us= 2.2 mm. The double-roughwalled model can more accurately simulate fluid flow through fractures than the single-rough-walled model. As usincreases to 2.6 mm, although the variations in Q/J and r1as shown in Fig. 2f have the same trends with those for us=2.2 mm(see Fig.2e),the value of Q/J is much larger due to the shear-induced larger void spaces.Fig.3 presents the variations in y-directional streamlines,Q/J,and r1for fractures at us=2.2 mm and 2.6 mm,in which a total of 942 and 1377 particles are used, respectively. They have the same variation trends with those as shown in Fig.2 and will be used for analyzing the anisotropic flow properties in the following section.

Fig. 3. (a, b) y-directional anomalous streamlines for fractures at us = 2.2 mm with J = 10-3 and 103, respectively; (c, d) y-directional anomalous streamlines for fractures at us =2.6 mm with J =10-3 and 103, respectively;(e, f) Variations in Q/J and r1 for y-directional flow with varying J from 10-3 to 103 at us = 2.2 mm and 2.6 mm,respectively.The different color lines in the figure represent the streamlines of different particles.

Fig. 4. Variations in Qy/Qx and r2 versus J for both single-rough-walled and double-rough-walled models at shear displacements of (a) 2.2 mm, and (b) 2.6 mm.

5. Shear-induced anisotropic hydraulic properties

Fig.4a and b exhibits the variations in Qy/Qx and r2 with varying J from 10-3 to 103 for fractures at us = 2.2 mm and 2.6 mm,respectively.Here,Qx and Qy are the flow rates in x-and y-direction,respectively, and r2 is defined as

where (Qy/Qx)drand (Qy/Qx)srare the values of Qy/Qxcalculated by solving fluid flow through double-rough-walled and single-roughwalled fractures, respectively.

The results show that the single-rough-walled model also overestimates the value of Qy/Qx,yet the overestimated degree(r2)decreases from 8.26%-12.01% to 0.91%-4.77% when usincreases from 2.2 mm to 2.6 mm.The increase in shear-induced void spaces decreases the effect of fracture surface roughness on the flow properties, which results in the close values of Qy/Qxfor the two models.In the linear flow regime(i.e.J ≤10-1),the effect of inertial force is negligible and Qy/Qxholds a constant value.In the nonlinear flow regime (i.e. J ＞10-1), the effect of inertial force plays a nonnegligible role and Qy/Qxvaries with increasing J. At us= 2.2 mm,with the increment of J from 10-1to 103,Qy/Qxdecreases from 1.11 to 1.04 and from 1.03 to 0.94 for the single-rough-walled and double-rough-walled models, respectively. This indicates that the inertial effect reduces Qymore significantly than Qxdue to the shear-induced contact and local aperture distributions. After increasing usfrom 2.2 mm to 2.6 mm,the contact and local aperture distributions are changed. As a result, Qy/Qxincreases from 1.13 to 1.24 and from 1.12 to 1.19 for the single-rough-walled and doublerough-walled models, respectively. In this case, the inertial effect reduces Qyless significantly than Qx,which plays an opposite role in that at us= 2.2 mm. For the double-rough-walled model that is similar to the natural cases,as usincreases from 2.2 mm to 2.6 mm,the value of Qy/Qxin the linear flow regime(i.e.J=10-3)increases from 1.04 to 1.12.This illustrates that the larger usgives rise to the stronger anisotropy of fluid flow, because during shearing, one surface will climb on the other one, resulting in collision between adjacent fracture asperities and thus generating stripy ridges of contact perpendicular to the shear direction.These stripy ridges of contacts block the flow along the shear direction but have little effects on the flow in the direction perpendicular to the shear direction (Huang et al., 2018). In the nonlinear flow regime, the anisotropy of fluid flow represented by Qy/Qxis either enhanced or decreased with respect to that in the linear flow regime,depending on the geometry of shear-induced void spaces.

6. Conclusions

This study proposed a 3D double-rough-walled fracture model that can reproduce the real geometrical features of fractures. The rough fracture surface in the model was generated using a modified SRA algorithm. The aperture distribution was determined by shearing the fracture walls utilizing a mechanistic model.A FLUENT module that was developed using a FVM code was used to simulate the shear-flow processes by directly solving the NS equations.The variations in the anomalous streamlines, the ratio of flow rate to hydraulic gradient, and the anisotropic flow properties are analyzed.

The results show that the double-rough-walled model can improve the accuracy of modeling fluid flow through 3D single fractures by approximately 14.99%-19.77% with respect to the commonly used single-rough-walled model. With increasing hydraulic gradient from 10-3to 103, the anomalous streamlines change from a steady state to a disturbed state, in which several eddies are formed and numerically visualized. As the shear displacement increases from 2.2 mm to 2.6 mm, the ratio of flow rate to hydraulic gradient increases by approximately one order of magnitude. The fluid flow is more anisotropic within fractures subjected to a larger shear displacement, because the large shear displacement generates some void spaces that significantly enhance the fluid flow in the direction perpendicular to the shear direction than that in the shear direction. In the nonlinear flow regime, as the shear displacement increases from 2.2 mm to 2.6 mm,the anisotropy of hydraulic properties(i.e.the ratio of flow rate in the y-direction to that in the x-direction) of fractures is either enhanced or decreased with respect to that in the linear flow regime,which depends on the variations in the void spaces induced by shearing.

The present study proposed a method to establish the doublerough-walled fracture model and to simulate fluid flow by numerically solving the NS equations, in which the rough fracture morphology and the whole dilation process during shearing can be taken into account by incorporating a mechanistic model. Using this method, some in-depth analyses can be subsequently performed by changing the parameters listed in Table 1,which can be obtained from in situ geological survey and/or according to the requirements of practical implications, such as the fracture morphology, the normal stiffness of surrounding rocks, the shear displacement, the initial normal stress, and so on. The doublerough-walled model can also provide a platform for accurate estimate of the channeling flow and solution of transport through rough fractures during shearing.However,the dilation of fractures during shearing is calculated using Eq.(4),which is obtained based on a great number of laboratory-scale experiments. In the future,we will verify whether Eq. (4) and the double-rough-walled fracture model are suitable for reservoir-scale simulations.Besides,the effects of shear-induced asperity degradation and roughness change of fracture surfaces, as well as “Yaw” effect, that make different local contact/aperture conditions, will be investigated.

Declaration of Competing Interest

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

Acknowledgments

This study has been partially funded by National Natural Science Foundation of China(Grant Nos.51979272 and 51709260),and the Natural Science Foundation of Jiangsu Province, China (Grant No.BK20170276). These supports are gratefully acknowledged.