Modeling linear and nonlinear fluid flow through sheared rough-walled joints taking into account boundary stiffness

Computers and Geotechnics 120 (2020) 103452

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Modeling linear and nonlinear fluid flow through sheared rough-walled joints taking into account boundary stiffness

T

⁎

Richeng Liua,b, Changsheng Wangb, Bo Lic, , Yujing Jiangb, Hongwen Jinga a

State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, PR China School of Engineering, Nagasaki University, Nagasaki 852-8521, Japan c Key Laboratory of Rock Mechanics and Geohazards of Zhejiang Province, Shaoxing University, Shaoxing 312000, PR China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Rock joint Shear Stiffness Fluid flow Permeability Nonlinear flow

This study modeled four hundred fluid flow cases through shear-induced aperture fields of a rough-walled joint with self-affine surface characteristics under different boundary stress and stiffness conditions by solving the Navier-Stokes equations. The results show that as the hydraulic gradient increases from 10−8 to 101, the fluid flow experiences changes characterized by a linear regime, a weak nonlinear regime and a strong nonlinear regime. The permeability is a constant in the linear regime and starts to decline when entering the weak nonlinear regime, and the declining rate gradually converges when reaching the strong nonlinear regime. Such transition is poorly indicated by the global Reynolds number alone; it is instead primarily governed by the magnitudes and distribution characters of local Reynolds numbers. The boundary stiffness has linear/nonlinear relations with the permeability, depending on the magnitude of the joint aperture, which is governed by surface roughness, joint compressive strength and mechanical boundary conditions. With increasing boundary stiffness from 0 to 2.0 GPa/m, the ratio of permeability in y-direction to that in x-direction increases from 2.49 to 26.49 by approximately a factor of 10, suggesting that anomalous flow may prevail in sheared zones with stiff surrounding confinements.

1. Introduction Fluid flow in rock joints is of special importance for assessing the performance of underground projects such as CO2 sequestration [1], geothermal energy development [2] and nuclear waste disposal [3]. Fluids flow through joints primarily via connected void spaces that are strongly influenced by mechanical boundary conditions, especially when a shear process presents [4,5]. Accurate modeling of fluid flow in such complex geometric structures of void spaces requires advanced modeling and computational techniques [6,7]. To date, simulating fluid flow through evolving void spaces of fractures during shearing taking into account realistic mechanical boundary conditions is still a challenging task [8–10]. Early works simplified a joint as an idealized smooth parallel plate and presumed that the fluid flow obeys cubic law, which has been confirmed to deviate from experimental measurements on natural joints mainly due to the presence of surface roughness [11,12]. Later, efforts have been devoted to modifying the cubic law on basis of smooth parallel-plate model by introducing different roughness parameters [8,11,13–18]. These modifications can effectively improve the

⁎

precision in assessment while keeping the simplicity of the governing equations (i.e., the framework of the cubic law) and therefore are computationally efficient. They can be easily incorporated into fracture network models and coupled with the governing equations of thermal, mechanical and chemical effects as well as mass transport [19,20]. The preferential flow paths (channeling flow) cannot be characterized if the rough-walled joints with variable local apertures are simplified to the smooth parallel-plate model with a uniform aperture. The inertial forces in the preferential flow paths can be much greater than those in the smooth parallel-plate model under a specific hydraulic gradient, resulting in the onset of nonlinear flow at relatively small Reynolds numbers [21]. Therefore, the exact flow field in the smooth parallelplate model is oversimplified, which significantly affects the reliability in flow and transport assessments [6,22–25]. With the development of computational and modeling techniques, high-resolution models of real void geometries in joints solving the Reynolds equation or NavierStokes equations (NS) are developed to obtain exact flow fields in linear and/or nonlinear flow regimes [8,26–30]. The responses of void spaces to boundary stresses and the subsequent influences on the fluid flow behavior and the associated other mechanisms can be realistically

Corresponding author. E-mail address: [email protected] (B. Li).

https://doi.org/10.1016/j.compgeo.2020.103452 Received 13 October 2019; Received in revised form 4 January 2020; Accepted 11 January 2020 0266-352X/ © 2020 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 1. (a–b) lower and upper surfaces of a fracture, (c–e) extraction process during shearing, (f) aperture field that provides void spaces for fluid flow and the hydraulic boundary condition.

which is only suitable for free rock blocks without external confinements, e.g., movement of free blocks on the surfaces of slopes [39–43]. For a rock joint below the surface, shear-induced dilation is restricted by surrounding rocks, which can be represented by a constant normal stiffness (CNS) boundary condition [44–48]. A large number of studies have investigated the influence of the boundary stiffness on the mechanical behavior of joints during shearing [48–50], yet little attention has been paid to hydraulic responses [8]. Although the linear flow condition is likely to hold in natural groundwater systems, many engineering practices, such as hydraulic

modeled [8,28,31]. To date, however, these models have seldom been coupled with shear processes under realistic mechanical boundary conditions. In a shear process, the joint dilates in the normal direction, and void spaces are constantly evolving due to the combined effects of normal dilation/contraction and relative movement of surfaces in the shear direction [15,32–35]. This process is governed by the mechanical properties of rocks as well as the mechanical boundary conditions and surface roughness [32,36–38]. Early works commonly apply a constant normal load (CNL) boundary condition for ease of experimentation, 2

3

100 100 100 100 100 100 100 100 100 100 150 150 150 150 150 150 150 150 150 150

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

0.50 0.50 0.50 0.50 0.50 2.00 2.00 2.00 2.00 2.00 0.50 0.50 0.50 0.50 0.50 2.00 2.00 2.00 2.00 2.00

σn0 (MPa)

0.00 0.50 1.00 1.50 2.00 0.00 0.50 1.00 1.50 2.00 0.00 0.50 1.00 1.50 2.00 0.00 0.50 1.00 1.50 2.00

Kn (GPa/m)

3.88 3.18 2.72 2.37 2.11 2.81 2.71 2.60 2.53 2.42 3.49 2.81 2.39 1.99 1.68 3.10 2.99 2.89 2.79 2.69

Em (mm)

3.96 3.24 2.74 2.38 2.10 2.83 2.73 2.63 2.53 2.45 4.31 3.50 2.95 2.55 2.25 3.15 3.03 2.91 2.81 2.71

δv (mm)

Note: Jx and Jy are hydraulic gradients along x- and y-directions, respectively.

JCS (MPa)

Case No.

Table 1 Basic parameters and calculated results for numerical simulation cases.

1.75 1.70 1.66 1.64 1.62 1.65 1.65 1.65 1.64 1.65 2.06 1.87 1.67 1.67 1.63 1.70 1.68 1.65 1.64 1.64

StDev (mm)

2.09 5.43 9.41 12.60 16.04 8.96 9.71 10.38 11.41 11.93 1.24 2.25 2.10 2.82 3.75 6.10 7.14 8.29 9.11 9.79

Rc (%)

−8

10 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8

Jx

~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

10 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101

1

−8

10 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8 10−8

Jy

~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

10 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101 101

1

2.16E−10 8.20E−11 2.82E−11 1.33E−11 6.70E−12 3.72E−11 2.95E−11 2.32E−11 1.74E−11 1.51E−11 3.26E−10 1.23E−10 4.45E−11 1.80E−11 1.02E−11 6.85E−11 5.56E−11 4.23E−11 3.45E−11 2.65E−11

1.20E−08 3.92E−09 1.22E−09 4.46E−10 1.40E−10 1.62E−09 1.24E−09 9.16E−10 6.90E−10 5.55E−10 1.85E−08 6.22E−09 2.08E−09 7.21E−10 2.98E−10 3.19E−09 2.51E−09 1.94E−09 1.50E−09 1.13E−09

4.59E−10 2.89E−10 1.87E−10 1.29E−10 9.05E−11 2.04E−10 1.84E−10 1.65E−10 1.50E−10 1.38E−10 5.54E−10 3.47E−10 2.28E−10 1.52E−10 1.10E−10 2.70E−10 2.43E−10 2.20E−10 1.99E−10 1.81E−10

Minimum

Minimum

Maximum

Ky (m2)

Kx (m2)

2.99E−08 1.56E−08 9.03E−09 5.67E−09 3.70E−09 1.01E−08 8.86E−09 7.82E−09 6.96E−09 6.22E−09 3.94E−08 2.01E−08 1.15E−08 7.12E−09 4.68E−09 1.41E−08 1.25E−08 1.10E−08 9.76E−09 8.66E−09

Maximum

2.13E+00 3.52E+00 6.63E+00 9.72E+00 1.35E+01 5.49E+00 6.22E+00 7.11E+00 8.60E+00 9.12E+00 1.70E+00 2.83E+00 5.12E+00 8.45E+00 1.07E+01 3.95E+00 4.38E+00 5.21E+00 5.77E+00 6.82E+00

Minimum

Ky/Kx

2.49E+00 3.98E+00 7.37E+00 1.27E+01 2.65E+01 6.22E+00 7.12E+00 8.54E+00 1.01E+01 1.12E+01 2.13E+00 3.23E+00 5.52E+00 9.88E+00 1.57E+01 4.44E+00 4.98E+00 5.69E+00 6.52E+00 7.65E+00

Maximum

R. Liu, et al.

Computers and Geotechnics 120 (2020) 103452

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

of consequences such as anomalous flow fields with coexisting linear and nonlinear flow components and the formation of eddies at locations with fierce geometric variations that evolve with changing boundary stresses and displacements [57–59]. Therefore, the hydraulic behavior of rock joints in both linear and nonlinear regimes undergoing shear needs to be characterized to improve the understanding on fluid flow in sheared zones, especially those involving human activities. Motivated by the problems mentioned above, this study aims to apply realistic boundary conditions to a rock joint undergoing shear and quantitatively investigate the responses of linear and nonlinear flow behaviors to the mechanical effects. First, we applied an analytical model for assessing the shear-induced dilation of rock joints under CNS boundary conditions to an artificially generated rough-walled rock joint. Then, the validity of the analytical model was verified through comparisons with direct shear test results and necessary parameters involved in the model were obtained. Finally, with the calculated dilation curves under different boundary conditions, a series of void space models were established and fluid flow simulations through these models were conducted by solving Navier-Stokes (NS) equations under a wide range of hydraulic gradient. Following this technical route, we were able to link the mechanical effects to the hydraulic behavior of sheared rock joints in both linear and nonlinear flow regimes. The obtained simulation results revealed the controlling effect of boundary stiffness on the onset of nonlinear flow and the permeability anisotropy. 2. Methodology 2.1. Joint surface generation The numerical model of rough-walled joints can be established either by scanning and digitizing natural joint surfaces or by numerically generating surfaces following specific mathematical functions that characterize the typical surface morphology [60–62]. Geological surveys indicate that surfaces of natural rock joints can exhibit statistical self-affine properties [63–65]. Different numerical methods have been developed to characterize such properties, based on which artificial joint surfaces are generated for various purposes [25]. One representative method among them is the successive random addition (SRA) with a fractional Brown motion (fBm), which is written as follows [66,67]:

〈Z (x + rl x , y + rl y ) − Z (x , y ) 〉 = 0,

(1)

σrl = r H σl,

(2)

where <∙> is the mathematical expectation, r is a constant, H is the Hurst exponent, Z is the asperity height and σ is the square root of the variance of Z. The distribution of the asperity height increment [Z (x + rl x , y + rl y ) − Z (x , y )] over distance l x2 + l y2 follows a Gaussian distribution [68]. Later, by adding additional random numbers to each existing point, Liu et al. [69] developed a new SRA algorithm to improve the poor scaling and correlation of traditional SRA algorithms. This modified SRA algorithm has been successfully applied to surface characterization of different rock types in previous studies [70–72]. For natural rock joints, H varies from 0.45 to 0.85 according to the rock types and geological formations where a smaller H corresponds to a rougher surface [73–75]. This algorithm was compiled on MATLAB, based on which the two surfaces of a joint with a length of 220 mm and a width of 100 mm were generated, as shown in Fig. 1(a) and (b). Here, H = 0.45 was selected to represent the highly rough joint surfaces existing in nature in order to highlight the difference between roughwalled joints and smooth parallel plates.

Fig. 2. (a) Joint surface profile with a JRC range of 12–14; (b) and (c) Comparison of tested and predicted evolutions of normal displacement on joint samples with two different strengths.

fracturing and water cycling in geothermal energy development have applied hydraulic pressures far greater than the natural level, which may trigger the onset of nonlinear flow in jointed rock masses [21,51–56]. With the presence of nonlinearity in flow, the impacts of shear on the fluid flow become more complex, which produces a series

2.2. Estimation of shear-induced dilation In the initial state before shear, two surfaces of the joint were well4

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 3. Evolution of normal displacement (dilation) during shearing of the artificial joint.

where v̇ is the dilation rate, v̇peak is the peak dilation rate, us − peak is the shear displacement at the peak shear stress, c0 is us / us − peak , at which dilation begins, and c1 and c2 are two decay constants. v̇peak can be determined by laboratory experiments and is correlated with H, us − peak , σn0, the joint compressive strength (JCS), the initial joint normal stiffness at zero normal stress (kni), and the maximum closure of joints (Vm). Here, JCS is estimated by the Schmidt rebound hammer test, and the initial normal stiffness and the maximum closure are determined by fitting the loading-unloading curves with Bandis’ hyperbolic function [77,78]. For details regarding the calculation and validation processes of this model, the reader is referred to Indraratna et al. [76].

mated representing a fresh tensile joint (see Fig. 1(c)). Under an initial normal stress σn0 and a boundary stiffness kn, the upper block was forced to move horizontally by a shear stress στ, while the lower block was fixed [32,45]. Then, cutting off the edges that were not in contact with each other, a joint with a size of 200 mm in length and 100 mm in width was extracted, as shown in Fig. 1(d) and (e). After shearing, the joint was dilated by a normal displacement δv, and the normal stress σn became greater than σn0 under CNS boundary conditions. Here, the normal displacement was calculated by an analytical model that characterizes dilation behavior of joints at different shearing stages as follows [76]:

δν =

∫0

us

ν ̇ dus

2.3. Direct shear test

for < 0(us / us − peak ) ⩽ c0 ⎧0 ⎪ u 1 ⎡1 − ⎪∫ s ν̇ for < c0 (us / us − peak ) ⩽ 1.0 peak ⎢ (c0 − 1)2 ⎪ 0 ⎣ ⎪ 2 us ⎪ − 1 ⎤ dus us − peak , = ⎥ ⎦ ⎨ ⎪ us for (us / us − peak ) > 1.0 ⎪ ∫0 ν ̇ peak exp ⎪ c2 ⎪ ⎧ ⎡ us ⎫ ⎪ ⎨− c1 us − peak − 1 ⎤ ⎬ dus ⎣ ⎦ ⎭ ⎩ ⎩

(

The permeability of sheared rock joints is largely governed by their dilation behavior during shear. Although the analytical model mentioned above has been validated by previous studies, it is necessary to confirm its suitability to the problems concerned here. The validity of the dilation behavior estimated by the analytical model was verified by comparing its predictions to direct shear test results. Cast plaster samples were selected over natural rock joints to facilitate the investigation on sensitivity of dilation behavior to stiffness boundary conditions with a certain surface morphology. The JRC profile defined by Barton and Choubey [39] with a value of 12–14 was selected to represent roughwalled joint surfaces. All samples are 100 mm in width, 200 mm in length and 100 mm in height. Samples with two different strengths

)

(

σn = σn0 + kn × δ v

)

(3) (4) 5

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 4. (a) ~ (d) frequency – local aperture curves under different boundary conditions, and (e) relation between Rc and kn.

were prepared to examine the effect of JCS. The first kind is a mixture of plaster, water and retardant with a weight ratio of 1:0.2:0.005 that has a JCS of 68 MPa. The second kind is a mixture of plaster, water, sand and retardant with a weight ratio of 1:1:0.28:0.005 that has a JCS of 25.6 MPa. Shear tests were performed on a servo-controlled direct shear apparatus capable of shear tests under various stiffness boundary conditions [46,49]. The testing procedure followed the ISRM suggested

method for laboratory determination of shear strength of rock joints [77]. Tests were conducted under an initial normal stress of 2 MPa, with three boundary normal stiffness values: 0 GPa/m (CNL), 3 GPa/m and 7 GPa/m. 2.4. Fluid flow simulation By applying a shear displacement of 20 mm and the calculated 6

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 5. Streamline distributions along x- and y-directions under different J and kn when JCS = 100 MPa and σn0 = 0.5 MPa. Each color line represents a streamline.

cross-sectional area A. The permeability K is defined as

dilation to the upper surface under a specific boundary condition, the aperture field in the joint was obtained, as shown in Fig. 1(f). Hydraulic pressures P1 and P2 were then applied to opposite boundaries while keeping the other sides impermeable. Applied hydraulic gradient J defined as (P1 - P2)/(ρgΔL) varied from 10−8 to 101, which corresponds to the hydraulic pressure gradient from 10−4 to 105 Pa/m. Here, ρ is the fluid density, g is the gravitational acceleration, and ΔL is the joint length. The range of J covers most values existing in nature and encountered in engineering practices. The fluid flow along both the x- and y-directions was calculated by solving the NS equations using the ANSYS Fluent software; the equations are written as follows [79]:

∂u + (u∙∇) u⎤ = −∇P + ∇∙T + ρ f ρ⎡ ⎣ ∂t ⎦

K=μ

μu Q ΔL = A ΔP ρgJ

(6)

where μ is the dynamic viscosity. The Reynolds number of the flow was calculated by [12]

Re =

ρue μ

(7)

where e is the mean aperture of a joint. The global Re of a joint was calculated using mean values of u and e of the entire joint, while the local Re at each element (see Fig. 1f) was calculated using the mean value of u and local value of e within the element. A square element with a constant local e was selected in the modeling of the aperture fields to form local parallel plates that allow precise calculation of Re at each element. Table 1 tabulates the values of variables and calculation cases. To

(5)

where u is the flow velocity tensor, T is the shear stress tensor, t is time, and f is the body force tensor. The volumetric flow rate Q in steady-flow state can be calculated using the mean velocity u multiplied by the 7

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 6. Variations in K when J = 10−8 ~ 101, kn = 0 ~ 2.0 GPa/m and σn0 = 0.5–2.0 MPa.

investigate the impact of the boundary stiffness on the fluid flow, five values of kn were applied to the joint, ranging from 0 (CNL boundary condition) to 2 GPa/m. Combined with two values of σn0 and two values of JCS, 20 numerical models with different aperture distributions were generated. Here, the values of 100 MPa and 150 MPa for JCS represent typical compressive strengths of course-grained and finegrained granites, respectively [80,81]. The value of σn0 = 2 MPa represents the in situ stress of shallow rocks, where engineering activities are most intensive and nonlinear flow is most likely to occur. Then, hydraulic gradients spanning ten orders of magnitude were applied to these models along either x or y direction, leading to 400 numerical flow tests in total. The fluid was regarded as incompressible distilled water at 20 °C with a density of 998.2 kg/m3 and a viscosity of 0.001 Pa∙s .

shown in Fig. 3. These results show that dilation is inhibited under a greater initial normal stress and a greater normal stiffness. Taking JCS = 100 MPa and shear displacement = 20 mm as an example, the normal displacement decreases from 3.96 to 2.10 mm by a rate of 46.97% when σn0 = 0.5 MPa and from 2.83 to 2.45 mm by a rate of 13.43% when σn0 = 2.0 MPa as kn increases from 0 to 2.0 GPa/m. Dilation is more sensitive to normal stiffness under a smaller initial normal stress, since a greater value of dilation under a small initial normal stress may trigger stronger reactive stresses from surrounding rocks as shear proceeds, which in turn restrict dilation in the residual stage. A harder rock can achieve a greater value of dilation under the same boundary condition compared to a softer rock since asperities of the former are more difficult to shear off, and climbing up along contacting asperities may be more significant.

3. Results and analysis

3.2. Aperture and contact area

3.1. Dilation of joints during shear

A shear displacement of 20 mm was selected to represent the residual shear stage (i.e., sheared joints), at which the resistance of joints to shear becomes almost constant. The shear process will impose damages on the contacting asperities, producing gauge particles of different sizes. The amount of particles depends on the rock strength, boundary stress and surface morphology [82]. Here, since relatively low normal stresses (< 2 MPa) were applied to hard rocks, the asperity damages were considerably small and most particles would stay in the original position as long as they were stressed [83]. Therefore, it is reasonable to presume that they have negligible influences on the fluid flow. Meanwhile, the contacting asperities can deform elastically or be crushed depending on the local stress status. The asperity deformation increases with the increasing stresses acting on the boundary, which can influence the magnitude of local apertures in some contents. Here, relatively small normal stresses (< 2 MPa) were applied to hard rocks (JCS = 100 MPa or 150 MPa), therefore, the resulting local aperture changes due to asperity deformation (< 0.01 mm) is much smaller than

Shear tests under boundary stiffness of 0 GPa/m (CNL), 3 GPa/m and 7 GPa/m were conducted on two kinds of artificial joint samples with different strengths. The tested surface profile and comparison of predicted and tested dilation curves are shown in Fig. 2. The comparison confirms that the shear model can precisely predict the dilation behavior of joints subject to different stress and stiffness boundary conditions from initial to residual shear stages. From the comparison, the values of parameters required to fit the dilation curves were obtained, as presented on the figure. Similar values were also recommended by Indraratna et al. [76], although the tested samples were different, indicating that the model is capable of estimating the dilation of rock joints of different types. Using the validated parameter values (i.e., c0 , c1 and c2 ), the dilation curves of typical granite were calculated by applying JCS = 100 MPa and 150 MPa and σn0 = 0.5 and 2.0 MPa to the analytical model, as 8

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 7. Frequencies of the value of Re when J = 10−8, 10−3 and 101, kn = 0 and 2.0 GPa/m and σn0 = 0.5 and 2.0 MPa for flow along x and y directions.

apparent surface area Rc increases linearly with increasing kn, and the increasing rate for the case of JCS = 100 MPa and σn0 = 0.5 MPa is the greatest (see Fig. 4(e)), because joints are more deformable when subject to a smaller σn with a smaller JCS. The increase rate is almost constant for the case of JCS = 150 MPa, while the initial value is smaller when the joint is subjected to a smaller σn0. These indicate that a greater kn gives rise to a greater σn and therefore a greater joint closure and contact area ratio for sheared joints under CNS boundary conditions, and Em is more sensitive to σn0 than JCS.

the shear-induced dilation (> 1 mm). Besides, the ratio of contacting area to the entire surface is typically at a level of 1% for stressed joints [83]. As a generic study, we primarily quantified the shear-induced dilation while the influences of elastic and plastic deformations of asperities on the magnitude of apertures were not considered. A geometrical model of void spaces at this state was established as presented in Fig. 1, and was meshed using tetrahedral cells with a maximum length of 1.0 mm on the help of an ANSYS module, ICEM CFD (The Integrated Computer Engineering and Manufacturing code for Computational Fluid Dynamics). The meshed files were then imported into another ANSYS module, FLUENT, for fluid flow simulation. The entire aperture field as shown in Fig. 1(f) is composed of local apertures of varying values, the frequency of which is plotted versus the local aperture in Fig. 4(a)–(d). In all cases, the local aperture exhibits a Gaussian distribution, which is consistent with the assumptions of SRA and those reported in the literature [10,84]. When σn0 = 0.5 MPa, the mean local aperture Em decreases from 3.88 mm to 2.11 mm by a rate of 45.62% when JCS = 100 MPa and from 3.49 mm to 1.68 mm by a rate of 51.86% when JCS = 150 MPa as kn increases from 0 to 2.0 GPa/m, whereas for σn0 = 2.0 MPa, Em decreases in a relatively small range, from 2.81 mm to 2.42 mm by a rate of 13.88% when JCS = 100 MPa and from 3.10 mm to 2.69 mm by a rate of 13.23% when JCS = 150 MPa. The frequency distributions are scattered when σn0 = 0.5 MPa, but are narrowed down to a compact curve when σn0 = 2.0 MPa. By treating the meshes that have a local aperture less than 0.2 mm as contacts, the ratio of the contacted area to the total

3.3. Permeability evolution from linear to nonlinear flow regimes Fig. 5 shows the streamline distributions under different J and kn. Contacted elements are represented by solid white squares. With increasing kn, both the number and area of contacted elements increase significantly, resulting in more tortuous streamlines and pronounced channelization. When J changes from 10−8 to 101, the streamlines vary negligibly in the regions where the local apertures do not change significantly as shown in the dashed yellow rectangular in Fig. 5(a) and (c), whereas the streamlines are robustly disturbed by a dramatic local aperture change such as contacts and the formation of eddies is identified as shown in the dashed white rectangular in Fig. 5(b) and (d). The evolutions of Kx and Ky with increasing J are shown in Fig. 6. When J is less than 10−4, the permeability of each model holds a constant value independent of J because the inertial force is far smaller than the viscous force in the entire flow field and a linear relationship between 9

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 8. Relation between K and kn when J = 10−8 and J = 101.

Fig. 9. Variations in Ky/Kx with J = 10−8 ~ 101, kn = 0 ~ 2.0 GPa/m and σn0 = 0.5–2.0 MPa.

10

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

Fig. 10. Relation between Ky/Kx and kn when J = 10−8 and J = 101.

hydraulic pressure and flow rate holds. The permeability starts to drop when J is greater than 10−4, and the dropping rate increases with increasing J. In this weak nonlinear regime, a small portion of eddy flow at locations with dramatic geometric variations starts to emerge, while the viscous force still prevails at most locations. After entering the strong nonlinear regime when J is greater than 10−2, the dropping rate of the permeability is almost a constant on a log-log plot of K vs. J, and eddy flow gradually dominates the entire flow field (Fig. 5). It is obvious that the permeability decreases with increasing values of kn, given that a greater value of kn represents stronger confinement on the joint during shearing and therefore results in a reduced dilation in the residual shear stage. The Reynolds number, which is the ratio of the inertial force to the viscous force, is extensively used in the characterization of flow regimes [26,85]. Here, the local Re at each square element (see Fig. 1) and the global Re for each joint model were calculated, and the fractions of Re when J = 10−8, J = 10−3 and J = 101, which represent linear, weak nonlinear and strong nonlinear regimes, are shown in Fig. 7. As shown in Fig. 6, the permeability is a constant when J is around 10−8, which is selected to represent the linear flow. The permeability changes dramatically when J is around 101 with an almost constant rate representing a strong nonlinear regime. A weak nonlinear flow regime is defined in the range of J = 10−4 ~ 10−2 to bridge the linear and strong nonlinear regimes, in which the decreasing rate of the permeability keeps increasing with the increasing J. Here, J = 10−3 is selected to represent the flow behavior in the weak nonlinear flow regime. Each point in Fig. 7 represents the proportion of elements in a range of an order of magnitude (such as 10−17 ~ 10−16, 10−16 ~ 10−15, 10−15 ~ 10−14, and so on) to the total elements. The maximum local value of Re and the global value of Re are on the order of 10−6 ~ 10−4 when J = 10−8, which is far smaller than typical value for the onset of nonlinear flow. Values of a small fraction of local Re become greater than 1 when J = 10−3, indicating that inertial force starts to transcend viscous force at a few local channels with greater flow rates than the average, which explains why calculated permeability starts to drop in the weak nonlinear regime. In the strong nonlinear flow regime, most

local Re numbers are greater than 1 and the magnitude of the maximum value is greater than 103. All local Re numbers have a single value that is equal to the global Re for a smooth parallel-plate model. In contrast, the Re numbers can distribute in a wide range spanning over 1012 orders of magnitude for rough-walled joint models. Dilation is suppressed under greater values of σn0 and kn, which increases the fraction of small Re numbers in the entire flow field. A global value of Re = 1 has long been used for detecting the onset of nonlinear flow through rock joints [26]. Recent studies have suggested the nonlinear flow may appear at much smaller values of Re, as low as an order of 10−3, for sheared joints [21]. A single global Re number cannot efficiently identify the onset of nonlinear flow. For example, global Re values are below 1 for all cases with J = 10−3 and flow in the x-direction; however, obvious decreases in permeability can be observed in Fig. 6 for these cases. Complex void spaces composed of variable local apertures in sheared joints produce anomalous flow fields that allow formation of eddies at small global Re values. Cases with a majority of local Re numbers far less than 1 could exhibit a nonlinear relationship between the hydraulic pressure and flow rate when large local Re numbers far greater than 1 exist in the system. The flow regime is controlled by the distribution characteristics of local Re numbers and their magnitudes. Particularly, the onset of nonlinear flow is controlled by the magnitude and fraction of large local Re numbers. J = 10−8 and 101 were selected to represent flow characteristics in linear and nonlinear regimes, respectively, and relations between K and kn in these distinct flow regimes are shown in Fig. 8. Permeability for the case of JCS = 150 MPa is greater than that for JCS = 100 MPa because asperities on stronger rock surfaces are more difficult to shear off, yielding greater dilation. Closure of joints proceeds as σn0 increases from 0.5 MPa to 2.0 MPa, thereby decreasing the permeability by approximately 0.5 ~ 1.0 orders of magnitude. When σn0 = 0.5 MPa, the permeability decreases with increasing kn, following exponential functions (see Fig. 8(a) and (b)), whereas when σn0 = 2.0 MPa, the permeability exhibits linear relations with kn (see Fig. 8(c) and (d)). A joint sheared under a relatively small initial normal stress could produce remarkable dilation, which may trigger greater reactive confinements 11

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

effect of roughness by producing more complex geometrical structures of void spaces at a smaller mean aperture, which subsequently complicates the distribution of local Re numbers. Under such circumstances (e.g., kn = 2.0 GPa/m and σn0 = 2.0 MPa), partial local Re numbers may have values far greater than the average, resulting in the onset of nonlinear flow at global Re numbers smaller (e.g., 10−1) than typically expected for joints (i.e., Re = 1). Due to the energy dissipation induced by inertial forces, the permeability of the models with J = 101 is approximately 2 orders of magnitude smaller than that with J = 10−8, indicating that optimum hydraulic gradient values need to be estimated and applied to maintain the conductivity of the joint, which is particularly important in the efficiency of resource extractions. Linear and nonlinear relations between permeability and boundary stiffness exist simultaneously, depending on the magnitude of the joint aperture. For a sheared joint, the permeability along the shear direction is less than that perpendicular to shear direction, and the permeability anisotropy is stronger in the linear regime than that in the strong nonlinear regime. The boundary stiffness significantly enhances permeability anisotropy, e.g., Ky/Kx exponentially increases from 2.49 to 26.49 by approximately a factor of 10, when the initial normal stress is relatively small. This indicates that anisotropic flow may prevail in hard rocks with strong surrounding confinements at shallow depths. Although 400 numerical cases were executed for investigating the mechanical effect on hydraulic behavior, only one joint surface was taken into account in the present study due to the heavy computational burden. In future studies, the effects of the surface roughness on the mechanical and hydraulic behaviors of joints subject to shear need to be investigated. Applying greater magnitudes of normal stress and stiffness is also important to fully unmask different mechanical and hydraulic responses between shallow and deep-seated rock joints. Besides, the hydraulic pressure with a comparable or greater magnitude to the confining stress can open the joint and enlarge the dilation, which in turn may alter the defined flow regime. Given the small magnitude of hydraulic pressure comparing to the normal stress, the impact of hydraulic pressure on the deformation of joint walls is not considered here. This fully coupled effect and its scale dependency need to be addressed in the future studies.

(i.e., increased normal stress during shearing) from surrounding rocks. Therefore, both initial normal stress and the boundary stiffness contribute to depressed dilation, and the permeability of sheared joints, which is governed by the magnitude of dilation, is more sensitive to kn at a smaller σn0. In the nonlinear flow regime, the inertial forces significantly contribute to energy dissipations. As a result, the permeability of the models with J = 101 is approximately 2 orders of magnitude smaller than that with J = 10−8. 3.4. Anisotropic permeability When reaching the residual stage (i.e., us = 20 mm), asperities on one surface slide over the corresponding asperities on the other surface accompanied by shearing off of asperity tops (Jiang et al., 2006). Since only facing slopes are in contact, stripe shaped contact areas with a long axis parallel to the y-direction are formed and tend to resist fluid flow along the x-direction (Fig. 5), while void spaces formed at front and back of contact areas provide preferable flow paths along the y-direction (Koyama et al., 2004). As a result, Ky is universally greater than Kx. Here, we use Ky/Kx to represent the permeability anisotropy, as plotted in Fig. 9. Ky/Kx has constant values in both the linear and strong nonlinear regimes and decreases in the weak nonlinear regime. In the linear flow regime, Ky/Kx is not correlated with J because both Ky and Kx are constants. In the strong nonlinear regime, both Kx and Ky undergo the same rate of decrease, as shown in Fig. 6, and therefore Ky/Kx is also a constant. Given the greater value of Ky than Kx, the flow along the y-direction may enter the weak nonlinear regime at smaller values of J than that along x-direction as indicated in Fig. 7, which explains the slight drops of Ky/Kx in the weak nonlinear regime. Ky/Kx increases exponentially with the increment of kn when σn0 = 0.5 MPa yet increases linearly when σn0 = 2 MPa (see Fig. 10). Taking J = 10−8 and JCS = 100 MPa as an example, Ky/Kx exponentially increases from 2.49 to 26.49 by approximately a factor of 10 when σn0 = 0.5 MPa and linearly increases from 6.22 to 11.21 by approximately a factor of 2 when σn0 = 2.0 MPa. The change in Ky/Kx has an opposite trend with δv, indicating that a greater value of joint aperture corresponds to a smaller permeability anisotropy. This is because the effects of roughness and contacts diminish with increasing magnitude of the aperture, resulting in more uniform flow fields [28]. The results show that Ky/Kx = 1.70–6.22 when kn = 0, which agree well with the ranges reported in the literature where the boundary stiffness is not considered [72,86,87]. In contrast, when kn is greater than 0 (i.e., CNS boundary condition), Ky/Kx may reach values as large as 26.49 at a small σn0. This level of anisotropy has not been reported before, which reveals the strong impact of the boundary stiffness on permeability anisotropy.

CRediT authorship contribution statement Richeng Liu: Conceptualization, Data curation, Writing - original draft. Changsheng Wang: Validation, Software. Bo Li: Methodology, Writing - review & editing. Yujing Jiang: Supervision, Writing - review & editing. Hongwen Jing: Visualization, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

4. Conclusions Fluid flow simulations of rock joints with self-affine surface characteristics sheared under constant normal stiffness boundary conditions were performed. The joint permeability evolution under a hydraulic gradient ranging from 10−8 ~ 101 was calculated by solving the NS equations. The effects of the joint compressive strength, boundary stress and stiffness on the linear and nonlinear flow behaviors were investigated. The results reveal that with an increasing hydraulic gradient spanning ten orders of magnitude, the fluid flow evolves from a linear regime when J ≤ 10−4 to a weak nonlinear regime when 10−4 < J ≤ 10−2 and finally reaches a strong nonlinear regime when J > 10−2. The onset of nonlinear flow is governed by the magnitude and distribution of local Re numbers rather than a single global Re number, which merely represents the average ratio of the inertial force to the viscous force. The greater the initial normal stress and the boundary stiffness, the smaller the mean aperture after shear (up to a decreasing ratio of 51.86%). These boundary conditions enhance the

Acknowledgments This study has been partially funded by Natural Science Foundation of China, China (Grant Nos. 51734009, 51979272, 51609136), Natural Science Foundation of Jiangsu Province, China (No. BK20170276), and State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining and Technology, China (No. SKLGDUEK1906). These supports are gratefully acknowledged. References [1] Juanes R, MacMinn CW, Szulczewski ML. The footprint of the CO2 plume during carbon dioxide storage in saline aquifers: storage efficiency for capillary trapping at the basin scale. Transp Porous Med 2010;82(1):19–30. [2] Mora P, Wang Y, Alonso-Marroquin F. Lattice solid/Boltzmann microscopic model

12

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

[3] [4]

[5] [6]

[7] [8]

[9] [10]

[11] [12] [13] [14]

[15] [16]

[17]

[18]

[19]

[20]

[21]

[22] [23] [24]

[25]

[26] [27]

[28] [29] [30]

[31] [32]

[33] [34] [35] [36]

to simulate solid/fluid systems—a tool to study creation of fluid flow networks for viable deep geothermal energy. J Earth Sci 2015;26(1):11–9. Tsang CF, Neretnieks I, Tsang Y. Hydrologic issues associated with nuclear waste repositories. Water Resour Res 2015;51(9):6923–72. Ameli P, Elkhoury JE, Morris JP, Detwiler RL. Fracture permeability alteration due to chemical and mechanical processes: a coupled high-resolution model. Rock Mech Rock Eng 2014;47(5):1563–73. Dang W, Wu W, Konietzky H, Qian J. Effect of shear-induced aperture evolution on fluid flow in rock fractures. Comput Geotech 2019;114:103152. Brush DJ, Thomson NR. Fluid flow in synthetic rough-walled fractures: NavierStokes, Stokes, and local cubic law simulations. Water Resour Res 2003;39(4),. https://doi.org/10.1029/2002WR001346. Javadi M, Sharifzadeh M, Shahriar K. A new geometrical model for non-linear fluid flow through rough fractures. J Hydrol 2010;389(1–2):18–30. Xiong X, Li B, Jiang Y, Koyama T, Zhang C. Experimental and numerical study of the geometrical and hydraulic characteristics of a single rock fracture during shear. Int J Rock Mech Min Sci 2011;48(8):1292–302. Pyrak-Nolte LJ, Nolte DD. Approaching a universal scaling relationship between fracture stiffness and fluid flow. Nat Commun 2016;7:10663. Huang N, Liu R, Jiang Y, Cheng Y, Li B. Shear-flow coupling characteristics of a three-dimensional discrete fracture network-fault model considering stress-induced aperture variations. J Hydrol 2019;571:416–24. Witherspoon PA, Wang JSY, Iwai K, Gale JE. Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res 1980;16(6):1016–24. Zimmerman RW, Bodvarsson GS. Hydraulic conductivity of rock fractures. Transp Porous Med 1996;23(1):1–30. Renshaw CE. On the relationship between mechanical and hydraulic apertures in rough-walled fractures. J Geophys Res 1995;100(B12):24629–36. Waite ME, Ge S, Spetzler H. A new conceptual model for fluid flow in discrete fractures: an experimental and numerical study. J Geophys Res 1999;104(B6):13049–59. Olsson R, Barton N. An improved model for hydromechanical coupling during shearing of rock joints. Int J Rock Mech Min Sci 2001;38(3):317–29. Rasouli V, Hosseinian A. Correlations developed for estimation of hydraulic parameters of rough fractures through the simulation of JRC flow channels. Rock Mech Rock Eng 2011;44(4):447–61. Wang L, Cardenas MB, Slottke DT, Ketcham RA, Sharp Jr. JM. Modification of the Local Cubic Law of fracture flow for weak inertia, tortuosity, and roughness. Water Resour Res 2015;51(4):2064–80. Xie LZ, Gao C, Ren L, Li CB. Numerical investigation of geometrical and hydraulic properties in a single rock fracture during shear displacement with the NavierStokes equations. Environ Earth Sci 2015;73(11):7061–74. Zhao Z, Li B, Jiang Y. Effects of fracture surface roughness on macroscopic fluid flow and solute transport in fracture networks. Rock Mech Rock Eng 2014;47(6):2279–86. Liu R, Jiang Y, Li B, Wang X. A fractal model for characterizing fluid flow in fractured rock masses based on randomly distributed rock fracture networks. Comput Geotech 2015;65:45–55. Javadi M, Sharifzadeh M, Shahriar K, Mitani Y. Critical Reynolds number for nonlinear flow through rough-walled fractures: the role of shear processes. Water Resour Res 2014;50(2):1789–804. Ge S. A governing equation for fluid flow in rough fractures. Water Resour Res 1997;33(1):53–61. Zimmerman RW, Yeo IW. Fluid flow in rock fractures: from the Navier-Stokes equations to the cubic law. Geophys Monogr-Am Geophys Union 2000;122:213–24. Koyama T, Neretnieks I, Jing L. A numerical study on differences in using NavierStokes and Reynolds equations for modeling the fluid flow and particle transport in single rock fractures with shear. Int J Rock Mech Min Sci 2008;45(7):1082–101. Wang L, Cardenas MB. Development of an empirical model relating permeability and specific stiffness for rough fractures from numerical deformation experiments. J Geophys Res 2016;121(7):4977–89. Zimmerman RW, Al-Yaarubi A, Pain CC, Grattoni CA. Non-linear regimes of fluid flow in rock fractures. Int J Rock Mech Min Sci 2004;41:163–9. Vilarrasa V, Koyama T, Neretnieks I, Jing L. Shear-induced flow channels in a single rock fracture and their effect on solute transport. Transport Porous Med 2011;87(2):503–23. Kang PK, Brown S, Juanes R. Emergence of anomalous transport in stressed rough fractures. Earth Planet Sci Lett 2016;454:46–54. Kong B, Chen S. Numerical simulation of fluid flow and sensitivity analysis in rough-wall fractures. J Petrol Sci Eng 2018;168:546–61. Marchand S, Mersch O, Selzer M, Nitschke F, Schoenball M, Schmittbuhl J, et al. A stochastic study of flow anisotropy and channelling in open rough fractures. Rock Mech Rock Eng 2019. https://doi.org/10.1007/s00603-019-01907-4. Zou L, Jing L, Cvetkovic V. Shear-enhanced nonlinear flow in rough-walled rock fractures. Int J Rock Mech Min 2017;97:33–45. Esaki T, Du S, Mitani Y, Ikusada K, Jing L. Development of a shear-flow test apparatus and determination of coupled properties for a single rock joint. Int J Rock Mech Min Sci 1999;36(5):641–50. Bahaaddini M, Sharrock G, Hebblewhite BK. Numerical direct shear tests to model the shear behaviour of rock joints. Comput Geotech 2013;51:101–15. Fathi A, Moradian Z, Rivard P, Ballivy G, Boyd AJ. Geometric effect of asperities on shear mechanism of rock joints. Rock Mech Rock Eng 2016;49(3):801–20. Thirukumaran S, Indraratna B. A review of shear strength models for rock joints subjected to constant normal stiffness. J Rock Mech Geotech Eng 2016;8(3):405–14. Li Y, Sun S, Tang C. Analytical prediction of the shear behaviour of rock joints with quantified waviness and unevenness through wavelet analysis. Rock Mech Rock Eng

2019. https://doi.org/10.1007/s00603-019-01817-5. [37] Li B, Li Y, Zhao Z, Liu R. A mechanical-hydraulic-solute transport model for roughwalled rock fractures subjected to shear under constant normal stiffness conditions. J Hydrol 2019;579:124153. [38] Liu R, He M, Huang N, Jiang Y, Yu L. Three-dimensional double-rough-walled modeling of fluid flow through self-affine shear fractures. J Rock Mech Geotech Eng 2019. https://doi.org/10.1016/j.jrmge.2019.09.002. [39] Barton N, Choubey V. The shear strength of rock joints in theory and practice. Rock Mech 1977;10(1–2):1–54. [40] Reitsma S, Kueper BH. Laboratory measurement of capillary pressure-saturation relationships in a rock fracture. Water Resour Res 1994;30(4):865–78. [41] Yeo IW, De Freitas MH, Zimmerman RW. Effect of shear displacement on the aperture and permeability of a rock fracture. Int J Rock Mech Min Sci 1998;35(8):1051–70. [42] Grasselli G, Egger P. Constitutive law for the shear strength of rock joints based on three-dimensional surface parameters. Int J Rock Mech Min Sci 2003;40(1):25–40. [43] Oh J, Cording EJ, Moon T. A joint shear model incorporating small-scale and largescale irregularities. Int J Rock Mech Min Sci 2015;76:78–87. [44] Johnston IW, Lam TSK, Williams AF. Constant normal stiffness direct shear testing for socketed pile design in weak rock. Geotechnique 1987;37(1):83–9. [45] Indraratna B, Haque A, Aziz N. Shear behaviour of idealized infilled joints under constant normal stiffness. Geotechnique 1999;49(3):331–55. [46] Jiang Y, Xiao J, Tanabashi Y, Mizokami T. Development of an automated servocontrolled direct shear apparatus applying a constant normal stiffness condition. Int J Rock Mech Min Sci 2004;41(2):275–86. [47] Lee YK, Park JW, Song JJ. Model for the shear behavior of rock joints under CNL and CNS conditions. Int Int J Rock Mech Min Sci 2014;70:252–63. [48] Shrivastava AK, Rao KS. Shear behaviour of rock joints under CNL and CNS boundary conditions. Geotech Geol Eng 2015;33(5):1205–20. [49] Jiang Y, Li B, Tanabashi Y. Estimating the relation between surface roughness and mechanical properties of rock joints. Int J Rock Mech Min Sci 2006;43(6):837–46. [50] Li B, Jiang YJ, Koyama T, Jing LR. Experimental study of the hydro-mechanical behavior of rock joints using a parallel-plate model containing contact areas and artificial fractures. Int J Rock Mech Min Sci 2008;45(3):362–75. [51] Ranjith PG, Darlington W. Nonlinear single-phase flow in real rock joints. Water Resou Res 2007;43(9):W09502. [52] Adler PM, Malevich AE, Mityushev VV. Nonlinear correction to Darcy’s law for channels with wavy walls. Acta Mech 2013;224(8):1823–48. [53] Zhang Z, Nemcik J. Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures. J Hydrol 2013;477:139–51. [54] Chen YF, Liu MM, Hu SH, Zhou CB. Non-Darcy's law-based analytical models for data interpretation of high-pressure packer tests in fractured rocks. Eng Geol 2015;199:91–106. [55] Briggs S, Karney BW, Sleep BE. Numerical modeling of the effects of roughness on flow and eddy formation in fractures. J Rock Mech Geotech Eng 2017;9(1):105–15. [56] Liu R, Li B, Yu L, Jiang Y, Jing H. A discrete-fracture-network fault model revealing permeability and aperture evolutions of a fault after earthquakes. Int J Rock Mech Min Sci 2018;107:19–24. [57] Cardenas MB, Slottke DT, Ketcham RA, JrJM Sharp. Navier-Stokes flow and transport simulations using real fractures shows heavy tailing due to eddies. Geophys Res Lett 2007;34(14):L14404. [58] Zou L, Jing L, Cvetkovic V. Roughness decomposition and nonlinear fluid flow in a single rock fracture. Int J Rock Mech Min Sci 2015;75:102–18. [59] Zhou JQ, Wang M, Wang L, Chen YF, Zhou CB. Emergence of nonlinear laminar flow in fractures during shear. Rock Mech Rock Eng 2018;51(11):3635–43. [60] Brown SR. Simple mathematical model of a rough fracture. J Geophys Res 1995;100(B4):5941–52. [61] Glover PWJ, Matsuki K, Hikima R, Hayashi K. Fluid flow in synthetic rough fractures and application to the Hachimantai geothermal hot dry rock test site. J Geophys Res 1998;103(B5):9621–35. [62] Hong ES, Lee IM, Lee JS. Measurement of rock joint roughness by 3D scanner. Geotech Test J 2006;29(6):482–9. [63] Power WL, Tullis TE. Euclidean and fractal models for the description of rock surface roughness. J Geophys Res 1991;96(B1):415–24. [64] Plouraboue F, Kurowski P, Boffa JM, Hulin JP, Roux X. Experimental study of the transport properties of rough self-affine fractures. J Contam Hydrol 2000;46(3–4):295–318. [65] Kulatilake P, Balasingam P, Park J, Morgan R. Natural rock joint roughness quantification through fractal techniques. Geotech Geol Eng 2006;24(5):1181. [66] Mandelbrot BB. The fractal geometry of nature. New York: W.H. Freeman; 1983. [67] Brown SR, Scholz CH. Broad bandwidth study of the topography of natural rock surfaces. J Geophys Res 1985;90(B14):12575–82. [68] Molz FJ, Liu HH, Szulga J. Fractional Brownian motion and fractional Gaussian noise in subsurface hydrology: a review, presentation of fundamental properties, and extensions. Water Resour Res 1997;33(10):2273–86. [69] Liu HH, Bodvarsson GS, Lu S, Molz FJ. A corrected and generalized successive random additions algorithm for simulating fractional Levy motions. Math Geol 2004;36(3):361–78. [70] Ye Z, Liu HH, Jiang Q, Zhou C. Two-phase flow properties of a horizontal fracture: the effect of aperture distribution. Adv Water Resour 2015;76:43–54. [71] Wang M, Chen YF, Ma GW, Zhou JQ, Zhou CB. Influence of surface roughness on nonlinear flow behaviors in 3D self-affine rough fractures: Lattice Boltzmann simulations. Adv Water Resour 2016;96:373–88. [72] Huang N, Liu R, JiangY. Numerical study of the geometrical and hydraulic characteristics of 3D self-affine rough fractures during shear. J Nat Gas Sci Eng 2017;45:127–42.

13

Computers and Geotechnics 120 (2020) 103452

R. Liu, et al.

[81] Rong G, Yang J, Chen L, Zhou C. Laboratory investigation of nonlinear flow characteristics in rough fractures during shear process. J Hydrol 2016;541:1385–94. [82] Zhao Z, Peng H, Wu W, Chen Y. Characteristics of shear-induced asperity degradation of rock fractures and implications for solute retardation. Int J Rock Mech Min Sci 2018;105:53–61. [83] Kling T, Vogler D, Pastewka L, Amann F, Blum P. Numerical simulations and validation of contact mechanics in a granodiorite fracture. Rock Mech Rock Eng 2018;51(9):2805–24. [84] Lapcevic PA, Novakowski KS, Sudicky EA. The interpretation of a tracer experiment conducted in a single fracture under conditions of natural groundwater flow. Water Resour Res 1999;35(8):2301–12. [85] Bear J. Dynamics of fluids in porous media. New York: Elsevier; 1972. [86] Koyama T, Fardin N, Jing L, Stephansson O. Numerical simulation of shear-induced flow anisotropy and scale-dependent aperture and transmissivity evolution of rock fracture replicas. Int J Rock Mech Min Sci 2006;43(1):89–106. [87] Koyama T, Fardin N, Jing L. Shear-induced anisotropy and heterogeneity of fluid flow in a single rock fracture with translational and rotary shear displacements: A numerical study. Int J Rock Mech Min Sci 2004;41(3). 426 426.

[73] Odling NE. Natural fracture profiles, fractal dimension and joint roughness coefficients. Rock Mech Rock Eng 1994;27(3):135–53. [74] Schmittbuhl J, Steyer A, Jouniaux L, Toussaint R. Fracture morphology and viscous transport. Int J Rock Mech Min Sci 2008;45(3):422–30. [75] Babadagli T, Ren X, Develi K. Effects of fractal surface roughness and lithology on single and multiphase flow in a single fracture: an experimental investigation. Int J Multiphas Flow 2015;68:40–58. [76] Indraratna B, Thirukumaran S, Brown ET, Zhu SP. Modelling the shear behaviour of rock joints with asperity damage under constant normal stiffness. Rock Mech Rock Eng 2015;48(1):179–95. [77] Ulusay R. The ISRM suggested methods for rock characterization, testing and monitoring: 2007–2014. Springer; 2014. [78] Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci 1983;20(6):249–68. [79] Neuville A, Flekkøy EG, Toussaint R. Influence of asperities on fluid and thermal flow in a fracture: a coupled Lattice Boltzmann study. J Geophys Res 2013;118(7):3394–407. [80] Singh KK, Singh DN, Ranjith PG. Laboratory simulation of flow through single fractured granite. Rock Mech Rock Eng 2015;48(3):987–1000.

14