Radial fluid flow regime in a single fracture under high hydraulic pressure during shear process

Journal of Hydrology 579 (2019) 124142

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Radial fluid flow regime in a single fracture under high hydraulic pressure during shear process

T

Cheng Caoa, Zengguang Xua, , Junrui Chaia,b, Yaqi Lia ⁎

a b

State Key Laboratory of Eco-hydraulics in Northwest Arid Region of China, Xi’an Univ. of Technology, Xi’an 710048, China School of Civil Engineering, Xijing Univ., Xi’an 710123, China

ARTICLE INFO

ABSTRACT

This manuscript was handled by Huaming Guo, Editor-in-Chief

The single-fracture seepage model is important for fracture network seepage models and underground engineering, but the existing knowledge is insufficient for the roles of the high hydraulic pressure, joint surface destruction, gouge material, and the anisotropy of flow regime in the single-fracture seepage model during shear process. In this study, coupled shear-flow tests were conducted at 0.2, 0.4, 0.6 and 0.8 MPa hydraulic pressure with radial fluid flow, and then 3D simulation models were established by using COMSOL Multiphysics 5.2a with the Reynolds-averaged Navier–Stokes k-ε turbulence model for initial and peak shear conditions. Results showed that the radial flow cubic law diverged from experimental results seriously, and non-Darcy flow behaviour was remarkable during tests. High hydraulic pressure decreased the accuracy of the radial flow cubic law dramatically by increasing the influence of inertia force, and the accuracy was sensitive to the increase in lower hydraulic pressure. The flow regime and velocity were varied and influenced by the anisotropic tortuosity ratio, maximum inclination angle, and the uneven distribution of aperture. The absolute values of velocity vertical to the radial distance reduced with decreasing roughness. The roughness, tortuosity ratio, and maximum inclination angle of the fracture joint reduced dynamically with destruction of the joint surface and formation of gouge material, leading to a significant change in the relationship between mechanical and hydraulic apertures. Moreover, the radial flow cubic law was modified on the basis of the relationship of mechanical and hydraulic aperture with hydraulic gradient. This study may serve as a basis to establish fracture network seepage model and other similar tests.

Keywords: Coupled shear-flow tests Radial flow cubic law Anisotropy of flow regime Joint surfaces destruction Gouge material High hydraulic pressure

1. Introduction The safety and stability of the fluid flow rock mass play a pivotal role in underground water supply, water reservoir engineering, and other fluid-flow-related engineering (Javadi et al., 2010; Rong et al., 2017; Wang et al., 2009; Olayiwola and Dejam, 2019). In a rock mass, fluid flow mainly concentrates in fractures; thus, the study of seepage in a single fracture is rapidly becoming a key instrument in rock mass engineering (Xiong et al., 2011; Mashayekhizadeh et al., 2012; Zhang et al., 2017; Zhang and Nemcik, 2013; Zhang et al., 2019a; Javadi et al., 2010), and the single-fracture seepage theory is fundamental to the fracture network seepage simulation model (Huang et al., 2019; Dejam et al., 2018; Zhang et al., 2019b; Tan et al., 2019). The Navier–Stokes (N–S) equations are widely used to describe fluid flow behaviour and analyse fluid–structure interaction (Javadi et al., 2014; Azarmanesh et al., 2019). However, the N–S equations are difficult to apply to rock mass with countless fractures, and the main case

⁎

of this problem is the inertia force item cannot be efficiently solved (Wang et al., 2015; Rong et al., 2016; Zhou et al., 2015; Cao et al., 2018). Thus, the use of simplified seepage equations, such as the Stokes and Reynolds equations, and the most widely used cubic law is proposed in the exploration of fluid flow evolution in a single fracture and evaluation of fracture hydraulic conductivity (Javadi et al., 2014; Xiao et al., 2016; Li et al., 2008). The cubic law shows that flow rate is linear with hydraulic pressure and the cube of aperture. As shown in Fig. 1(a), the cubic law assumes the joint surfaces are smooth and parallel and the fluid flow in fracture obeys Darcy’s law. However, for natural fracture (Fig. 1(b)) the joint surfaces are always rough, the contact areas are ubiquitous in fracture, the effective seepage channels are narrowed by vortexes and gouge material, and the flow path is usually tortuous. In addition, the flow regime is easy to be turbulent with high hydraulic pressure. Thus, the cubic law cannot provide adequate support when applied to natural aperture (Zhou et al., 2015). Many researchers have shown increasing interest in the flow behaviour in a single fracture with

Corresponding author. E-mail addresses: [email protected] (Z. Xu), [email protected] (J. Chai).

https://doi.org/10.1016/j.jhydrol.2019.124142 Received 13 May 2019; Received in revised form 10 September 2019; Accepted 11 September 2019 Available online 12 September 2019 0022-1694/ © 2019 Published by Elsevier B.V.

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

Nomenclature a b CNL Cε1, Cε2, e e0 eh em F g I J JRC k Pk Q Qc Qe Re ri

ro RV u un VVR α and β

[ML−5T−1], Model coefficients of Forchheimer law [ML−8], Model coefficients of Forchheimer law [–], Constant normal load and Cμ [–], Calculated coefficients of Reynolds-averaged Navier-Stokes k-ε turbulence model [L], Fracture aperture in cubic law [L], Initial mechanical aperture [L], Hydraulic aperture [L], Mechanical aperture [MLT−2], Body force tensor [LT−2], Gravitational acceleration [–], Unit tensor [–], Hydraulic gradient [–], Joint roughness coefficient [ML2T−2], Turbulent kinetic energy [ML2T−2], Turbulent energy generated by average velocity gradient [L3T−1], Flow rate in cubic law [L3T−1], Flow rate calculated from cubic law [L3T−1], Flow rate measured from tests [–], Reynolds number [L], Joint surface outer radius of specimen

δ δp Δh Δun ε θmax μ μT ρ σapert σk σε τ υ P

complex conditions and in improving the accuracy of the single-fracture seepage model. Experiments and numerical simulation play important roles in addressing the difficult issues. Recent experiments have included confining stress-flow tests (Chen et al., 2019; Zhang et al., 2017; Zhang and Nemcik, 2013; Zhou et al., 2015; Rong et al., 2017) and shear-flow tests (Esaki et al., 1999; Javadi et al., 2010; Olsson and Barton, 2001; Rong et al., 2016; Xiong et al., 2011; Yin et al., 2017). Numerical simulation models were usually established with the Reynolds equation, the Stokes equation, the N–S equations, and turbulence models (Zhang et al., 2019a; Javadi et al., 2010; Wang et al., 2015; Xiong et al., 2011; Zhang et al., 2017). Experimental and numerical studies aimed at investigating the influence of the factors in Fig. 1(b) in

[L], Water inlet radius of specimen [LT−1], Radial flow velocity [LT−1], Flow velocity tensor [L], Normal displacement of aperture due to normal stress [LT−1], Flow velocity vertical to the radial distance [–], Coefficients in the fitting formula of mechanical and hydraulic aperture [L], Shear displacement [L], Peak shear displacement [L], Difference in the pressure head of water inlet and outlet [L], Normal displacement of aperture due to shear deformation [–], Dissipation rate [–], Maximum inclination angle [ML−1T−1], Viscous coefficient [ML−1T−1], Turbulent viscous coefficient [ML−3], Fluid density [L], Standard deviation of aperture [–], Turbulent Prandtl number of turbulent kinetic energy [–], Turbulent Prandtl number of dissipation rate [–], Tortuosity ratio [L2T−1], Fluid kinematic viscosity [ML−1T−2], Pressure gradient

the cubic law and modifying the cubic law. Other non-linear seepage models such as Forchheimer or Lzbash laws, have also been widely used in recent years. Li et al. (2008) reported that increased contact ratio evidently reduces the hydraulic conductivity of fracture and the distribution of contact areas is another factor that influences fluid flow in a fracture. The contact ratio of fracture varies during shear. Xiao et al. (2016) established a contact algorithm of artificial marble fracture during shear by performing laboratory tests and using charge-coupled device cameras and found remarkable changes in contact ratio and mean aperture. The algorithm helps build the seepage model but neglects the influence of gouge materials from crushed asperities and may

Joint surface

Flow direction

(a) Vortex Contact area

Joint surface

Flow direction Gouge material

(b)

Fig. 1. Fluid flow in a single fracture: (a) smooth and parallel joint surface, (b) natural rough joint surface. 2

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

overestimate the hydraulic conductivity of a fracture. When fluid flow bypasses contact areas and asperities or climbs over asperities, the seepage path becomes tortuous. Thus, the actual path for fluid flow is always larger than the straight-line distance, and flow velocity changes frequently. Flow velocity, the hydraulic conductivity of fracture, and energy loss are diverse results from anisotropic roughness. Zhang et al. (2019a) reported how the shape and density of roughness elements affect fluid flow in a regular single roughness fracture, and found the remarkable influence of vortex caused when fluid flow bypasses the roughness elements. However, the anisotropy of flow regime cannot be considered because the simulation models in the study of Zhang et al. (2019a) were 2D. After using numerical simulations, Thompson and Brown (1991) found that the important directional characteristic of joint surfaces is degree roughness and that the direction of fluid flow and tha the groove of asperities are influencing factors. Employing experiments and simulation method, Auradou et al. (2006) found that fluid flow velocity and fracture hydraulic conductivity perpendicular to shear are larger than those parallel to shear resulting from the change in aperture and the contact of void space, and the displacement of joint surfaces considerably contributes to the effect. By solving Reynolds equation, Vilarrasa et al. (2011) indicated that the flow rate and particle travel velocity in the direction perpendicular to shear direction greatly increase because the connectivity of void space changes during shear, and fluid flow easily passes through the region of fracture with excellent connectivity and low energy loss. Matsuki et al. (2010) argued the difference in angles between shear macroscopic fluid flow directions and provided the formula of mechanical aperture (em) and hydraulic aperture (eh) with angle. Most experimental results showed that Darcy flow is infrequent in fracture with rough joint surfaces, especially under complex conditions. Zhang and Nemcik (2013) found that Darcy’s law applies only to the mated rough fracture, because the small aperture and high tortuosity of the flow pathway resist fluid flow severely and reduce flow velocity, and the behaviour of non-Darcy flow for all non-mated fractures is remarkable. In the experiments by Zhang and Nemcik (2013) the nonDarcy flow behaviour is more obvious in high-roughness fracture. Zhang et al. (2017) showed that the local Darcy and non-Darcy flow are concomitant in some conditions because of the inhomogeneous aperture and flow distribution, and the linear and non-linear seepage models need to be applied simultaneously. The external stress exerts a strong influence on flow behaviour, especially when the confining stress is increased. Chen et al. (2019) found that the non-Darcy behaviour of fluid flow increases with the reduction of the hydraulic pressure and the various distributions of the interconnected void areas. Some studies (Javadi et al., 2014; Rong et al., 2016) showed that the non-Darcy behaviour of fluid flow decrease with shear process because of the decreasing influence of roughness, but the other factors are unclear. For the boundedness of the cubic law above, several researchers aimed to modify the cubic law, and the most widely used approach is establishing the relation between em and eh by considering different conditions (Li et al., 2016). Olsson and Barton (2001) proposed the empirical engineering model for the hydromechanical coupling of rock joints by considering joint roughness coefficient (JRC) and mobilized JRC during shear. The model indicated that the relation between em and eh is nonlinear and separated by the peak shear displacement. Considering aperture as 3D distribution, Zimmerman and Bodvarsson (1996) established the relation between em and the standard deviation (σapert) of aperture, and then calculated eh with em, σapert, and contact ratio during shear. Xiong et al. (2011) specified that results are only approximate when the σapert/em ratio in Zimmerman and Bodvarsson (1996)’s model approaches 0.816. Xiong et al. (2011) added the standard deviation of local fracture slope and Reynolds number (reflecting the effect of inertia force from the side) from Zimmerman and Bodvarsson (1996)’s model to describe the relationship between em and eh. Xie et al. (2015) also used σapert to describe the relation between em and eh and explained the change in σapert. The distribution of the

absolute velocity and pressure drop under constant pressure gradient during shear is then obtained. The results of the abovementioned studies indicated that em is always larger than eh, which implies that the cubic law overstates the fracture hydraulic conductivity and can be modified by replacing eh with em. By contrast, Wang et al. (2015) modified the cubic law directly by considering local tortuosity and roughness. After indicating the validity of the modified cubic law, Wang et al. (2015) found that the arithmetic mean of effective errors for the studied fracture is reduced from 40.8% to 3.7% compared with local cubic law. In this study, we investigated the flow regime in a single fracture with joint surface destruction, gouge material, high hydraulic pressure, and anisotropic roughness, and analysed the difference between the cubic law and experimental results during shear process. Firstly, the coupled shear-flow tests were presented in artificial single fractures at a constant normal stress of 1.91 MPa and different high hydraulic pressures of 0.2, 0.4, 0.6 and 0.8 MPa with radial fluid flow. Then, 3D simulation models were established by using COMSOL Multiphysics 5.2a (COMSOL Multiphysics, 2016) with the Reynolds-averaged N–S k-ε turbulence model for initial and peak shear conditions. The radial flow cubic law was modified by establishing the relationship of mechanical and hydraulic aperture with hydraulic gradient. 2. Materials and methods 2.1. Specimens and fracture sample preparation To guarantee the repeatability of the tests and reduce costs, we built the test specimens with plaster, water, and retardant at a weight ratio of 1:0.25:0.005. Table 1 shows the physico-mechanical properties of the specimen. The upper and lower specimens are shown in Fig. 2 (Cao et al., 2018). All the specimens were 200 mm in diameter and 75 mm in height, and the surfaces of the specimens were regular dentate with 30° along the shear direction. As shown in Fig. 2(b), an injection hole 8 mm in diameter in the centre of the lower specimen was installed. Thus, the fluid flow was transported radially through the fracture between the upper and lower specimens. 2.2. Experimental test method The test apparatus TJXW-600 microprocessor control coupled shearflow test system was used in this paper. The side view of the apparatus is shown in Fig. 3 (Cao et al., 2018). The apparatus can provide the shear and normal loads of up to 600 kN by the hydraumatic servo oil source, and the maximum hydraulic pressure can reach 3 MPa. Normal and shear displacement were measured by using a cable displacement meter. The data of normal displacement, normal stress, shear displacement, shear stress, and aperture flow rate can be collected during shear. The initial setting parameter was stable, and the shear box sealing property was assured during test. As shown in Fig. 2(c), the boundary of the side and bottom of the upper and lower specimens was no-slip boundary, and the shear process was achieved by fixing the lower specimen and moving the upper specimen. For the tests, normal stress condition was constant normal load (CNL) boundary at 1.91 MPa. The boundary of the water inlet and Table 1 Physico-mechanical properties of the upper and lower specimens. Physico-mechanical properties Density Compressive strength Modulus of elasticity Poission's ratio Cohesion Internal friction angle

3

Unit

Value 3

g/cm MPa GPa MPa °

2.066 38.80 28.70 0.23 5.3 60

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

200 mm

200 mm

Water inlet

(a)

(b)

75 mm

Shear direction

Fracture 200 mm

(c) Fig. 2. Experimental specimens: (a) upper specimen, (b) lower specimen, (c) initial shear status of the upper and lower specimens.

outlet was at constant hydraulic pressure boundaries of 0.2, 0.4, 0.6 and 0.8. Normal stress started to increase at the beginning of the test, and the inlet hydraulic pressure initially increased when normal stress reached 0.6 MPa, when the hydraulic pressure increased to the set value steadily and the shear process was operated, and the tests terminated when shear displacement (δ) reached 32 mm.

(u· ) =

µT = Cµ

I + (µ + µ T )( u + ( u )T)] + F

(u· ) k =

(1) (2)

·(u ) = 0

·

µ+

µT k

k + Pk

µT

k2

+ C 1 Pk k

C2

2

k

(4) (5) (6)

where ρ is the fluid density, u is the flow velocity tensor, I is the unit tensor, μ is the viscous coefficient, μT is the turbulent viscous coefficient, F is the body force tensor, k is the turbulent kinetic energy, σk is the turbulent Prandtl number of turbulent kinetic energy, Pk is the turbulent energy generated by average velocity gradient, ε is the dissipation rate, σε is the turbulent Prandtl number of dissipation rate, and Cε1, Cε2, and Cμ are the calculated coefficients. As the default constant in COMSOL Multiphysics 5.2a, Cε1 = 1.44, Cε2 = 1.92, and Cμ = 0.09. The model range contained the interspace between the upper and lower specimens in initial and peak shear conditions. The boundary condition of water inlet was flow velocity boundary, and the values were calculated with the flow rate from the experimental results and the area of the simulation models. Water outlet was the pressure boundary, and the values were 0 in all simulation models. The joint surfaces were the wall boundary. Fig. 4 shows the schematic diagram of the distribution of aperture of simulation models, in which e0 is the initial mechanical aperture and Δun is the normal displacement of aperture due to shear stress. Considering that the contact areas between the lower and upper joint surfaces and e0 cannot be measured during

To investigate the distribution of streamlines and flow velocity before shear was broken, we established 3D seepage simulation models for initial and peak shear conditions with COMSOL Multiphysics 5.2a (COMSOL Multiphysics, 2016). Referring to Zhang et al. (2019a)’s study we simulated the fluid flow in small aperture with the k-ε turbulence model, COMSOL Multiphysics 5.2a provides the Reynoldsaveraged N–S k-ε modules and the governing equations are shown as follows (COMSOL Multiphysics, 2016):

·[

µ+

Pk = µT [ u: ( u + ( u )T )]

2.3. Numerical simulation method

(u· ) u =

·

(3) 4

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

3 2 6

4

5 7 8

1 10

9

11

1. Control system, 2. Hydraumatic servo oil source, 3. Hydraulic pressure system, 4. Normal loads, 5. Shear loads, 6. Displacement meter, 7. Upper shear box, 8. Lower shear box, 9. Water inlet, 10. Water outlet, 11. Flow rate measure system Fig. 3. Side view of TJXW-600 microprocessor control coupled shear-flow test system.

10 mm

Peak part

10 mm

Peak part Shear direction

Shear direction

Aperture Δ un

Aperture Valley part

Valley part

(a)

(b)

Fig. 4. Schematic diagram of aperture distribution of simulation models: (a) initial shear condition, (b) peak shear condition.

tests, we neglected the contact areas and assumed that e0 was distributed equally everywhere as shown in Fig. 4(a). By comparing the flow rates of initial shear condition models and experimental results, e0 can be back-calculated through the least square method and the value was 0.201 mm. For coupled shear-flow tests, em can be confirmed using the following equation (Esaki et al., 1999):

em = e0

un + un

calculation parameters of the initial and peak shear condition models are shown in Table 2. As shown in Fig. 5, the flow rate results of the simulation models were close to experiments. Thus, the value of e0 we calculated was usable. 3. Results versus cubic law

(7)

3.1. Cubic law

where un is the normal displacement of aperture due to normal stress. For Eq. (7) un was equal to 0 under CNL boundary (Xiong et al., 2011). Thus, the aperture of the peak shear condition model can be determined by adding e0 and Δun (0.481, 0.492, 0.474 and 0.475 respective) at the peak shear displacement (δp) (Fig. 4(b)). The

The radial flow cubic law assumes that the joint surfaces of fracture are smooth and paralleled, the fluid flow state is laminar and distributes axisymmetric, and can be described as follows (Zhang and Nemcik, 2013): 5

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

0.8

Table 2 COMSOL model settings for the flow in different conditions.

0.2 0.4 0.6 0.8

Mechanical aperture (mm)

Flow velocity of water inlet (m/s)

Initial

Peak shear

Initial

Peak shear

0.201

0.682 0.693 0.675 0.676

2.11 3.13 4.26 5.22

1.95 2.98 4.07 5.09

0.6 0.5

Δu (mm)

Hydraulic pressure (MPa)

0.2 MPa 0.4 MPa 0.6 MPa 0.8 MPa

0.7

0.4 0.3 0.2

Q=

2 g 12 ln

() ri ro

0.1

· h· e3 (8)

0.0

where Q is the flow rate in fracture, g is the gravitational acceleration, υ is the fluid kinematic viscosity, ro is the water inlet radius, ri is the joint surface outer radius, Δh is the difference in the pressure head of water inlet and outlet and e is the fracture aperture. When Eq. (8) was used to calculate the flow rate of a fracture, e = em was applied and the result was denoted as Qc. When Eq. (8) was used to identify eh, Q = Qe was obtained by applying e = eh.

20

25

30

35

As shown in Fig. 7(b), higher hydraulic pressure obviously contributed to the accuracy loss of the radial flow cubic law. The increment in Qc/Qe ratio was different when the hydraulic pressure increased by 0.2 MPa and the accuracy of the radial flow cubic law was more sensitive to the increase in lower hydraulic pressure. The radial flow cubic law is a linear equation simplified from the N–S equations after disregarding inertial forces and assumes that the fluid flow in the fracture is Darcy flow. In our study, as shown in Fig. 8, when we added the origin of coordinates, which represented that Qe = 0 when the hydraulic pressure equalled to 0, the relationship between pressure gradient and Qe was nonlinear and fitted with Forchheimer law as follows (Zhang and Nemcik, 2013):

(9)

P = aQe + bQe2

(10)

where P is the pressure gradient, and a and b are model coefficients. As shown in Fig. 8, the relationship between P and Qe could be described by non-linear Forchheimer law accurately. Thus, the nonDarcy flow behaviour was remarkable in this study, and the flow regime did not satisfy the assumption of the radial flow cubic law. The influence intension of inertial forces of fluid flow can be estimated by using the Reynolds number, which represents the specific value of inertial and viscous forces. Thus, a high Reynolds number would cause non-Darcy flow, and Reynolds number increases proportionally with flow velocity. Hence, the high hydraulic pressure increases the nonlinear behaviour of fluid flow and simultaneously reduces the accuracy of the radial flow cubic law, and the fluid flow with high hydraulic pressure does not obey Darcy’s law.

60

Flow rate (cm3·s-1)

15

3.3. Dependency of flow rate on hydraulic pressure

where the length unit in Eq. (9) is mm. As shown in Fig. 7(b), the Qc/Qe ratio was much larger than 1 during shear process with different hydraulic pressure, which indicated that the radial flow cubic law seriously deviated from the experiments and overstated the hydraulic conductivity of the fracture at the entire process of shear. The change in Qc/Qe ratio indicated that the accuracy of the radial flow cubic law reduced dramatically before δp. Then, the difference between the radial flow cubic law and experimental results stabilised because of the joint surface destruction. Simulation value (initial) Experimental value (initial) Simulation value (peak) Experimental value (peak)

3.4. Nonaxisymmetry of streamlines and flow velocity

40

The distribution of streamlines with 0.6 MPa hydraulic pressure of initial and peak shear condition models is shown in Fig. 9. The streamlines distribution with other hydraulic pressure is shown in Appendix A (Figs. A1 and A2). As shown in Figs. 9, A1 and A2, results showed that the distribution of streamlines was nonaxisymmetric with water inlet, which deviated from the assumption of the radial flow cubic law. In addition, the streamlines were tortuous when bypassed by the dentate aperture and close to groove through the water inlet. Given the uneven distribution of aperture, the stream line of the peak shear condition model was untidy and sparse. To investigate the absolute values of flow velocity distribution of different flow directions, we divided the joint surface into 11 flow directions, as shown in Fig. 10(a). Fig. 10(b) shows the results of

30

20

10

10

Fig. 6. Change in the normal displacement of aperture due to shear stress with shear displacement.

The change in Δun, measured flow rate (Qe) and Qc during shear tests are shown in Figs. 6 and 7. As shown in Figs. 6 and 7(a), Δun increased signally before δp and then changed slightly, Qe had a remarkable increment before δp and hydraulic pressure had a strong influence in Qe but did not change Δun. Combining the calculated e0, we can change Eq. (7) to the following equation to confirm em:

50

5

δ (mm)

3.2. Estimated and measured flow rates for shear tests

em = 0.201 + un

0

0.2

0.4

0.6

0.8

P (MPa) Fig. 5. Flow rates observed in the experiments and estimated with the COMSOL models. 6

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

1000

0.2 MPa 0.4 MPa 0.6 MPa 0.8 MPa

100

800

Qc/Qe ratio

Qe (cm3/s)

80

60

600

400

40

0.2 MPa 0.4 MPa 0.6 MPa 0.8 MPa

200

20 0

5

10

15

20

25

30

35

δ (mm)

0

0

5

10

15

20

25

30

35

δ (mm)

(b)

(a)

Fig. 7. Results of flow rate: (a) measured flow rate under different hydraulic pressure, (b) change in Qc/Qe ratio with shear displacement.

tortuosity ratio (τ), which is the actual seepage path relative to the radius of joint surface (Walsh and Brace, 1984), and the maximum inclination angle (θmax) of asperities in different flow directions. As shown in Fig. 10(b), τ distributed symmetrically with the 90° flow direction and the flow directions with the same τ had different θmax values. Both τ and θmax were strongly affected by the shape and density of roughness elements. The absolute values of flow velocity of the initial and peak shear condition models with 0.6 MPa hydraulic pressure are shown in Fig. 11. The flow velocity results of other hydraulic pressure are shown in Appendix A (Figs. A3 and A4). As shown in Figs. 11, A3 and A4, the flow velocity of different flow directions was diverse in the initial and peak shear condition models, and the values of flow velocity reduced with increased radial distance because of the energy loss and fluid flow expansion. In addition, the values of flow velocity reduced with increasing

τ at the same radius because the seepage path was extended with increase in τ and the energy loss was remarkable when fluid flow impacted joint surface and bypassed asperities. The flow velocity evolution curve became more fluctuant, and the frequency and amplitude of the curves were augmented with the increase in τ. The minimal value points of the curves appeared in the peak and valley parts in Fig. 4. In the flow directions with the same τ, large θmax advanced the phase of the flow velocity evolution curve. The results agreed with the study of Zhang et al. (2019a), where the increased roughness elements density reduced the hydraulic conductivity of fracture and the flow velocity, and the energy loss became drastic. The shape of roughness elements exerted strong influence in the velocity and regime of flow fluid. The waveform of the flow velocity evolution curves with the same flow direction were more undulating and complex in peak shear condition models resulting from the uneven distribution of aperture (Fig. 4(b)), and the flow velocity was reduced drastically in the void space with large aperture. Moreover, high hydraulic pressure increased the flow velocity of all flow directions and slightly influenced the distribution of streamlines. Thus, the inertia force of fluid flow increased in all sides. When the flow velocity was divided to the radial velocity (RV) and the velocity vertical to the radial distance (VVR), taking the 0°, 90° and 180° flow directions of the initial and peak shear condition models with 0.6 MPa for example, the absolute values of RV and VVR results are shown in Fig. 12. As shown in Fig. 12, the values of RV were larger than the values of VVR in both initial and peak shear condition models. For the 90° flow direction the values of VVR were dinky comparing with RV. Thus, the values of VVR reduced with decreasing roughness. Moreover, the value difference between RV and VVR was strongly influenced by τ and θmax.

Pressure gradient (MPa/m)

8

6

4

2

0

0

20

40

60

δ=1.84 mm δ=2.95 mm δ=4.91 mm δ=9.16 mm δ=13.20 mm δ=16.53 mm δ=20.36 mm δ=26.28 mm δ=32.00 mm

R2=0.9918 R2=0.9949 R2=0.9937 R2=0.9926 R2=0.9908 R2=0.9899 R2=0.9887 R2=0.9901 R2=0.9910

80

100

4. Discussion 4.1. Impact of joint surface destruction As shown in Fig. 7(b), the change in Qc/Qe ratio stabilized after δp, which indicated that the destruction of joint surfaces strongly influenced the fluid flow regime and the accuracy of the radial flow cubic law. The destruction of joint surfaces and the forms of gouge material in the peak and final shear conditions is shown in Fig. 13. As shown in Fig. 13, the joint surfaces were destructed completely after δp and the

Qe (cm3/s) Fig. 8. Regression analysis of pressure gradient as a function of measured flow rate using Forchheimer’s law during shear. 7

Journal of Hydrology 579 (2019) 124142

Flow velocity

Flow velocity

C. Cao, et al.

(a)

(b)

Fig. 9. Simulated results of streamlines distribution with 0.6 MPa hydraulic pressure: (a) initial shear condition model, (b) peak shear condition model.

20°

0°

1.40 60

40°

1.35

50

1.30

40

1.25

60°

200 mm

100°

1.20

30

τ

90°

θmax (°)

80°

1.15

20

θmax τ

1.10

10

120°

1.05

0

140° 160°

1.00 0

180°

20

40

60

80

100 120 140 160 180

Angle (°)

(b)

(a)

Fig. 10. Angle-dependent properties of the join surface: (a) Schematic diagram of different flow directions divided in specimen, (b) change in tortuosity ratio and maximum inclination angle in different flow directions. 1E+01 0° 20°

1E+00

40° 60° 80°

1E-01

90° 100° 120° 140°

1E-02

160° 180°

1E-03

0

20

40

60

80

Absolute value of velocity (m/s)

Absolute value of velocity (m/s)

1E+01

Radial distance (mm)

20° 40° 60°

1E-01

80° 90° 100°

1E-02

120° 140° 160°

1E-03

1E-04

100

0°

1E+00

180°

0

20

40

60

80

100

Radial distance (mm)

(a)

(b)

Fig. 11. Results of the absolute values of flow velocity in different flow directions with 0.6 MPa hydraulic pressure: (a) initial shear condition model, (b) peak shear condition model.

gouge material buried the regular dentate asperities completely in the final shear condition. In addition, the forms of the gouge material changed during the shear tests. At the peak shear condition, the gouge material was massive, and then the volume and distribution of the

gouge material were reduced and rearranged with shear destruction. The gouge material was ground to small granular and powdery in the final shear condition. The JRC of the fracture reduced immediately after δp because of the 8

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

1E+00

Absolute value of velocity (m/s)

Absolute value of velocity (m/s)

1E+01

1E-01

1E-02

RV 0° RV 90° RV 180°

1E-03

1E-04

0

20

VVR 0° VVR 90° VVR 180°

40

1E+01

1E+01

1E+00

1E+00

1E-01

1E-01

1E-02

1E-02

1E-03

1E-03 RV 0° RV 90° RV 180°

1E-04

60

80

1E-05

100

Radial distance (mm)

0

20

VVR 0° VVR 90° VVR 180°

40

1E-04

60

80

1E-05 100

Radial distance (mm)

(b)

(a)

Fig. 12. Results of the absolute values of radial velocity and the velocity vertical to the radial distance in 0°, 90°, and 180° flow directions with 0.6 MPa hydraulic pressure: (a) initial shear condition model, (b) peak shear condition model.

destruction of joint surfaces and the distribution of gouge material with different forms. The contact areas, τ and θmax of the upper and lower joint surfaces were changed dynamically at the same time. Fig. 14 shows the schematic of the distribution of the gouge material in the peak and final shear conditions. After δp, the seepage channels of fracture increased and became complex, and the τ and θmax of different flow directions became more diverse than before δp. The contact areas of fracture were decided by the forms of the gouge material, and the distribution of aperture was greater even when compared with the

200 mm

10 mm

200 mm

(a) 10 mm

Peak shear condition

Final shear condition

(b)

Fig. 13. Joint surface destruction of lower specimen and forms of gouge material during shear.

Fig. 14. Schematic diagram of aperture distribution, joint surface destruction and gouge material: (a) peak shear condition, (b) final shear condition.

9

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

0.095

0.115

0.2 MPa R2=0.9770 0.4 MPa R2=0.9494 0.6 MPa R2=0.9647 0.8 MPa R2=0.9916

0.090

0.110 0.105

0.085

0.100

eh (mm)

eh (mm)

0.2 MPa R2=0.9454 0.4 MPa R2=0.9602 0.6 MPa R2=0.9567 0.8 MPa R2=0.9574

0.080

0.095 0.090

0.075

0.085 0.080

0.070

0.075

0.065

0.2

0.3

0.4

0.5

0.6

0.72

0.7

0.76

0.80

0.84

0.88

em (mm)

em (mm)

(b)

(a)

Fig. 15. Regression analysis of mechanical and hydraulic aperture: (a) before peak shear displacement, (b) after peak shear displacement.

initial to peak shear conditions. Moreover, fluid flow could transport through the small granular and powdery gouge material, which may be treated as a porous medium. These results could be used to explain the reason for the steady trend of Qc/Qe ratio after δp and the non-Darcy flow behaviour decrease in the studies by Javadi et al. (2014) and Rong et al. (2016). Thus, the influence of surfaces destruction and gouge material in single-fracture seepage model during shear, especially the continuous destruction of joint surfaces and gouge material needs to be investigated.

the radial flow cubic law can be modified by combining Eqs. (8) and (11) and the relation among α, β and J. 4.3. Comparison to similar existing studies Many studies investigated the seepage model of a single fracture with complex conditions. Comparing to similar existing studies, we focused on the influence of high hydraulic pressure, joint surface destruction, gouge material and anisotropy of flow regime. Fractures with tens to hundreds of metre of water head is common in nature and engineering, because the limit of test apparatus in most coupled shearflow tests is reported with low J: 0–117 (Rong et al., 2016); 2–16.7 (Javadi et al., 2014), 0–120 (Yin et al., 2017), 0.5–10 (Li et al., 2008), and 0–21 (Olsson and Barton, 2001). Other experiments (Zhang and Nemcik, 2013; Zhou et al., 2015; Zhang et al., 2017; Chen et al., 2019) with high J were carried out with confining stress and concentrated on the nonlinear relationship between hydraulic gradient and flow rate by using Forchheimer or Lzbash laws, which cannot reflect the relationship between aperture and flow rate and were hard to apply to fracture network seepage models. In the present study, we investigated the influence of high hydraulic pressure in the accuracy of radial flow cubic law and the relationship between eh and em during shear process and then modified the radial flow cubic law with J. Single-fracture flow tests were mainly conducted under confining or shear stress. For the confining stress flow tests (Chen et al., 2019; Zhang et al., 2017; Zhang and Nemcik, 2013; Zhou et al., 2015; Rong et al., 2017), the aperture of fracture was only influenced by confining stress, the relative position and geometrical morphology of joint surfaces. The aperture did not change during tests. For the shear-flow tests (Esaki et al., 1999; Javadi et al., 2010; Olsson and Barton, 2001; Rong et al., 2016; Xiong et al., 2011; Yin et al., 2017) the fracture aperture and the relative position of joint surfaces changed with shear process. For most stress-flow tests, the fracture joints were unbroken, but in nature it is common for the destruction of fracture joints with complex stress condition and gouge material is widespread in some fracture. Thus, we analysed the influence of joint surface destruction and gouge material in the geometrical morphology of fracture, accuracy of radial flow and the relationship between em and eh. Although the analysis of joint surface destruction and gouge material in our study was primary, we thought it was potential for other similar studies.

4.2. Modification of the cubic law According to results above, the radial flow cubic law deviated from the experiments seriously and must be modified when applied to experiments or the natural fractures, particularly with high hydraulic pressure. For the tests in this paper, we modified the radial flow cubic law by replacing eh with em in two stages separated by δp. According to the experimental results the relation between em and eh with different hydraulic pressures is shown in Fig. 15 and can be described with the following equation:

eh =

+ em,

<

p;

eh = e

em ,

>

p

(11)

where α and β are the fitting coefficients. As shown in Fig. 15, em was much larger than eh because the effective seepage channels were reduced and the reduction of em was limited by the gouge material. In addition, eh had a linear relationship with em before δp and then translated to an exponential relationship resulting from the joint surface destruction and gouge material. Eq. (11) shows that eh > 0 when em = 0 before and after δp, which indicated that no matter how small the fracture aperture was, the fluid flow could transport though the fracture in the tests of this paper. Many connected seepage channels appeared in the closed fracture because the fracture with rough joint surfaces could not completely close. Durham and Bonner (1994) reported that fluid flow could transport though the closed fracture even when normal stress reaches 160 MPa. The change in α and β with different hydraulic gradients (J) is shown in Fig. 16, where J represents the theoretical value calculated by the corresponding hydraulic and the radius of fracture and is a dimensionless number. Before δp, α and β changed with J linearly, after δp, α and β had power function relationships with J. Thus, in this paper, 10

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

0.076

0.0146

a R2=0.9421 a=0.0793-2.2425E-5J b R2=0.9689 b=0.0132+1.775E-6J

1.52 0.0144

1.50

0.030

0.072

b

a 0.068

a

0.0142

0.070

1.48

a R2=0.9806 a=0.1134J -0.2320 b R2=0.9631 b=0.8896J 0.0788

0.028

0.0140

b

0.074

1.54

0.032

1.46 1.44

0.066

0.026

0.0138

1.42

0.064

200

300

400

500

600

700

1.40

0.024

0.0136

0.062 800

200

J

300

400

500

600

700

800

1.38

J

Fig. 16. Regression analysis of fitting coefficients of Eq. (11): (a) before peak shear displacement, (b) after peak shear displacement.

We also investigated the distribution of streamlines and flow velocity in different flow directions and found obvious anisotropy of flow regime caused by anisotropic joint surfaces, indicating that the seepage model for a settled fracture should be adjusted according to fluid flow direction. Although the anisotropy of flow regime and flow velocity was analysed with artificial fracture in this study, we thought it was common in natural fracture. However, few studies discussed about this issue, thus, further experiments or simulation works are suggested.

Given these limitations, relevant works are suggested to perform in single-fracture shear-flow studies with experimental and numerical methods that consider the influence of joint surface destruction, gouge material, anisotropy of flow regime and high hydraulic pressure. 5. Summary and conclusions This study aimed to investigate the fluid flow regime in single fracture during shear process influenced by high hydraulic pressure, joint surface destruction, gouge material and anisotropy of flow regime. Coupled shear flow tests were conducted with 0.2, 0.4, 0.6 and 0.8 hydraulic pressure, and 3D simulation models were established by using COMSOL Multiphysics 5.2a with Reynolds-averaged N–S k-ε turbulence model. The main summaries and conclusions are drawn as follows:

4.4. Limitation remarks Despite the meticulous care we have taken, the following limitations were not avoided in our study: (1) Coupled shear-flow tests: Similar to most shear-flow tests, the shear box was sealed and invisible in our tests. We cannot observe the destruction of joint surfaces, the distribution of gouge material with different forms, and the change in joint surface JRC, τ and θmax during shear process at any time. The change in these factors changed dynamically during shear tests and determined the establishment of seepage model. (2) Simulation model: The simulation models we used cannot consider the destruction of joint surfaces and the form change of the gouge material, especially the coupling between fluid flow and gouge material. Further numerical analysis is suggested to investigate the dynamic changes in the existing factors in the single-fracture seepage model. (3) Joint surfaces of fracture: Coupled shear-flow tests were conducted on artificial joint with regular dentate asperities, because the influence of joint surfaces destruction and gouge material were remarkable and transparently with regular dentate joint in our study, but the distribution of roughness and aperture of natural joint was spatial and arbitrarily. Thus, when the results of our study were applied to natural fracture, the geometrical morphological relationship between the different joint surfaces needs be established. According to the suggestion of Mirzaghorbanali et al. (2014), we can extend the results of this study to natural fracture by mapping the joint surface geometry with Fourier transformation or fractal methods.

(1) The accuracy of radial flow cubic law decreased dramatically before peak shear displacement, and then tended to be steady resulting from shear broken. High hydraulic pressure decreased the accuracy by increasing the influence of inertia force, and the accuracy was more sensitive to the increase in lower hydraulic pressure. (2) Joint surfaces destruction and gouge material exerted remarkable influence on the fracture geometrical morphology and flow regime. Because of joint surface destruction and gouge material the roughness, tortuosity ratio and maximum inclination angle of joint surfaces decreased dynamically during shear process, and the distribution of aperture and effective seepage channels became unpredictable. The difference between cubic law and experimental results and flow regime became steady because of the joint surface destruction and gouge material, especially the forms of gouge material. (3) The non-Darcy flow behaviour was obvious in our study because of the high hydraulic pressure and complex geometrical morphology of joint surfaces. The anisotropic behaviour of fluid flow was serious, and the streamlines and flow velocity were varied in different flow directions resulting from the anisotropic roughness of the joint. Both tortuosity ratio and maximum inclination angle had strong influence in the change in flow velocity, and the absolute values of velocity vertical to the radial distance reduced with

11

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

decreasing roughness. (4) The radial flow cubic law was modified pirical equation between mechanical and different hydraulic gradient. Mechanical showed a linear relationship before peak tended to be non-linear because of joint gouge material.

interests or personal relationships that could have appeared to influence the work reported in this paper.

by establishing the emhydraulic aperture with and hydraulic aperture shear displacement and surface destruction and

Acknowledgements This work was supported by National Natural Science Foundation of China (51922088, 51679197) and Project 2019JLM-57 supported by Joint Foundation of Shaanxi, China.

Declaration of Competing Interest The authors declare that they have no known competing financial Appendix A

Flow velocity

Flow velocity

The simulated results of streamlines distribution with 0.2, 0.4 and 0.8 MPa hydraulic pressure from simulation models are shown in Figs. A1 (initial shear condition) and A2 (peak shear condition). The change in the absolute flow velocity of different flow directions with 0.2, 0.4 and 0.8 MPa hydraulic pressure from simulation models are shown in Figs. A3 (initial shear condition) and A4 (peak shear condition).

(a)

Flow velocity

(b)

(c) Fig. A1. Simulated results of streamlines distribution at initial shear condition: (a) simulation model with 0.2 MPa hydraulic pressure, (b) simulation model with 0.4 MPa hydraulic pressure, (c) simulation model with 0.8 MPa hydraulic pressure.

12

Journal of Hydrology 579 (2019) 124142

Flow velocity

Flow velocity

C. Cao, et al.

(a)

Flow velocity

(b)

(c) Fig. A2. Simulated results of streamlines distribution at peak shear condition: (a) simulation model with 0.2 MPa hydraulic pressure, (b) simulation model with 0.4 MPa hydraulic pressure, (c) simulation model with 0.8 MPa hydraulic pressure.

13

Journal of Hydrology 579 (2019) 124142

C. Cao, et al. 1E+01 0°

1E+00

Absolute value of velocity (m/s)

Absolute value of velocity (m/s)

1E+01

20° 40° 60°

1E-01

80° 90° 100°

1E-02

120° 140° 160°

1E-03

1E-04

180°

0

20

40

60

80

0°

Radial distance (mm)

40° 60° 80°

1E-01

90° 100° 120° 140°

1E-02

160° 180°

1E-03

100

20°

1E+00

0

20

40

60

80

100

Radial distance (mm)

(a)

(b)

Absolute value of velocity (m/s)

1E+01 0° 20°

1E+00

40° 60° 80°

1E-01

90° 100° 120° 140°

1E-02

160° 180°

1E-03

0

20

40

60

80

100

Radial distance (mm)

(c)

Fig. A3. Simulated results of the absolute values of flow velocity in different flow directions at initial shear condition: (a) simulation model with 0.2 MPa hydraulic pressure, (b) simulation model with 0.4 MPa hydraulic pressure, (c) simulation model with 0.8 MPa hydraulic pressure.

14

Journal of Hydrology 579 (2019) 124142

C. Cao, et al.

1E+01

1E+00

0°

1E-01

40°

20° 60° 80°

1E-02

90° 1E-03

100° 120°

1E-04

140° 160°

1E-05 1E-06

180° 0

20

40

60

80

Absolute value of velocity (m/s)

Absolute value of velocity (m/s)

1E+01

100

0°

1E+00

20° 40°

1E-01

60° 80°

1E-02

90° 100° 120°

1E-03

140° 160°

1E-04 1E-05

Radial distance (mm)

180°

0

20

40

60

80

100

Radial distance (mm)

(b)

(a)

Absolute value of velocity (m/s)

1E+01 0°

1E+00

20° 40° 60°

1E-01

80° 90° 100°

1E-02

120° 140° 160°

1E-03

1E-04

180°

0

20

40

60

80

100

Radial distance (mm)

(c) Fig. A4. Simulated results of the absolute values of flow velocity in different flow directions at peak shear condition: (a) simulation model with 0.2 MPa hydraulic pressure, (b) simulation model with 0.4 MPa hydraulic pressure, (c) simulation model with 0.8 MPa hydraulic pressure.

Resour. Res. 50 (2), 1789–1804. https://doi.org/10.1002/2013WR014610. Li, B., Jiang, Y., Koyama, T., Jing, L., Tanabashi, Y., 2008. 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. 45 (3), 362–375. https://doi.org/10.1016/j.ijrmms.2007.06.004. Li, B., Liu, R., Jiang, Y., 2016. Influences of hydraulic gradient, surface roughness, intersecting angle, and scale effect on nonlinear flow behavior at single fracture intersections. J. Hydrol. 538, 440–453. https://doi.org/10.1016/j.jhydrol.2016.04. 053. Huang, N., Liu, R., Jiang, Y., Cheng, Y., Li, B., 2019. Shear-flow coupling characteristics of a three-dimensional discrete fracture network-fault model considering stress-induced aperture variations. J. Hydraul. Eng. 571, 416–424. https://doi.org/10.1016/j. jhydrol.2019.01.068. Mashayekhizadeh, V., Kharrat, R., Ghazanfari, M.H., Dejam, M., 2012. An experimental investigation of fracture tilt angle effects on frequency and stability of liquid bridges in fractured porous media. Petrol. Sci. Technol. 30 (8), 807–816. https://doi.org/10. 1080/10916466.2010.492370. Matsuki, K., Kimura, Y., Sakaguchi, K., Kizaki, K., Giwelli, A.A., 2010. Effect of shear displacement on the hydraulic conductivity of a fracture. Int. J. Rock. Mech. Min. 47 (3), 436–449. https://doi.org/10.1016/j.ijrmms.2009.10.002. Mirzaghorbanali, A., Nemcik, J., Aziz, N., 2014. Effects of cyclic loading on the shear behaviour of infilled rock joints under constant normal stiffness conditions. Rock Mech. Rock Eng. 47 (4), 1373–1391. https://doi.org/10.1007/s00603-013-0452-1. Olayiwola, S.O., Dejam, M., 2019. A comprehensive review on interaction of nanoparticles with low salinity water and surfactant for enhanced oil recovery in sandstone and carbonate reservoirs. Fuel 241, 1045–1057. https://doi.org/10.1016/j. fuel.2018.12.122. Olsson, R., Barton, N., 2001. An improved model for hydromechanical coupling during shearing of rock joints. Int. J. Rock. Mech. Min. 38 (3), 317–329. https://doi.org/10. 1016/S1365-1609(00)00079-4. Rong, G., Hou, D., Yang, J., Cheng, L., Zhou, C., 2017. Experimental study of flow

References Auradou, H., Drazer, G., Boschan, A., Hulin, J., Koplok, J., 2006. Flow channeling in a single fracture induced by shear displacement. Geothermics 35 (5–6), 576–588. https://doi.org/10.1016/j.geothermics.2006.11.004. Azarmanesh, M., Dejam, M., Azizian, P., Yesiloz, G., Mohamad, A.A., Sanati-Nezhad, A., 2019. Passive microinjection within high-throughput microfluidics for controlled actuation of droplets and cells. Sci. Rep. UK 9, 6723. https://doi.org/10.1038/ s41598-019-43056-2. Cao, C., Xu, Z., Chai, J., Qin, Y., Ran, T., 2018. Mechanical and hydraulic behaviors in a single fracture with asperities crushed during shear. Int. J. Geomech. 18 (11), 04018148. https://doi.org/10.1061/(ASCE)GM.1943-5622.0001277. Chen, Y., Lian, H., Liang, W., Yang, J., Nguyen, V.P., Bordas, S.P.A., 2019. The influence of fracture geometry variation on non-Darcy flow in fractures under confining stresses. Int. J. Rock. Mech. Min. 113, 59–71. https://doi.org/10.1016/j.ijrmms. 2018.11.017. COMSOL Multiphysics, 2016. 5.2a User’s Guide. Stockholm, Sweden. Dejam, M., Hassanzadeh, H., Chen, Z., 2018. Shear dispersion in a rough-walled fracture. SPE J. 23 (5), 1669–1688. https://doi.org/10.2118/189994-PA. Durham, W.B., Bonner, B.P., 1994. Self-propping and fluid flow in slightly offset joints at high effective pressures. J. Geophys. Res. Sol. Ea. 99 (B5), 9391–9399. https://doi. org/10.1029/94JB00242. Esaki, T., Du, S., Mitani, Y., Ikusada, K., Jing, L., 1999. Development of a shear-flow test apparatus and determination of coupled properties for a single rock joint. Int. J. Rock. Mech. Min. 36 (5), 641–650. https://doi.org/10.1016/S0148-9062(99)00044-3. Javadi, M., Sharifzadeh, M., Shahriar, K., 2010. A new geometrical model for non-linear fluid flow through rough fractures. J. Hydrol. 389 (1–2), 18–30. https://doi.org/10. 1016/j.jhydrol.2010.05.010. Javadi, M., Sharifzadeh, M., Shahriar, K., Mitani, Y., 2014. Critical Reynolds number for nonlinear flow through rough-walled fractures: the role of shear processes. Water

15

Journal of Hydrology 579 (2019) 124142

C. Cao, et al. characteristics in non-mated rock fractures considering 3D definition of fracture surfaces. Eng. Geol. 220, 152–163. https://doi.org/10.1016/j.enggeo.2017.02.005. Rong, G., Yang, J., Cheng, L., Zhou, C., 2016. Laboratory investigation of nonlinear flow characteristics in rough fractures during shear process. J. Hydrol. 541 (Part B), 1385–1394. https://doi.org/10.1016/j.jhydrol.2016.08.043. Tan, R., Chai, J., Cao, C., 2019. Experimental investigation of the permeability measurement of radial flow through a single rough fracture under shearing action. Adv. Civ. Eng. 2019, 671729. https://doi.org/10.1155/2019/6717295. Thompson, M.E., Brown, S.R., 1991. The effect of anisotropic surface roughness on flow and transport in fractures. J. Geophys. Res. 96 (B13), 21923–21932. https://doi.org/ 10.1029/91JB02252. Vilarrasa, V., Koyama, T., Neretnieks, I., Jing, L., 2011. Shear-induced flow channels in a single rock fracture and their effect on solute transport. Transp. Porous Med. 87 (2), 503–523. https://doi.org/10.1007/s11242-010-9698-1. Walsh, J.B., Brace, W.F., 1984. The effect of pressure on porosity and the transport properties of rock. J. Geophys. Res. 89 (B11), 9425–9431. https://doi.org/10.1029/ JB089iB11p09425. Wang, L., Cardenas, M.B., Slottke, D.T., Ketcham, R.A., Sharp Jr., J.M., 2015. Modification of the local cubic law of fracture flow for weak inertia, tortuosity, and roughness. Water Resour. Res. 51 (4), 2064–2080. https://doi.org/10.1002/ 2014WR015815. Wang, S.Y., Sun, L., Au, A.S.K., Yang, T.H., Tang, C.A., 2009. 2D-numerical analysis of hydraulic fracturing in heterogeneous geo-materials. Constr. Build. Mater. 23, 2196–2206. https://doi.org/10.1016/j.conbuildmat.2008.12.004. Xiao, W., Xia, C., Wang, W., Du, S., Deng, R., 2016. Contact algorithm for determining aperture evolution of rock fracture during shearing. Int. J. Geomech. 16 (3), 04015070. https://doi.org/10.1061/(ASCE)GM.1943-5622.0000581. Xie, L.Z., Gao, C., Ren, L., Li, B.C., 2015. Numerical investigation of geometrical and

hydraulic properties in a single rock fracture during shear displacement with the Navier-Stokes equations. Environ. Earth Sci. 73 (11), 7061–7074. https://doi.org/10. 1007/s12665-015-4256-3. Xiong, X., Li, B., Jiang, Y., Koyama, T., Zhang, C., 2011. Experimental and numerical study of the geometrical and hydraulic characteristics of a single rock fracture during shear. Int. J. Rock. Mech. Min. 48 (8), 1292–1302. https://doi.org/10.1016/j.ijrmms. 2011.09.009. Yin, Q., Ma, G., Jing, H., Wang, H., Su, H., Wang, Y., 2017. Hydraulic properties of 3D rough-walled fractures during shearing: an experimental study. J. Hydrol. 555, 169–184. https://doi.org/10.1016/j.jhydrol.2017.10.019. Zhang, Q., Luo, S., Ma, H., Wang, X., Qian, J., 2019a. Simulation on the water flow affected by the shape and density of roughness elements in a single rough fracture. J. Hydrol. 573, 456–468. https://doi.org/10.1016/j.jhydrol.2019.03.069. Zhang, W., Dao, B., Liu, Z., Zhou, C., 2017. A pore-scale numerical model for non-Darcy fluid flow through roughwalled fractures. Comput. Geotech. 87, 139–148. https:// doi.org/10.1016/j.compgeo.2017.02.011. Zhang, X., Chai, J., Qin, Y., Cao, J., Cao, C., 2019b. Experimental study on seepage and stress of single-fracture radiation flow. KSCE J. Civ. Eng. 23 (3), 1132–1140. https:// doi.org/10.1007/s12205-019-1519-7. Zhang, Z., Nemcik, J., 2013. Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures. J. Hydrol. 477, 139–151. https://doi.org/10.1016/j. jhydrol.2012.11.024. Zhou, J., Hu, S., Fang, S., Chen, Y., Zhou, C., 2015. Nonlinear flow behavior at low Reynolds numbers through rough-walled fractures subjected to normal compressive loading. Int. J. Rock. Mech. Min. 80, 202–218. https://doi.org/10.1016/j.ijrmms. 2015.09.027. Zimmerman, R.W., Bodvarsson, G.S., 1996. Hydraulic conductivity of rock fractures. Transport. Porous. Med. 23 (1), 1–30. https://doi.org/10.1007/BF00145263.

16