Empirical prediction of hydraulic aperture of 2D rough fractures: a systematic numerical study

Abstract

This study aims to propose an empirical prediction model of hydraulic aperture of 2D rough fractures through numerical simulations by considering the influences of fracture length, average mechanical aperture, minimum mechanical aperture, joint roughness coefficient (JRC) and hydraulic gradient. We generate 600 numerical models using successive random additions (SRA) algorithm and for each model, seven hydraulic gradients spanning from 2.5 × 10⁻⁷ to 1 are considered to fully cover both linear and nonlinear flow regimes. As a result, a total of 4200 fluid flow cases are simulated, which can provide sufficient data for the prediction of hydraulic aperture. The results show that as the ratio of average mechanical aperture to fracture length increases from 0.01 to 0.2, the hydraulic aperture increases following logarithm functions. As the hydraulic gradient increases from 2.5 × 10⁻⁷ to 1, the hydraulic aperture decreases following logarithm functions. When a relatively low hydraulic gradient (i.e., 5 × 10⁻⁷) is applied between the inlet and the outlet boundaries, the streamlines are of parallel distribution within the fractures. However, when a relatively large hydraulic gradient (i.e., 0.5) is applied between the inlet and the outlet boundaries, the streamlines are disturbed and a number of eddies are formed. The hydraulic aperture predicted using the proposed empirical functions agree well with the calculated results and is more reliable than those available in the preceding literature. In practice, the hydraulic aperture can be calculated as a first-order estimation using the proposed prediction model when the associated parameters are given.

Graphical abstract

Keywords

fluid flow / rough fracture surface / mechanical aperture / hydraulic aperture / predictive model

1 Introduction

An appropriate estimation of the hydraulic characteristics of fluid flow through fractured rock masses is of great importance in many practical applications such as oil/gas extraction and storage (Kim et al., 2007; Wang et al., 2015b; Wang et al., 2017; Ju et al., 2019), geothermal energy development (Hou et al., 2018; Kumari and Ranjith, 2019; Gong et al., 2020; Guo et al., 2021), remediation of contaminated groundwater (Zhao et al., 2016; Liu et al., 2021a), and nuclear waste disposal (Awual et al., 2020).

It has been confirmed by previous works that the hydraulic aperture (e_h) exerts a dominate control on hydraulic characteristics of fractured rock masses, which is potentially sensitive to mechanical aperture (e_m), surface roughness and hydraulic gradient (J) (Zimmerman and Bodvarsson, 1996; Xiong et al., 2013; Javadi et al., 2014; Zou et al., 2017; Dang et al., 2019; Liu et al., 2020; Tan et al., 2020). The mechanical aperture is defined as the arithmetic average distance between the opposite fracture walls (Barton et al., 1985; Renshaw, 1995; Zhao et al., 2017), which is independent of inertial effect, has been widely used as an important parameter to characterize the geometric properties of single fractures in previous studies. However, the distributions of e_m and the surface roughness of natural fractures are random and irregular. There are rarely smooth fractures in natural rock masses and the surface roughness gives rise to the increment of the flow paths and the decrement of the permeability/transmissivity (Li et al., 2021; Liu et al., 2022). Many methods have been used to characterize fracture surface roughness, such as root-mean-square of first derivative of asperity height, Hurst exponent, fractal dimension, and joint roughness coefficient (JRC) (Myers, 1962; Barton, 1974; Babadagli et al., 2015; Huang et al., 2019), where JRC is an important factor for estimating the fracture surface roughness. Many researchers found that the relationship between flow rate (Q) and J in rough-walled fractures is nonlinear deviating from the cubic law (Zimmerman and Bodvarsson, 1996; Cao et al., 2019; Sun et al., 2020; Tan et al., 2020). The e_h is typically derived by back-calculating the cubic law assuming a smooth parallel-plate fracture model (Snow, 1970; Tsang and Witherspoon, 1981; Barton et al., 1985; Nowamooz et al., 2009; Wang et al., 2015a; Yin et al., 2017).

The natural rough fracture as shown in Fig.1(a) usually features lots of contacts, which will reduce the equivalent permeability of the fracture. While maintaining equal permeability, rough fracture containing contacts can be equivalent to the fracture with smaller mechanical aperture without contacts (Liu et al., 2015; Li et al., 2019). As shown in Fig.1(b), the fracture length l' of equivalent model generally equals to the fracture length l of original rough fracture, while the mechanical aperture e' of equivalent model will be smaller than the mechanical aperture e of original rough fracture. At this point, the fracture can be divided into a series of 2D fracture profiles without contacts for permeability estimation. Besides, in practical engineering, it is difficult to measure the 3D void spaces inside fractures, but predicting the fluid flow behaviors through 2D outcrop profiles is an effective method (Bisdom et al., 2016; Yu et al., 2022). So, the simplified model as shown in Fig.1(b) has considered the influence of contacts and can be used for a fast calculation of permeability. A series of typical fluid flow tests and numerical simulations had been carried out to study the evaluation of e_h under different e_m, JRC, and J, and predicted models of e_h had been proposed (Yeo et al., 1998; Liu et al., 2016a; Huang et al., 2017; Sun et al., 2020; Tan et al., 2020; Zhang and Chai, 2020; Liu et al., 2021b). Ge et al. (2019) investigated the relationship between JRC and e_h by numerical simulation using the lattice Boltzmann method, and their study focuses on artificially created 2D fractures with random roughness following Gaussian distributions. Sun et al. (2020) reported a new method for predicting the e_h of rough rock fractures, by considering both the effects of the e_m and fracture morphology regardless of the influence of J. The predicted models of e_h reported in the previous studies are listed in Tab.1. It can be seen that some studies only established the relationship between JRC and the e_m (Patir and Cheng, 1978; Barton et al., 1985; Olsson and Barton, 2001; Scesi and Gattinoni, 2007; Chen et al., 2017), while some studies only revealed the relationship between e_m and J (Witherspoon et al., 1980; Liu, 2005; Rasouli and Hosseinian, 2011; Xie et al., 2015). The above studies mainly focused on the influences of JRC, e_m, and J on e_h of fractures. In real situation, the fluid flows through very small openings and continues its flow, whereas in a 2D channel there would be no fluid flow if the opening becomes zero. Therefore, the effects of the minimum (e_min) on hydraulic characteristics should be fully understood. Nevertheless, the influence of e_min and fracture length (l), on the evolution of e_h have rarely been considered.

This study systematically investigates the influences of e_m, l, e_min, JRC, and J on e_h and proposes a model for predicting e_h. First, we generated 600 2D rough rock fracture models (as shown in Tab.2) with the l = 100−400 mm, e_m = 4−20 mm, JRC = 1.82−11.44 and e_min = 0−20 mm, for each model spanning 7 orders of inlet pressure magnitude to fully cover both the linear and nonlinear flow regimes. For the models of l = 100 mm, the J varies from 10⁻⁶ to 10⁰; for the models of l = 200 mm, the J varies from 5 × 10⁻¹ to 5 × 10⁻¹; for the models of l = 300 mm, the J varies from 3.33 × 10⁻¹ to 3.33 × 10⁻¹; for the models of l = 400 mm, the J varies from 2.5 × 10⁻⁷ to 2.5 × 10⁻¹. Then, the numerical simulations were performed to investigate the hydraulic properties by solving the Navier-Stokes equations based on a finite volume method (FVM) code. Five parameters (e_m, e_min, JRC, l, J) are adopted to characterize the e_h. Finally, based on the back propagation neural network (BP-NN) algorithm, a training database containing 4200 fluid flow cases were established and the empirical prediction model of e_h was proposed.

2 Governing equations and numerical models

2.1 Governing equations

Our models are designed to represent the first-order roughness of fractures formed by primary asperities, which play dominant roles in the nonlinear fluid flow properties of fractures (Zou et al., 2015). The calculation formula proposed by Tse and Cruden (1979) has been suggested by International Society for Rock Mechanics (ISRM) (Brown, 1981) as a representative JRC measurement method, and has been widely accepted and applied by many studies (Brown, 1981; Yang et al., 2001; Rasouli and Hosseinian, 2011; Chen et al., 2015; Wang et al., 2023; Yin et al., 2023). Their equation is obtained based on the measurement results of the ten profiles presented by Barton and Choubey (1977), and the JRC values for the ten profiles are ranging from 0 to 20, covering the JRC range of this study. Thus, we consider this function as a representative. The expression of the equation is as follows:

(1)\mathit{J}\mathit{R}\mathit{C}=32.2+32.47\mathrm{l}\mathrm{g}{\mathit{Z}}_2,

where Z₂ is the root mean square of the first deviation of a profile (Myers, 1962; Tse and Cruden, 1979):

(2){\text{Z}}_{\text{2}} ={\left[\frac{\text{1}}{M}\sum {\left(\frac{z_{\mathrm{i}-1}-z_i}{x_{i-1}-x_i}\right)}^{\text{2}}\right]}^{1/2},

where x_i and z_i represent the x- and z-coordinate of the fracture surface profile i, respectively; M is the number of sampling points along the l (i.e., x-coordinate) of a fracture.

For incompressible and steady-state Newtonian fluid, fluid flow in fractures is governed by the Navier-Stokes equations, which are derived based on Newton’s second law (Foias et al., 2002):

(3)\rho [\mathrm{\partial }\mathit{u}/\mathrm{\partial }\mathit{t}+(\mathit{u}\cdot \mathrm{\nabla }) \mathit{u}]=- \mathrm{\nabla }\mathit{P}+\mathrm{\nabla }\cdot \mathit{T}+\rho \mathit{f},

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

For fractures that are conceptualized as two smooth parallel plates, the cubic law can be derived from the Navier-Stokes equations and applied in describing fluid flow in the fractures when the flow in the fractures is assumed to be laminar and in the steady-state regime (Tse and Cruden, 1979). In such cases, the nonlinear terms (u·∇) u can be deleted. The Navier-Stokes equations can be simplified to Eq. (4), which is the cubic law (Tsang and Witherspoon, 1981; Renshaw, 1995; Wang and Cardenas, 2014; Wang et al., 2015a; Yin et al., 2017):

(4)Q=-\frac{w\rho g{\text{e}}_{\text{h}}^3}{12\mu }J,

where Q is the flow rate, g is the gravitational acceleration, μ is the dynamic viscosity, w is the width of a fracture that equals to 1 m for 2D models. J is the hydraulic gradient, which can be expressed as

(5)J=-\Delta P/(\rho gl )，

where −ΔP is the difference in fluid pressure between inlet and outlet.

Based on the Eq. (4), the e_h in 2D fractures can be expressed as

(6)e_{\mathrm{h}}={(\frac{12Q\mu }{\rho gJ})}^{1/3}.

Equation (6) shows the standard method for calculating e_h, which has been widely accepted in the fields of rock mechanics and rock engineering (Witherspoon et al., 1980; Brown, 1981; Matsuki et al., 2010; Cardona et al., 2021). The Q can be calculated by solving cubic law or Navier-Stokes equations, according to requirements. For smooth parallel-plate models, the cubic law can be used for calculating Q when the flow velocity is sufficiently small. In this study, the fractur surface roughness and nonlinearity of fluid flow have been considered, so that the cubic law is not suitable and the Navier-Stokes equations are solved. For a rougher fracture or a larger hydraulic gradient that corresponds to a stronger nonlinearity, the Q is smaller, resulting in that the e_h is smaller by substituting Q into Eq. (3).

2.2 Generation of 2D numerical models

Natural fracture surfaces are rough and show self-affine characteristics (Brown and Scholz, 1985; Power and Tullis, 1991; Odling, 1994; Kulatilake et al., 2006; Ge et al., 2014). A large number of numerical methods have been developed to artificially generate fracture surface such as the successive random additions (SRA), the randomization of the Weierstrass-Mandelbrot function and the Fourier transformation (Wang et al., 2016). In the present study, the widely accepted SRA with a fractional Brown motion (fBm) is adopted, written as (Liu et al., 2004; Ye et al., 2015):

(7)

(8){\sigma }_{\mathit{r}\mathit{l}}{\mathit{=}\mathit{r}}^{\mathit{H}}{\sigma }_{\mathit{l}},

where < > is the mathematical expectation, r is a constant value, H is the Hurst coefficient that typically ranges between 0.45 and 0.85 for rock fracture surfaces (Odling, 1994; Schmittbuhl et al., 2008; Babadagli et al., 2015), Z is the asperity height, and σ is the square root of variance of Z. The asperity height increment [Z (x+rl_x, y+rl_y) – Z (x, y)] over distance \sqrt{l_x^2+l_y^2} follows Gaussian distribution (Liu and Molz, 1996). The Hurst exponent H is related to fractal dimension D, using the equation (Odling, 1994):

(9)D=E-H ,

where E is the Euclidean dimension of the embedding medium (3 for surfaces and 2 for surface profiles).

2.3 Numerical simulation processes

Using the SRA algorithm developed by (Liu et al., 2004), rough surfaces were generated using different random number seeds with H = 0.5. The 3D surfaces of fractures were divided into four groups and processed to the sizes of 100 mm × 100 mm, 100 mm × 200 mm, 100 mm × 300 mm and 100 mm × 400 mm, respectively.

A total of 20 3D models of the surface morphology were established and one example is shown in Fig.2. The 3D fracture models are composed of the upper rough surface and the lower rough surface. The color of the upper surface represents the height of the asperities.

The 2D models were generated using cutting planes that are perpendicular to the xy plane along z-direction with a sampling interval of 1 mm. Therefore, l of the 2D models are 100 mm, 200 mm, 300 mm, and 400 mm corresponding to the size of the 3D models. And JRC is in the range of 1.82−11.44, which is controlled by the asperity height of the 3D fracture models. During the formation of rough fractures, the protrusions on the rough walls of the fractures will be broken and worn under the action of shear force or chemical corrosion, resulting in different top and bottom wall morphologies of nonmatching fractures. The arithmetic mean of the two-walled roughness (JRC_ave) is often used to characterize the fracture morphology, and the influence of different morphologies can be analyzed by the minimum mechanical aperture (e_min), which is a parameter in this study (Rasouli and Hosseinian, 2011; Sun et al., 2020). In addition, expressing the JRC values of the upper and lower surfaces with average JRC values can help make one JRC value corresponding to one single fracture, which will help simplify the expression on the roughness of fractures and can be easily obtained through mapping outcrops of rock masses. Thus, the JRC of each 2D model is the average value of the JRC of the upper surface and that of the lower surface in this study. For each JRC, e_m was set in the range of 4−20 mm with an interval of 4 mm to investigate the effect of e_m on hydraulic characteristics. The minimum e_min is 0 mm and the maximum e_min is 20 mm. −ΔP varies from 10⁻³ to 10³, spanning 7 orders of magnitude to fully cover both linear and nonlinear flow regimes. One 2D model is shown in Fig.3. The upper and lower surfaces of the fracture are modeled as impermeable boundaries and assumed to be nondeformable. The fluid flows through the 2D model from the inlet at the left boundary to the outlet at the right boundary. The inlet and outlet were modeled as pressure-inlet and pressure-outlet boundaries, in which the inlet pressure and outlet pressure are P₁ and P₂, respectively. Each model was meshed in the mechanical aperture direction (perpendicular to the flow direction) using 20 layers, which is the same number of layers used by Xiong et al. (2013) and is verified to be sufficiently validate during the fluid flow simulation process (Xiong et al., 2011; Liu et al., 2017). An enlarged view of the meshes performed using ANSYS ICEM is shown in Fig.3.

The fractures were plotted in AutoCAD and exported as SAT files. The SAT files were imported into ANSYS ICEM for meshing. The quadrilateral meshes with a maximum side length of 0.2 mm were adopted to mesh the fractures. Thus, there are more than 20 layers along the mechanical aperture direction. The meshed fractures models were saved as MESH files and then imported into ANSYS FLUENT for calculating the Q through each model. The numerical simulations were performed to characterize nonlinear flow in rock fractures by solving the Navier-Stokes equations based on an FVM code, which has been verified in previous studies (Li et al., 2016a; Li et al., 2016b; Liu et al., 2016b, 2016c). No-slip boundary conditions were assigned to the two fracture walls. The density and viscosity of water are 998.2 kg/m³ and 0.001003 Pa·s, respectively, when the temperature is assumed to be 25°C. We have also checked the effect of number of iterations and found that the calculated results are stable after 1500 iterations. So, the number of iterations was determined as 1500.

When the hydraulic pressures are applied at the boundaries, Q can be calculated by solving Eq. (3). And by substituting Q into Eq. (6), the value of e_h can be obtained.

3 Results and discussion

3.1 Validity of proposed numerical simulation code

To verify the validity of the proposed numerical simulation code, 10 parallel plate models with l = 100 mm and e_m increases from 1 mm to 10 mm are established. A sufficiently small hydraulic pressure drop ( −ΔP/Δl) of 0.1 Pa/m is applied on the fracture models to guarantee linear flow between inlet and outlet. By solving the Navier-Stokes equations, Q through the model can be calculated. The theoretical values of Q are derived according to Eq. (4).

The numerically calculated results agree well with the theoretical results predicted by Eq. (4) as shown in Fig.4. The numerically calculated results are slightly smaller than the theoretical results when e_m ranges from 7 mm to 10 mm. This is because when e_m increases, the plane areas of the opening space imported into ANSYS ICEM for meshing also increase, the grid density will decrease relatively. The computational accuracy of the FLUENT software is precisely affected by the mesh density. Relative errors between numerical results calculated by proposed numerical method and theoretical results predicted by Eq. (4) are less than 0.5% for all cases, so the deviations are acceptable. These results efficiently verify the validity of the proposed numerical simulation code.

3.2 Influences of J, l, e_m, e_min, and JRC on hydraulic characteristics

For the models with e_m = 4 mm, e_m = 8 mm, e_m = 12 mm, e_m = 16 mm, and e_m = 20 mm, the variations in Q/J versus J are shown in Fig.5−Fig.6. The variations in Q/J exhibit three-stage characteristics with increasing J: a linear regime (weak inertial regime), a transition region, and a nonlinear region (strong inertial regime). J is used rather than Reynolds number because J is a macroscopic dimensionless parameter that generally has a known value in practices during hydraulic pump tests with prescribed hydraulic pressures (Li et al., 2016a). When e_m = 4 mm and l = 100 mm, the linear regime of the variation of Q/J exists in J = 10⁻⁶−10⁻⁵. When e_m = 4 mm, l = 200 mm, 300 mm, 400 mm, the linear regime of the variation of Q/J exists in J = 10⁻⁶−10⁻⁴. Therefore, the larger l, the larger J corresponding to fluid flow in the linear regime when e_m = 4 mm. When e_m = 8 mm, the linear regime of the variation of Q/J are much smaller than those of e_m = 4 mm. When the e_m is equal to 12 mm, 16 mm and 20 mm respectively, the Q/J has only an obvious nonlinear region, even when J = 10⁻⁷. When e_m is relatively small (i.e., 4 mm and 8 mm), Q/J is significantly influenced by e_min and JRC. When e_m is relatively large (i.e., 12 mm, 16 mm, 20 mm), the effect of the e_min and the JRC on the Q/J−J relationships is not obvious as shown in Fig.7 and Fig.6. This is because when e_m is sufficiently small, the flow rate is significantly influenced by the viscous force, which is predominantly controlled by e_min and JRC. When e_m is relatively large, the inertial force is much larger than the viscous force, and e_m, l, and J play the leading role in the evolution of Q/J in these cases.

3.3 Streamline distributions

To visualize the flow paths, a number of particles are injected at the inlet and the streamlines are recorded according to the particle variables as shown in Fig.10. According to Eq. (5), J can be calculated when −ΔP is given when l = 200 mm. When −ΔP is sufficiently small (i.e., −ΔP = 0.001 Pa and J = 5 × 10⁻⁷), the Q is relatively small and the effect of inertial forces can be neglected, resulting in that the fluid flow is in the linear regime. When −ΔP is large (i.e., − ΔP = 1000 Pa and J = 0.5), the fluid flow converts to the nonlinear flow regime and the effect of inertial forces cannot be negligible with respect to viscous forces. Here, J = 5 × 10⁻⁷ and J = 0.5 are chosen to represent the linear and nonlinear flow regimes, and the corresponding streamline distributions and velocity distributions are presented in Fig.10 and Fig.11, respectively. In the linear flow regime, the particles smoothly flow through the void spaces formed by the tortuous lower and upper surfaces. Since the viscous force is much larger than the inertial force, no eddies are formed, and the direction of the flow velocity vector is uniformly distributed from the inlet direction to the outlet direction (i.e., Fig.10(a), 10(c), 10(e), Fig.11(a), 11(c), 11(e)). In the nonlinear flow regime, the inertial force cannot be negligible with respect to the viscous force. Many eddies are located at different locations with different sizes and shapes. These eddies give rise to energy losses, decreasing the transmissivity/permeability of fractures. When JRC is smaller (i.e., JRC = 5.42), the eddies exist in the place where local aperture changes significantly (i.e., the partial enlarged view of Fig.10(b)) and there are almost no eddies in the place where local aperture does not change robustly. Whereas when JRC is larger (i.e., JRC = 10.74), the eddies are distributed within the total aperture fields, due to the influences of local aperture variations and rough surfaces of lower and upper walls (i.e., the partial enlarged view of Fig.10(f)). Therefore, the energy losses more significantly with a larger JRC, resulting in smaller transmissivity/permeability. Meanwhile, obvious backflow phenomenon can be observed at the eddies as shown in the partially enlarged view of Fig.11. The backflow in the nonlinear flow is caused by the inertial effect, which has been verified to decrease the permeability of fractured rock masses to a large extent.

3.4 Influences of J, l, e_m, and e_min on e_h

It is found that JRC has little effect on fluid flow characteristics within the mechanical aperture range of 4 mm to 20 mm. When the aperture is relatively large, nonlinear flow is mainly controlled by the aperture rather than JRC, which is consistent with the results in previous studies (Zou et al., 2015; Liu et al., 2016a). Thus, we choose the models with JRC randomly distributed from 1.82 to 5.71. In these cases, the e_h values are mainly controlled by l, e_m, and J. Fig.12 shows the variations in e_h with e_m/l varying from 0.01 to 0.2 and J varying from 2.5 × 10⁻⁷−10⁰. In all cases with different l, the variations of e_h with e_m/l follow the logarithm functions. For a smaller l (i.e., 100 mm), the fitting curve of e_h is closer to straight.

Fig.13 shows the effect of J on e_h under different e_min and l. When e_min varies from 10⁻⁷ to 10⁻¹ and l varies from 100 mm to 400 mm, and e_h decreases following the logarithm function with the increment of J. The values of e_h are close to e_min when J = 10⁻⁷, but e_h decreases with the increment of J. As shown in Fig.13(a), when J = 10⁻⁶, e_min and e_h are 19.15 mm and 16.95 mm, respectively, and when J = 10⁻⁴, e_h is 9.73 mm, which is exactly twice as small as e_min. This is because when J is sufficiently small, the inertial force is far smaller than the viscous force in the entire flow field and the determination of flow path is significantly influenced by e_min. Increment of J would give rise to the number and volume of eddies and backflows which narrow down the effective flow paths and consequently lead to the decrease of e_h.

3.5 Empirical prediction model of e_h

The e_h is significantly correlated with JRC, l, e_m, e_min, and J. Based on the database of 4200 fluid flow cases, an empirical prediction model of e_h is proposed based on the five parameters (l, e_m, JRC, e_min, and J). The database was imported into the MATLAB software, and the artificial neural network (ANN) method was applied to calculate the input and output of neurons in each layer. ANNs can promise models to account for implicit relationships between variables because their topology structure is similar to multilayer perceptrons (McCulloch and Pitts, 1943). As a traditional ANN only contains a forward-propagation stage, the back propagation neural network (BP-NN) is designed to reduce fitting errors by adding a back-propagation stage to adjust weights and thresholds online (Rumelhart et al., 1986). A three-layer structure, which includes input layer, hidden layer and output layer, is applied to present the information transmission, as shown in Fig.14.

In Fig.14, i nodes in the input layer, j nodes in the hidden layer and m nodes in the output layer are set. x, h, and y denote input-layer input, hidden-layer output and output-layer output, respectively. W₁ denotes weights between input and hidden layers, and W₂ denotes weights between hidden and output layers. B₁ and B₂ stand for hidden-layer and output-layer thresholds, respectively. W_1^{\prime }, W_2^{\prime }, B_1^{\prime }, and B_2^{\prime } denote corresponding adjusted weights and thresholds, and σu, σv, σb, and σd are the corresponding corrections, respectively. The process of BP-NN is illustrated as follows.

1) According to x_i, u_ij, and ε_j, the information flow from the input layer to the hidden layer is shown as

(10)\mathit{H}=f( {\mathit{W}}_1\mathit{X}+{\mathit{B}}_1),

where X and H denote the input and hidden vectors respectively; f() denotes the transfer function of the hidden layer, which always goes with logsig, tansig and the purelin function (Lawrence, 1993). Logsig function was adopted in this paper.

2) According to H, W₂, and B₂, mapping the information flow between the hidden layer and output layer is expressed as

(11)\mathit{Y}=\varphi ( {\mathit{W}}_2\mathit{X}+{\mathit{B}}_2),

where Y denotes output vector; \varphi () denotes the transfer function of the output layer, and a linear function was applied in our experimental setup.

3) We define the term \overline{\mathit{Y}} as the expected value, and then the modeling error reads as

(12)E(\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r})=\mathrm{s}\mathrm{u}\mathrm{m}(\mathit{Y}-\overline{\mathit{Y}}{)}^2.

4) The gradient descent method for finding optimal weights and thresholds to minimize the modeling error is presented as follows:

(13)\{\begin{array}{l}\mathit{\sigma }\mathit{u}=\kappa {\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }{\mathit{W}}_1}},\mathit{\sigma }\mathit{v}=\kappa {\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }{\mathit{W}}_2}}\\ \mathit{\sigma }\mathit{b}=\kappa {\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }{\mathit{B}}_1}},\mathit{\sigma }\mathit{d}=\kappa {\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }{\mathit{B}}_2}}\end{array},

where \kappa is the learning rate, representing the step length of the process involved in approximation for the optimal value.

5) Steps 1)−4) are executed cyclically until the modeling error reaches the setting threshold or the iteration number reaches the setting value. The inputs are the six parameters (l, e_m, JRC, J, e_min, e_h) from the database of 4200 fluid flow cases. The functional relationship between each parameter and the e_h is fitted. The BP-NN-based equation can be obtained as follows:

(14){\mathit{e}}_{\mathrm{h}}={\mathit{W}}_2\times \mathrm{l}\mathrm{g}\mathrm{s}\mathrm{i}\mathrm{g}( {\mathit{W}}_1\times \mathit{X}+{\mathit{B}}_1)+{\mathit{B}}_2,

where sig(x) = 1/(1 + exp(−x)); W₁, W₂, B₁, and B₂ are all coefficients in matrix form, and X is a matrix related to the five parameters, as follows:

(15)\mathit{X}=\left(\begin{array}{ccc}l& e_{\mathrm{m}}& JRC\end{array}\begin{array}{cc}& J\end{array}\begin{array}{cc}& e_{\mathrm{m}\mathrm{i}\mathrm{n}}\end{array}\right),

(16){\mathit{W}}_1=\left(\begin{array}{ccccc}1.12& 0.34& \text{0.11}& -\text{7.87}& \text{0.67}\\ \text{4.01}& \text{0.44}& \text{0.05}& -\text{3.9}& \text{0.79}\\ -\text{0.1}& -\text{0.34}& \text{0.05}& \text{4.1}& \text{0.93}\\ \text{0.66}& \text{0.3}& \text{0.09}& -\text{22}& \text{0.48}\\ -\text{0.09}& \text{0.05}& \text{0.02}& \text{2.39}& \text{0.81}\end{array}\right),

(17){\mathit{W}}_2=\left(\begin{array}{ccccc}-\text{6.84}& \text{1.68}& -\text{1.53}& \text{14.45}& \text{4.54}\end{array}\right),

(18){\mathit{B}}_1={\left(\begin{array}{ccccc}-\text{9.98}& -\text{7.94}& \text{3.48}& -\text{24.76}& \text{4.02}\end{array}\right)}^{\text{T}},

(19){\mathit{B}}_2=- \mathrm{3.85.}

The predictive models of e_h deduced in previous literature are presented in Tab.3. Barton et al. (1985) and Scesi and Gattinoni (2007) presented the prediction models of e_h by considering the effects of e_m and JRC. Rasouli and Hosseinian (2011) presented the models of e_h by considering the effects of e_m, e_min, and JRC. The models employed in their research are all 2D rough fractures based on the Barton standard curves, and the fluid flow in their models is assumed to be incompressible, laminar, isothermal and in the steady-state regime of a viscous Newtonian fluid, which is the same condition with this study. Therefore, the parameters used in the numerical simulation calculation of this study can be substituted into the previous prediction models for calculation. The results of the predicted model in the present study and the results of numerical simulation are presented in Fig.15. The e_h calculated by the prediction model proposed by Barton et al. (1985) is generally smaller than that of numerical simulation, in which the correlation coefficient R² between results calculated by Barton’s model and simulated values is 0.35. This is because the units of e_m in their study is microns and the fluid flow behaviors are governed to a large extent by surface characteristics such as JRC. The e_h calculated by the prediction model proposed by Scesi and Gattinoni (2007) are generally distributed on the both sides of the numerical simulation results of e_h, but the R² between results calculated by Scesi᾽s model and simulated values is less than 0.28, and the fitting degree is not as good as it between the results of the prediction model proposed in this study and the simulation results. Considerable errors arise from both the consideration of only two parameters (e_m and JRC) and a small number of models used in their investigations. The e_h calculated by the prediction model proposed by Rasouli and Hosseinian (2011) is generally larger than that of the numerical simulation, and R² between results calculated by Rasouli᾽s model and simulated values is 0.02. The reason for the error is that the prediction equation considers only geometric parameters to characterize e_h, ignoring the influence of J. The R² between the predicted results by the prediction model in this study and the numerical simulation results is 0.84. Taken together, these results suggest that compared with the other predictive methods, the e_h predicted by the proposed model agrees better with the simulation values, which correspondingly verifies the validity of the predictive models of e_h. Two main reasons account for the accuracy of the predictive models reported in this study. First, the five parameters (e_m, e_min, l, JRC, and J) can essentially determine the geometric characteristics and hydraulic characteristics of a fracture (Zimmerman and Bodvarsson, 1996). Second, BP-NN, as a comprehensive approach to addressing both regression and classification problems and has numerous applications in many scientific fields (Schmidhuber, 2015; Wu et al., 2021), is suitable to solve the high-dimension relationship between e_h.

To verify the applicability of the model proposed in this study, we take the model data of Sun et al. (2020) into calculation, and substituted the same parameters as their numerical model into the models proposed by this study and previous studies for calculation. The e_min ranges from 0.15 mm to 1.13 mm. The detailed parameter settings are listed in Tab.4, where e_h1 denotes the calculation results of e_h in Sun’s study, e_h2 denotes the e_h predicted by the empirical model proposed in this study, and e_h3,e_h4, and e_h5 correspond to the e_h predicted by the models proposed by Barton et al. (1985), Scesi and Gattinoni (2007), and Rasouli and Hosseinian (2011), respectively. Finally, we compared all of the predicted results with Sun’s calculated results and found that the empirical model proposed by our study still has a good accuracy (see Fig.16). This indicates that the empirical model proposed in the present study is reliable for apertures with varying values.

4 Conclusions

In this study, the Navier–Stokes equations were solved to characterize fluid flow in 600 rough rock models generated by successive random additions (SRA) algorithm. Based on the FLUENT numerical simulations and the artificial neural network (ANN) algorithm, a training database containing 4200 fluid flow cases was established. Average mechanical aperture (e_m), minimum mechanical aperture (e_min), joint roughness coefficient (JRC), fracture length (l) and hydraulic gradient (J) were adopted to train the empirical prediction model of the hydraulic aperture (e_h). The distribution of the streamlines and the evolution of flow rate (Q) were investigated.

When e_m is relatively small (i.e., 4 mm and 8 mm), the increment of l could increase the range of the linear regime of fluid flow. The hydraulic characteristics of fractures with relatively small e_m are significantly influenced by e_min and JRC. When e_m is relatively large (i.e., 12 mm, 16 mm, and 20 mm), the fluid flow shows obvious nonlinearity, and the effect of the e_min and the JRC on the Q/J−J relationships is not obvious. When e_m = 20 mm, the corresponding ranges of JRC and e_min are 0−11.26 and 16−20 mm respectively, the Q/J−Q curves in this range are almost coincident. This is because when e_m is sufficiently small, the flow rate is significantly influenced by the viscous force, which is predominantly controlled by e_min and JRC. When e_m is relatively large, the inertial force is much larger than the viscous force, e_m, l, and J play the leading role in the evolution of Q/J in these cases. The nonlinearly decreasing relationship between e_h and J follows logarithm functions, and e_h is close to e_min at a sufficiently small J, but e_h decreases downwards at relatively large J, i.e., when J = 10⁻⁶, e_min and e_h are 19.15 mm and 16.95 mm, respectively, and when J = 10⁻⁴, e_h is 9.73 mm, which is exactly twice as small as e_min. This is because when J is sufficiently small, the inertial force is far smaller than the viscous force in the entire flow field and the determination of flow path is significantly influenced by e_min. Increment of J would give rise to the number and volume of eddies and backflows, narrow down the effective flow paths, enhance the nonlinearity of fluid flow, and then consequently lead to the decrease of e_h.

An empirical model for predicting the e_h of rough rock fractures was proposed by considering the effects of e_m, e_min, JRC, l, and J. Five parameters (e_m, e_min, JRC, l, J) are adopted to characterized the hydraulic characteristics and e_h. The equation was compared with three existing empirical equations using the parameters from the database in this study, showing that the proposed equation has better and robust predictive performance.

The predictive model reported in this study could be easily extended for more complex cases, such as 3D rough fractures and fracture networks, and can even be extended to handle many engineering geology problems related to hydraulic attributes of the fractured rock. Noted that the fractures utilized in this study are simple lines cut from 3D models, the hydraulic properties of which might be different from those of natural cases. We will focus on 3D rough fractures in future works.

References

[1] Amadei B, Illangasekare T (1992). Analytical solutions for steady and transient flow in non-homogeneous and anisotropic rock joints.Int J Rock Mech Min Sci Geomech Abstr, 29(6): 561–572

[2] Awual M R, Yaita T, Kobayashi T, Shiwaku H, Suzuki S (2020). Improving cesium removal to clean-up the contaminated water using modified conjugate material.J Environ Chem Eng, 8(2): 103684

[3] Babadagli T, Ren X, Develi K (2015). Effects of fractal surface roughness and lithology on single and multiphase flow in a single fracture: an experimental investigation.Int J Multiph Flow, 68: 40–58

[4] Barton N (1974). Review of a new shear strength criterion for rock joints.Eng Geol, 11(11): 287–332

[5] Barton N, Bandis S, Bakhtar K (1985). Strength, deformation and conductivity coupling of rock joints.Int J Rock Mech Min Sci Geomech Abstr, 22(3): 121–140

[6] Barton N, Choubey V (1977). The shear strength of rock joints in theory and practice.Rock Mech, 10(1–2): 1–54

[7] Bisdom K, Bertotti G, Nick H M (2016). The impact of in-situ stress and outcrop-based fracture geometry on hydraulic aperture and upscaled permeability in fractured reservoirs.Tectonophysics, 690: 63–75

[8] Brown E (1981). Rock characterization, testing & monitoring: ISRM suggested methods.Biospectroscopy, 18(6): 109

[9] Brown S R, Scholz C H (1985). Broad bandwidth study of the topography of natural rock surfaces.J Geophys Res, 90(B14): 12575–12582

[10] Cao C, Xu Z, Chai J, Li Y (2019). Radial fluid flow regime in a single fracture under high hydraulic pressure during shear process.J Hydrol (Amst), 579: 124142

[11] Cardona A, Finkbeiner T, Santamarina J C (2021). Natural rock fractures: from aperture to fluid flow.Rock Mech Rock Eng, 54(11): 5827–5844

[12] Chen Y D, Liang W G, Lian H J, Yang J F, Nguyen V P (2017). Experimental study on the effect of fracture geometric characteristics on the permeability in deformable rough-walled fractures.Int J Rock Mech Min Sci, 98: 121–140

[13] Chen Y F, Zhou J Q, Hu S H, Hu R, Zhou C B (2015). Evaluation of forchheimer equation coefficients for non-darcy flow in deformable rough-walled fractures.J Hydrol (Amst), 529: 993–1006

[14] Cruz P, Quadros E, Fo D, Marrano A (1982). Evaluation of opening and hydraulic conductivity of rock discontinuities. The 23rd U.S Symposium on Rock Mechanics (USRMS), Berkeley: California

[15] Dang W, Wu W H K, Konietzky H, Qian J (2019). Effect of shear-induced aperture evolution on fluid flow in rock fractures.Comput Geotech, 114: 103152

[16] Foias C, Manley O, Rosa R, Temam R (2002). Navier-Stokes equations and turbulence.Phys Today, 55(10): 54–56

[17] Ge Y, Kulatilake P H S W, Tang H, Xiong C (2014). Investigation of natural rock joint roughness.Comput Geotech, 55: 290–305

[18] Ge Z, Yuan T, Li Y (2019). Numerical study on the mechanism of fluid flow through single rough fractures with different JRC. Physica, Mechanica & Astronomica, 49(1): 1–10

[19] Gong F C, Guo T K, Sun W, Li Z M, Yang B, Chen Y M, Qu Z Q (2020). Evaluation of geothermal energy extraction in Enhanced Geothermal System (EGS) with multiple fracturing horizontal wells (MFHW).Renew Energy, 151: 1339–1351

[20] Guo T K, Zhang Y L, He J Y, Gong F C, Chen M, Liu X Q (2021). Research on geothermal development model of abandoned high temperature oil reservoir in north China oilfield.Renew Energy, 177: 1–12

[21] Hakami E (1995). Aperture distribution of rock fractures. Dissertation for Doctoral Degree. Stockholm: Royal Institute of Technology

[22] Hou J C, Cao M C, Liu P K (2018). Development and utilization of geothermal energy in China: current practices and future strategies.Renew Energy, 125: 401–412

[23] Huang N, Liu R C, Jiang Y J (2017). Numerical study of the geometrical and hydraulic characteristics of 3D self-affine rough fractures during shear.J Nat Gas Sci Eng, 45: 127–142

[24] Huang N, Liu R C, Jiang Y J, Cheng Y, Li B (2019). Shear-flow coupling characteristics of a three-dimensional discrete fracture network-fault model considering stress-induced aperture variations.J Hydrol (Amst), 571: 416–424

[25] 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 Resour Res, 50(2): 1789–1804

[26] Ju Y, Dong J B, Gao F, Wang J G (2019). Evaluation of water permeability of rough fractures based on a self-affine fractal model and optimized segmentation algorithm.Adv Water Resour, 129: 99–111

[27] Kim J, Cho W, Chung I M, Heo J H (2007). On the stochastic simulation procedure of estimating critical hydraulic gradient for gas storage in unlined rock caverns.Geosci J, 11(3): 249–258

[28] Kulatilake P H S W, Balasingam P, Park J, Morgan R (2006). Natural rock joint roughness quantification through fractal techniques.Geotech Geol Eng, 24(5): 1181–1202

[29] Kumari W G P, Ranjith P G (2019). Sustainable development of enhanced geothermal systems based on geotechnical research – a review.Earth Sci Rev, 199: 102955

[30] LawrenceJ (1993). Introduction to Neural Networks. California Scientific Software

[31] Li B, Jiang Y (2013). Quantitative estimation of fluid flow mechanism in rock fracture taking into account the influences of JRC and Reynolds number. J Min Mater Process Inst Jpn, 129 (7): 479–484 (in Japanese)

[32] Li B, Li Y, Zhao Z, Liu R (2019). A mechanical-hydraulic-solute transport model for rough-walled rock fractures subjected to shear under constant normal stiffness conditions.J Hydrol (Amst), 579: 124153

[33] Li B, Liu R, Jiang Y (2016a). Influences of hydraulic gradient, surface roughness, intersecting angle, and scale effect on nonlinear flow behavior at single fracture intersections.J Hydrol (Amst), 538: 440–453

[34] Li J L, Li X H, Zhang B, Sui B, Wang P C, Zhang M (2021). Effect of lower surface roughness on nonlinear hydraulic properties of fractures.Geofluids, 2021: 6612378

[35] Li X, Feng Z J, Han G, Elsworth D, Marone C, Saffer D, Cheon D S (2016b). Breakdown pressure and fracture surface morphology of hydraulic fracturing in shale with H₂O, CO₂ and N₂.Geomechanics Geophys Geo-Energy Geo-Resour, 2(2): 63–76

[36] Liu E (2005). Effects of fracture aperture and roughness on hydraulic and mechanical properties of rocks: implication of seismic characterization of fractured reservoirs.J Geophys Eng, 2(1): 38–47

[37] Liu H H, Bodvarsson G S, Lu S L, Molz F J (2004). A corrected and generalized successive random additions algorithm for simulating fractional levy motions.Math Geol, 36(3): 361–378

[38] Liu H H, Molz F J (1996). Discrimination of fractional Brownian movement and fractional Gaussian noise structures in permeability and related property distributions with range analyses.Water Resour Res, 32(8): 2601–2605

[39] Liu J W, Wei K H, Xu S W, Cui J, Ma J, Xiao X L, Xi B D, He X S (2021a). Surfactant-enhanced remediation of oil-contaminated soil and groundwater: a review.Sci Total Environ, 756: 144142

[40] Liu R C, He M, Huang N, Jiang Y J, Yu L Y (2020). Three-dimensional double-rough-walled modeling of fluid flow through self-affine shear fractures.J Rock Mech Geotech Eng, 12(1): 41–49

[41] Liu R C, Jiang Y J, Li B (2016c). Effects of intersection and dead-end of fractures on nonlinear flow and particle transport in rock fracture networks.Geosci J, 20(3): 415–426

[42] Liu R C, Jing H W, Li X Z, Yin Q, Xu Z G, He M, Taylor T (2021b). An experimental study on fractal pore size distribution and hydro-mechanical properties of granites after high temperature treatment.Fractals, 29(4): 1–13

[43] Liu R C, Li B, Jiang Y J (2016a). Critical hydraulic gradient for nonlinear flow through rock fracture networks: the roles of aperture, surface roughness, and number of intersections.Adv Water Resour, 88: 53–65

[44] Liu R C, Wang Y S, Li B, Jing H W, Li S C, Yang H Q (2022). Linear and nonlinear fluid flow responses of connected fractures subject to shearing under constant normal load and constant normal stiffness boundary conditions.Comput Geotech, 141: 104517

[45] Liu R, Jiang Y, Li B, Wang X (2015). A fractal model for characterizing fluid flow in fractured rock masses based on randomly distributed rock fracture networks.Comput Geotech, 65: 45–55

[46] Liu R, Li B, Jiang Y (2016b). Critical hydraulic gradient for nonlinear flow through rock fracture networks: the roles of aperture, surface roughness, and number of intersections.Adv Water Resour, 88(Feb): 53–65

[47] Liu R, Yu L, Jiang Y (2017). Quantitative estimates of normalized transmissivity and the onset of nonlinear fluid flow through rough rock fractures.Rock Mech Rock Eng, 50: 1063–1071

[48] Lomize G (1951). Flow in Fractured Rock. Moscow: Gosemergoizdat, 127–129 (in Russian)

[49] Louis C, Maini Y (1969). Determination of in-situ hydraulic parameters in jointed rock.International Society of Rock Mechanics Proceedings, 1: 1–19

[50] Matsuki K, Kimura Y, Sakaguchi K, Kizaki A, Giwelli A (2010). Effect of shear displacement on the hydraulic conductivity of a fracture.Int J Rock Mech Min Sci, 47(3): 436e449

[51] McCulloch W S, Pitts W (1943). A logical calculus of the ideas immanent in nervous activity.Bull Math Biophys, 5(4): 115–133

[52] Myers N O (1962). Characterization of surface roughness.Wear, 5(3): 182–189

[53] Nowamooz A, Radilla G, Fourar M (2009). Non-Darcian two-phase flow in a transparent replica of a rough-walled rock fracture.Water Resour Res, 45(7): W07406

[54] Odling N E (1994). Natural fracture profiles, fractal dimension and joint roughness coefficients.Rock Mech Rock Eng, 27(3): 135–153

[55] Olsson R, Barton N (2001). An improved model for hydromechanical coupling during shearing of rock joints.Int J Rock Mech Min Sci, 38(3): 317–329

[56] Patir N, Cheng H S (1978). An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication.Int J Rock Mech Min Sci, 100(1): 12–17

[57] Power W L, Tullis T E (1991). Euclidean and fractal models for the description of rock surface roughness.Int J Rock Mech Min Sci Geomech Abstr, 28(6): A344

[58] Rasouli V, Hosseinian A (2011). Correlations developed for estimation of hydraulic parameters of rough fractures through the simulation of JRC flow channels.Rock Mech Rock Eng, 44(4): 447–461

[59] Renshaw C E (1995). On the relationship between mechanical and hydraulic apertures in rough-walled fractures.J Geophys Res, 100(B12): 24629–24636

[60] Rumelhart D E, Hinton G E, Williams R J (1986). Learning representations by back propagating errors.Nature, 323(6088): 533–536

[61] Scesi L, Gattinoni P (2007). Roughness control on hydraulic conductivity in fractured rocks.Hydrogeol J, 15(2): 201–211

[62] Schmidhuber J (2015). Deep learning in neural networks: an overview.Neural Netw, 61: 85–117

[63] Schmittbuhl J, Steyer A, Jouniaux L, Toussaint R (2008). Fracture morphology and viscous transport.Int J Rock Mech Min Sci, 45(3): 422–430

[64] Snow D T (1970). The frequency and apertures of fractures in rock.Int J Rock Mech Min Sci Geomech Abstr, 7(1): 23–40

[65] Sun Z H, Wang L Q, Zhou J Q, Wang C S (2020). A new method for determining the hydraulic aperture of rough rock fractures using the support vector regression.Eng Geol, 271: 105618

[66] Tan J, Rong G, Zhan H B, He R H, Sha S, Li B W (2020). An innovative method to evaluate hydraulic conductivity of a single rock fracture based on geometric characteristics.Rock Mech Rock Eng, 53(10): 4767–4786

[67] TsangY W, Witherspoon P A (1981). Hydromechanical behavior of a deformable rock fracture subject to normal stress. J Geophys Res, 86(B10 NB10): 9287–9298 10.1029/JB086iB10p09287

[68] Tse R, Cruden D M (1979). Estimating joint roughness coefficients.Int J Rock Mech Min Sci Geomech Abstr, 16(5): 303–307

[69] Waite M E, Ge S, Spetzler H (1999). A new conceptual model for fluid flow in discrete fractures: an experimental and numerical study.J Geophys Res, 104: 13049–13059

[70] Walsh J (1981). Effect of pore pressure and confining pressure on fracture permeability.Int J Rock Mech Min Sci Geomech Abstr, 18(5): 429–435

[71] Wang L C, Cardenas M B (2014). Non-Fickian transport through two-dimensional rough fractures: assessment and prediction.Water Resour Res, 50(2): 871–884

[72] Wang L C, Cardenas M B, Slottke D T, Ketcham R A, Sharp J M Jr (2015a). Modification of the Local Cubic Law of fracture flow for weak inertia, tortuosity, and roughness.Water Resour Res, 51(4): 2064–2080

[73] Wang L, Tian Y, Yu X Y, Wang C, Yao B W, Wang S H, Winterfeld P H, Wang X, Yang Z Z, Wang Y H, Cui J Y, Wu Y S (2017). Advances in improved/enhanced oil recovery technologies for tight and shale reservoirs.Fuel, 210: 425–445

[74] Wang M, Chen Y F, Ma G W, Zhou J Q, Zhou C B (2016). Influence of surface roughness on nonlinear flow behaviors in 3D self-affine rough fractures: Lattice Boltzmann simulations.Adv Water Resour, 96: 373–388

[75] Wang Y K, Zhang Z Y, Liu X Q, Xue K S (2023). Relative permeability of two-phase fluid flow through rough fractures: the roles of fracture roughness and confining pressure.Adv Water Resour, 175: 104426

[76] Wang Z C, Li S C, Qiao L P, Zhang Q S (2015b). Finite element analysis of the hydro-mechanical behavior of an underground crude oil storage facility in granite subject to cyclic loading during operation.Int J Rock Mech Min Sci, 73: 70–81

[77] Witherspoon P A, Wang J S Y, Iwai K, Gale J E (1980). Validity of cubic law for fluid-flow in a deformable rock fracture.Water Resour Res, 16(6): 1016–1024

[78] Wu Z H, Pan S R, Chen F W, Long G D, Zhang C Q, Yu P S (2021). A comprehensive survey on graph neural networks.IEEE Trans Neural Netw Learn Syst, 32(1): 4–24

[79] Xie L Z, Gao C, Ren L, Li C B (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

[80] 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 Sci, 48: 1292–1302

[81] Xiong X, Li B, Jiang Y, Koyama T, Zhang C (2013). Experimental and numerical study of the geometrical and hydraulic characteristics of a single rock fracture during shear.Seismic Safety Evaluation of Concrete Dams, 22: 513–538

[82] Yang Z Y, Lo S C, Di C C (2001). Reassessing the joint roughness coefficient (JRC) estimation using Z₂.Rock Mech Rock Eng, 34(3): 243–251

[83] Ye Z, Liu H H, Jiang Q, Zhou C (2015). Two-phase flow properties of a horizontal fracture: the effect of aperture distribution.Adv Water Resour, 76(feb): 43–54

[84] Yeo I W, de Freitas M H, Zimmerman R W (1998). Effect of shear displacement on the aperture and permeability of a rock fracture.Int J Rock Mech Min Sci, 35(8): 1051–1070

[85] Yin Q, Ma G W, Jing H W, Wang H D, Su H J, Wang Y C, Liu R C (2017). Hydraulic properties of 3D rough-walled fractures during shearing: an experimental study.J Hydrol (Amst), 555: 169–184

[86] Yin Q, Nie X X, Wu J Y, Wang Q, Bian K Q, Jing H W (2023). Experimental study on unloading induced shear performances of 3D saw-tooth rock fractures.Int J Min Sci Technol, 33(4): 463–479

[87] Yu L, Zhang J, Liu R C, Li S C, Liu D, Wang X L (2022). Estimation of the representative elementary volume of three-dimensional fracture networks based on permeability and trace map analysis: a case study.Eng Geol, 309: 106848

[88] Zhang Y, Chai J R (2020). Effect of surface morphology on fluid flow in rough fractures: a review.J Nat Gas Sci Eng, 79: 103343

[89] Zhao X, Liu W, Cai Z Q, Han B, Qian T W, Zhao D Y (2016). An overview of preparation and applications of stabilized zero-valent iron nanoparticles for soil and groundwater remediation.Water Res, 100: 245–266

[90] Zhao Y L, Zhang L Y, Wang W J, Tang J Z, Lin H, Wan W (2017). Transient pulse test and morphological analysis of single rock fractures.Int J Rock Mech Min Sci, 91: 139–154

[91] Zimmerman R W, Bodvarsson G S (1996). Hydraulic conductivity of rock fractures.Transp Porous Media, 23(1): 1–30

[92] Zou L C, Jing L, Cvetkovic V (2015). Roughness decomposition and nonlinear fluid flow in a single rock fracture.Int J Rock Mech Min Sci, 75: 102–118

[93] Zou L C, Jing L, Cvetkovic V (2017). Shear-enhanced nonlinear flow in rough-walled rock fractures.Int J Rock Mech Min Sci, 97: 33–45

RIGHTS & PERMISSIONS

Higher Education Press