# Discrepancy between Forward and Reverse Seepage Characteristics in a Single Rough Fracture

Authors: Yulong Niu; Yuan Wang; Zhiyu Sun; Jinghua Li; Xin Xiang; Zhikui Wang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Under the uniform seepage boundary condition, the forward and reverse seepage flow rates will vary widely. In response to this phenomenon, this paper studies the mechanism of the difference in seepage characteristics between the forward and reverse directions from the fracture profile’s roughness, Darcy and non-Darcy seepage, and distribution of eddy currents. First of all, wavelet transform was used to decompose and reconstruct the single rough fracture for 8 times, and then 9 fracture models with various roughness degrees were constructed. Then, Finite Volume Method was used to simulate the seepage in the original and approximate models of the Reynolds number varying from 0.001 to 1000. The results show that the nonlinear seepage is necessary for the difference between the forward and reverse seepage characteristics of the single rough fracture. The mechanism of the discrepancy between forward and reverse nonlinear seepages is the diverse shapes and distribution of the eddy current generated by the forward and reverse seepage. The secondary roughness provides space for the existence and development of the eddy current. The secondary roughness is the geometric basis of the discrepancy between the forward and reverse seepage characteristics.


---

## Body

## 1. Introduction

In fractured rock masses, fractures act as the dominant seepage channels controlling the entire rock mass’s hydraulic properties compared to the matrix with low permeability [1]. The study of hydraulic properties of rock fractures involves the transport of pollutants [2], geological storage of CO2 [3], geothermal and petroleum resources development and utilization [4, 5], prediction of seepage flow rate [6], and many other engineering fields. The single rough fracture is the basic unit of the rock fracture network. Thoroughly studying the fluid flow in the single rough fracture is the basis of understanding the fluid flow and solute transport in fractured rock masses [7, 8].Both the porous medium and the fractured rock mass have the feature of seepage anisotropy [1], and so does the single rough fracture. For example, the studies of Thompson and Brown [9], Grasselli et al. [10], and Cao et al. [11] show that the anisotropy of the fracture surface roughness is the main cause of the seepage anisotropy (X-axis and Y-axis). Yeo et al. [12], Auradou et al. [13], and Huang et al. [14] studied the effect of fracture space variation on seepage anisotropy. However, besides the phenomenon that the seepage characteristics of the X-axis and the Y-axis are different, the seepage characteristics in the forward and reverse directions of the X-axis (the +X-direction and the –X-direction) are also different. For example, Grasselli et al. [10] performed a single rough fracture seepage experiment in the +X-direction and –X-direction, respectively, with Reynolds numbers varying from 5 to 150. The experimental results show that the equivalent hydraulic aperture decreases as the Reynolds number increases. Moreover, the equivalent hydraulic aperture in the +X-direction is greater than that of –X-direction under the same Reynolds number. Boutt et al. [15] studied the transport of colloids in +X-direction and –X-direction in single rough fractures through numerical simulation, which confirmed the discrepancy between the forward and reverse seepages in single rough fractures. Xie et al. [16] found that the flow rate of +X-direction and –X-direction was different under shear displacement through numerical simulations related to the local eddy current. Cardenas et al. [17] studied eddy currents’ sensitivity to the fluid flow in the +X-direction and –X-direction and the effect of eddy current on solute transport. They also pointed out that the anisotropy of 3D fractures and the directional transport of 2D fractures are fracture fluid mechanics studies’ weaknesses. The laboratory experiments and numerical simulations show that the seepage characteristics in the forward and reverse directions are different in a single rough fracture, and it is related to the eddy current in the fracture. Nevertheless, how the eddy current affects the forward and reverse seepage characteristics and how the fracture’s roughness affects eddy currents need to be further studied.Patton [18], Lee et al. [19], Xie et al. [20], and other authors hold that the roughness of rough fractures exists as a multiple-scale and shows a characteristic of self-similarity. The International Society for Rock Mechanics (ISRM) divides the roughness of rock surfaces (joint or fracture) into large-scale fluctuations (also known as primary roughness) and small-scale fluctuations (also known as secondary roughness) [21]. Jing et al. [22], Kana et al. [23], and Yang et al. [24] studied the role of primary roughness and secondary roughness in the shear process of fractures. Within small displacement, the shear process is controlled by secondary roughness, and, within large displacement, the shear process is controlled by primary roughness. The control effect of roughness on the fluid flow in fractures is widely recognized. But the influence of the primary and secondary roughness on the seepage of the fracture is still in the stage of exploration; only Sharifzadeh [25], Zou et al. [26], Zou et al. [27], Wang et al. [28], Dou et al. [7], and other few scholars carried out related researches. Among them, Zou et al. [26] and Zou et al. [27] found that high-frequency secondary roughness is the main reason (except for Reynolds number) for the dynamic evolution of the eddy current, and the higher the high-frequency roughness, the more obvious and complex the eddy current region. Wang et al. [28] show that the primary roughness mainly controls the flow path and large-scale pressure distribution, while the secondary roughness controls the nonlinear seepage characteristics of the local-scale fluid flow. Meanwhile, as the pressure gradient increases, the secondary roughness increases the complexity of the local velocity distribution by creating and expanding eddy currents and recirculation regions. Dou et al. [7] found that the secondary roughness significantly enhanced the nonlinear flow (i.e., eddies and tortuous streamlines) and the non-Fickian transport. It can be observed in the existing research results that secondary roughness (high-frequency roughness) has a decisive effect on the occurrence and development of nonlinear seepage and eddy current regions. But whether the secondary roughness (high-frequency roughness) is related to the discrepancy between the forward and reverse seepage characteristics still needs further exploration.In summary, the discrepancy between the forward and reverse seepage characteristics of the single rough fracture remains to be discussed and studied. In order to study the above problems, the wavelet transform is used to decompose and reconstruct the single rough fracture profile generated by Brazil splitting test. The fracture profiles at different scales (different decomposition levels) are obtained, and the corresponding models are built. Then FVM (Finite Volume Method) was used to simulate seepage in forward and reverse directions under different Reynolds numbers. The mechanism of the discrepancy between forward and reverse seepage characteristics of the single rough fracture is analyzed from the aspects of equivalent hydraulic aperture, non-Darcy coefficient, and local seepage field. The correlation between inconsistent roughness in the forward and reverse directions and the inconsistent seepage characteristics is also discussed.


## 2. Fracture Profiles at Different Scales

For multiscale analysis of 2D rough fracture profiles, the Fourier transform [29] and the wavelet transform [30] are the most frequently used tools in signal processing disciplines. The geometry of the natural fracture profile can be viewed as a nonstationary signal. In the Fourier transform, there are some limitations to the analysis of nonstationary signals because only the stationary signal can be processed. Meanwhile, the wavelet transform can better handle the nonstationary signal, and the wavelet transform method is advantageous to analyze the multiscale effect of surface roughness on fluid flow through rock fractures [26–28, 31]. Therefore, the wavelet transform is used to build the profile of the single rough fracture at different scales.


### 2.1. Building Single Rough Fracture Profiles

In this paper, the rough fracture surface of coarse-grained granite is used as the original surface to create the fracture models at different scales through the wavelet transform. The size of the standard cylindrical rock sample is 100 mm × 50 mm. A binocular 3D laser scanning system was used to scan the Brazilian splitting test’s fracture surface. The binocular 3D laser scanning system’s accuracy is ±20μm, and the scanning interval is set at 0.1 mm on X-axis and 0.2 mm on Y-axis to obtain surface elevation data of the sample. We can get the digital fractured surface the same as the entire fracture surface through the points scanned. The vertical section of the middle of the specimen was taken as the sample’s typical fracture model’s original profile. The entire process of constructing an accurate rough fracture profile is shown in Figure 1.Figure 1
Generating process of the single fracture profile of coarse-grained granite: (a) specimen after Brazilian splitting test; (b) fracture surface formed by splitting; (c) 3D laser scanning; (d) fracture surface generated by scanning; (e) single rough fracture profile.
(a)(b)(c)(d)(e)


### 2.2. Basic Theory of the Wavelet Transform

The single rough fracture profile has characteristics of multiple scales and is similar to signals. Therefore, the 2D single rough fracture profile can be regarded as a signal functionft to be processed, where f represents the vertical height of fluctuations and t represents the horizontal distance along the fracture. When applying the wavelet transform, equation (1) is used.(1)Wfa,b=f,ψa,b=1a∫−∞+∞ftψt−ba¯dt,(2)ψta,b=1aψt−ba,a,b∈R;a≠0,where ψt is the mother wavelet function; ψa,bt is a wavelet sequence function; a is the scale factor; b is the translation factor. Therefore, the wavelet transform can be seen as a convolution of the signal function ft and a wavelet sequence function ψa,bt with the bandpass filter. A detailed discussion of the wavelet transform can be found in Zou et al.’s work [27].Before performing the wavelet transform, the appropriate mother wavelet function needs to be chosen first. Which wavelet to be chosen as the mother wavelet mainly depends on the application requirements [32]. Mirzaeian [33] pointed out that high-order Daubechies wavelets are more suitable for processing surface roughness than other wavelets in wavelet libraries. Sharifzadeh [25], Zou et al. [26], Zou et al. [27], and Wang et al. [28] used Db8 wavelets in the Daubechies wavelet series to process 2D rough single fracture profiles or 3D single rough fracture surfaces. Therefore, the Db8 wavelet was used as the mother wavelet function in the wavelet transform.


### 2.3. Multiscale Decomposition of Single Rough Fracture Profiles

In this section, the wavelet transform and reconstruction will be performed on the single rough fracture in Figure1(e) as the original fracture model S. The profile of the single rough fracture is defined as fx,x∈0,L, where L represents the length of the profile along the X-axis. In this paper, L = 100 mm. F is the height of the measured profile at the position of x. For the fracture profile, the Mallat algorithm is used for eight-level decomposition and reconstruction [34]. The approximate models are labeled as A1–A8, and the detailed models are labeled as D1–D8. The specific decomposition process labeled can be found in Zou et al.’s work [27].Figure2 shows the approximate models Ai and the detailed models Di of each decomposition level. The approximate models Ai represents the geometry that plays a dominant role in the fracture profile with a low frequency, large wavelength, and large amplitude. The detailed models Di represents small-scale ups and downs with a high frequency, small wavelength, and small amplitude compared with the approximate profile. As the level of decomposition increases, it can be seen that the wavelength and amplitude of Di gradually increase, varying from 0.05 mm in the first level to more than 0.5 mm in the eighth level. So, it gradually shows a large-scale characteristic (low frequency and large wavelength), which is similar to the entire profile. At the same time, as the level of decomposition increases, the approximate profile (blue line) gradually loses the original profile information compared with the entire profile, leading to the larger wavelength and smaller amplitude, and it slowly turns smooth and parallel.Figure 2
The approximate profiles Ai and detailed profiles Di under different decomposition levels: (a) Level 1. (b) Level 2. (c) Level 3. (d) Level 4. (e) Level 5. (f) Level 6. (g) Level 7. (h) Level 8.
(a)(b)(c)(d)(e)(f)(g)(h)According to Figure2, the approximate profile A4 shape remains the same as the original profile S. The shape of the approximate profile A5 is different from the original profile S. The fluctuation difference Δ (the highest vertical height minus the lowest vertical height) under different decomposition levels is counted, as shown in Figure 3. The inflection point of the curve of fluctuation difference Δ and decomposition levels is at level four. The value of fluctuation difference Δ before and after the fourth level changes significantly. Therefore, in this paper, the approximate profile A4 is taken as the primary roughness (wavy fluctuation structure). The original profile S minus the approximate profile A4 (D1 + D2 + D3 + D4) is considered as the secondary roughness (fine roughness structure).Figure 3
The variation of fluctuations difference Δ with decomposition level.The wavelet transform is used to decompose and reconstruct the rough fracture profiles, separating the fractured profiles into the generally approximate profiles and the locally detailed profiles in a unique way. It establishes the one-to-one relationship between locally detailed profiles, approximate profiles, and the original profile, making the roughness of the fracture profile unique and accurate at different scales [7, 26–28]. It also provides a new solution for studying the geometric characteristics of the fracture roughness at different scales. At the same time, it also lays the geometric foundation for analyzing the effect of the detailed profile of rough fractures significantly when the secondary roughness affects the discrepancy between the forward and reverse hydraulic characteristics.


## 2.1. Building Single Rough Fracture Profiles

In this paper, the rough fracture surface of coarse-grained granite is used as the original surface to create the fracture models at different scales through the wavelet transform. The size of the standard cylindrical rock sample is 100 mm × 50 mm. A binocular 3D laser scanning system was used to scan the Brazilian splitting test’s fracture surface. The binocular 3D laser scanning system’s accuracy is ±20μm, and the scanning interval is set at 0.1 mm on X-axis and 0.2 mm on Y-axis to obtain surface elevation data of the sample. We can get the digital fractured surface the same as the entire fracture surface through the points scanned. The vertical section of the middle of the specimen was taken as the sample’s typical fracture model’s original profile. The entire process of constructing an accurate rough fracture profile is shown in Figure 1.Figure 1
Generating process of the single fracture profile of coarse-grained granite: (a) specimen after Brazilian splitting test; (b) fracture surface formed by splitting; (c) 3D laser scanning; (d) fracture surface generated by scanning; (e) single rough fracture profile.
(a)(b)(c)(d)(e)


## 2.2. Basic Theory of the Wavelet Transform

The single rough fracture profile has characteristics of multiple scales and is similar to signals. Therefore, the 2D single rough fracture profile can be regarded as a signal functionft to be processed, where f represents the vertical height of fluctuations and t represents the horizontal distance along the fracture. When applying the wavelet transform, equation (1) is used.(1)Wfa,b=f,ψa,b=1a∫−∞+∞ftψt−ba¯dt,(2)ψta,b=1aψt−ba,a,b∈R;a≠0,where ψt is the mother wavelet function; ψa,bt is a wavelet sequence function; a is the scale factor; b is the translation factor. Therefore, the wavelet transform can be seen as a convolution of the signal function ft and a wavelet sequence function ψa,bt with the bandpass filter. A detailed discussion of the wavelet transform can be found in Zou et al.’s work [27].Before performing the wavelet transform, the appropriate mother wavelet function needs to be chosen first. Which wavelet to be chosen as the mother wavelet mainly depends on the application requirements [32]. Mirzaeian [33] pointed out that high-order Daubechies wavelets are more suitable for processing surface roughness than other wavelets in wavelet libraries. Sharifzadeh [25], Zou et al. [26], Zou et al. [27], and Wang et al. [28] used Db8 wavelets in the Daubechies wavelet series to process 2D rough single fracture profiles or 3D single rough fracture surfaces. Therefore, the Db8 wavelet was used as the mother wavelet function in the wavelet transform.


## 2.3. Multiscale Decomposition of Single Rough Fracture Profiles

In this section, the wavelet transform and reconstruction will be performed on the single rough fracture in Figure1(e) as the original fracture model S. The profile of the single rough fracture is defined as fx,x∈0,L, where L represents the length of the profile along the X-axis. In this paper, L = 100 mm. F is the height of the measured profile at the position of x. For the fracture profile, the Mallat algorithm is used for eight-level decomposition and reconstruction [34]. The approximate models are labeled as A1–A8, and the detailed models are labeled as D1–D8. The specific decomposition process labeled can be found in Zou et al.’s work [27].Figure2 shows the approximate models Ai and the detailed models Di of each decomposition level. The approximate models Ai represents the geometry that plays a dominant role in the fracture profile with a low frequency, large wavelength, and large amplitude. The detailed models Di represents small-scale ups and downs with a high frequency, small wavelength, and small amplitude compared with the approximate profile. As the level of decomposition increases, it can be seen that the wavelength and amplitude of Di gradually increase, varying from 0.05 mm in the first level to more than 0.5 mm in the eighth level. So, it gradually shows a large-scale characteristic (low frequency and large wavelength), which is similar to the entire profile. At the same time, as the level of decomposition increases, the approximate profile (blue line) gradually loses the original profile information compared with the entire profile, leading to the larger wavelength and smaller amplitude, and it slowly turns smooth and parallel.Figure 2
The approximate profiles Ai and detailed profiles Di under different decomposition levels: (a) Level 1. (b) Level 2. (c) Level 3. (d) Level 4. (e) Level 5. (f) Level 6. (g) Level 7. (h) Level 8.
(a)(b)(c)(d)(e)(f)(g)(h)According to Figure2, the approximate profile A4 shape remains the same as the original profile S. The shape of the approximate profile A5 is different from the original profile S. The fluctuation difference Δ (the highest vertical height minus the lowest vertical height) under different decomposition levels is counted, as shown in Figure 3. The inflection point of the curve of fluctuation difference Δ and decomposition levels is at level four. The value of fluctuation difference Δ before and after the fourth level changes significantly. Therefore, in this paper, the approximate profile A4 is taken as the primary roughness (wavy fluctuation structure). The original profile S minus the approximate profile A4 (D1 + D2 + D3 + D4) is considered as the secondary roughness (fine roughness structure).Figure 3
The variation of fluctuations difference Δ with decomposition level.The wavelet transform is used to decompose and reconstruct the rough fracture profiles, separating the fractured profiles into the generally approximate profiles and the locally detailed profiles in a unique way. It establishes the one-to-one relationship between locally detailed profiles, approximate profiles, and the original profile, making the roughness of the fracture profile unique and accurate at different scales [7, 26–28]. It also provides a new solution for studying the geometric characteristics of the fracture roughness at different scales. At the same time, it also lays the geometric foundation for analyzing the effect of the detailed profile of rough fractures significantly when the secondary roughness affects the discrepancy between the forward and reverse hydraulic characteristics.


## 3. Numerical Simulation of the Forward and Reverse Seepage in the Single Rough Fracture

### 3.1. Governing Equations and Numerical Methods

TheN-S equation and the continuity equation (momentum conservation and mass conservation) are often used to describe fluid flow in fractures quantitatively [35–37]. The equation of momentum conservation and mass conservation of the incompressible Newtonian fluid flow in a stable state can be written as follows:(3)ρu⟶⋅∇u⟶=−∇P+μ∇2u⟶,∇⋅u⟶=0,where ρ is the density of the fluid, u⟶ is the velocity vector, μ is viscosity, P is the total pressure, ∇ is Hamiltonian operator, and ∇2 is Laplace operator.Numerical methods are widely used to solve theN-S equation to study the flow characteristics in rough fractures because it cannot be solved directly [15, 17, 35, 37–39]. In this paper, Fluent 16.0 developed based on FVM code was adopted to simulate the seepage in rough fractures. Javadi et al. [36], Qian et al. [40], Qian et al. [41], and Liu et al. [42] have verified the reliability of the software and studied the Darcy or non-Darcy fluid flow in the fracture. In this study, DNS (direct numerical simulation) and RANS (Reynolds average Navier–Stokes equation) are used for the simulation of laminar and turbulent fluid flow, respectively, and the standard k−ε model is used to simulate turbulent fluid flow. The detailed theory can be seen in Javadi et al.’s work [36].


### 3.2. Calculation Area and Boundary Conditions

In order to study the influence of the detailed profile on the forward and reverse seepage characteristics of single rough fractures, the original profile model S and the eight approximate profile models Ai are moved upward by 0.5 mm, respectively, building nine rough fracture models at different scales (different decomposition levels). Furthermore, the model is discretized by a quadrilateral mesh. The average side length of the element is 0.02 mm, which is the same as the fracture surface’s scanning accuracy.For comparative studies, the same hydraulic boundary conditions are assigned to each fracture model. The entrance boundary is set as the velocity inlet boundary. The flow rate values corresponding to 13 sets of Reynolds number values (0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100, 500, and 1000) are taken. These Reynolds numbers are chosen because they are common in laboratory conditions and water conservancy project [27]. The Reynolds number is defined as the ratio of inertial force to viscous force. For fluid flow in rough fractures, it is given by Zimmerman et al. [43], Javadi et al. [44], and Zou et al. [27] as follows:(4)Re=ρQμw=ρveinμ,where Q is the flow rate, w is the width perpendicular to the direction of the fracture model surface (in 2D conditions,w = 1m), v is the mean value of the velocity along the +X-direction, ein is the aperture of the inlet boundary, the density of water ρ = 1000 kg/m3, the coefficient of dynamic viscosity μ = 1 × 10−3Pa·s, and the effect of the gravity is ignored.In order to compare and analyze the influence of flow direction on the seepage field, the left side is set as the inlet and the right side is set as the outlet when simulating forward seepage. On the contrary, when simulating reverse flow, the right side is set as the inlet and the left side as the outlet. In this paper, the fracture profile is vertically translated upward to create a fracture model with seepage channels, soein of all fracture models equals 0.5 mm. The velocity of entrance boundary is set at a series of values between 2.00962×10−6m/s and 2.00962 m/s, making the Reynolds number equal 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100, 500, and 1000. At the same time, the exit boundary of all fracture models is set at zero pressure. The fracture models’ upper and lower boundaries are set as boundary conditions without fluid flow and displacement.


## 3.1. Governing Equations and Numerical Methods

TheN-S equation and the continuity equation (momentum conservation and mass conservation) are often used to describe fluid flow in fractures quantitatively [35–37]. The equation of momentum conservation and mass conservation of the incompressible Newtonian fluid flow in a stable state can be written as follows:(3)ρu⟶⋅∇u⟶=−∇P+μ∇2u⟶,∇⋅u⟶=0,where ρ is the density of the fluid, u⟶ is the velocity vector, μ is viscosity, P is the total pressure, ∇ is Hamiltonian operator, and ∇2 is Laplace operator.Numerical methods are widely used to solve theN-S equation to study the flow characteristics in rough fractures because it cannot be solved directly [15, 17, 35, 37–39]. In this paper, Fluent 16.0 developed based on FVM code was adopted to simulate the seepage in rough fractures. Javadi et al. [36], Qian et al. [40], Qian et al. [41], and Liu et al. [42] have verified the reliability of the software and studied the Darcy or non-Darcy fluid flow in the fracture. In this study, DNS (direct numerical simulation) and RANS (Reynolds average Navier–Stokes equation) are used for the simulation of laminar and turbulent fluid flow, respectively, and the standard k−ε model is used to simulate turbulent fluid flow. The detailed theory can be seen in Javadi et al.’s work [36].


## 3.2. Calculation Area and Boundary Conditions

In order to study the influence of the detailed profile on the forward and reverse seepage characteristics of single rough fractures, the original profile model S and the eight approximate profile models Ai are moved upward by 0.5 mm, respectively, building nine rough fracture models at different scales (different decomposition levels). Furthermore, the model is discretized by a quadrilateral mesh. The average side length of the element is 0.02 mm, which is the same as the fracture surface’s scanning accuracy.For comparative studies, the same hydraulic boundary conditions are assigned to each fracture model. The entrance boundary is set as the velocity inlet boundary. The flow rate values corresponding to 13 sets of Reynolds number values (0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100, 500, and 1000) are taken. These Reynolds numbers are chosen because they are common in laboratory conditions and water conservancy project [27]. The Reynolds number is defined as the ratio of inertial force to viscous force. For fluid flow in rough fractures, it is given by Zimmerman et al. [43], Javadi et al. [44], and Zou et al. [27] as follows:(4)Re=ρQμw=ρveinμ,where Q is the flow rate, w is the width perpendicular to the direction of the fracture model surface (in 2D conditions,w = 1m), v is the mean value of the velocity along the +X-direction, ein is the aperture of the inlet boundary, the density of water ρ = 1000 kg/m3, the coefficient of dynamic viscosity μ = 1 × 10−3Pa·s, and the effect of the gravity is ignored.In order to compare and analyze the influence of flow direction on the seepage field, the left side is set as the inlet and the right side is set as the outlet when simulating forward seepage. On the contrary, when simulating reverse flow, the right side is set as the inlet and the left side as the outlet. In this paper, the fracture profile is vertically translated upward to create a fracture model with seepage channels, soein of all fracture models equals 0.5 mm. The velocity of entrance boundary is set at a series of values between 2.00962×10−6m/s and 2.00962 m/s, making the Reynolds number equal 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100, 500, and 1000. At the same time, the exit boundary of all fracture models is set at zero pressure. The fracture models’ upper and lower boundaries are set as boundary conditions without fluid flow and displacement.


## 4. Evolution of the Discrepancy between Forward and Reverse Seepages

### 4.1. Evolution of Equivalent Hydraulic Aperture

The simplest and most well-known model for describing fluid flow in a single rock fracture is the smooth and parallel plate model, which simplifies the rock fracture into two smooth parallel plates with a gap between the upper and lower surfaces. When the flow rate is low, the state of the fluid flow in the smooth and parallel plate model is laminar. In the incompressible Newtonian fluid flow, the quantitative relationship (equation (5)) between the flow rate Q and the cubic of the aperture e can be obtained under the constant pressure gradient, which is the famous cubic law:(5)Q=−we312μ∇P.Because the natural fracture surface is far from the smooth plane, equation (5) needs to be revised before the cubic law is applied to rough fractures, where the aperture will be replaced by the equivalent hydraulic aperture eh. The equivalent hydraulic aperture eh is acquired based on the Darcy flow state’s back-calculation, reflecting the rough fracture’s discharge capacity in the Darcy flow state. The cubic law is used to calculate the equivalent hydraulic aperture eh under all selected Reynolds number, and it can also be used to quantitatively describe the discharge capacity of the rough fracture under different Reynolds numbers [26, 28, 45]. In order to reflect the effect of fracture roughness and nonlinear seepage on the water conductivity of fracture, the equivalent hydraulic aperture is normalized (the ratio of the equivalent hydraulic aperture eh to the average mechanical aperture em). Li and Jiang [46], Qian et al. [47], and Zhang et al. [48] have established the relationship between em/eh, roughness, and Reynolds number Re. Besides, under the same roughness, em/eh will increase with the increase of Re. Cao et al. [11] also proved by experiments that increasing hydraulic pressure decreased hydraulic conductivity and em/eh discrepancy. Briggs et al. [45] define eh/em as the relative effective aperture, and the average mechanical aperture em in this paper is 0.5 mm.When the Reynolds number varies from 0.001 to 1000, the trend of the relative effective apertureeh/em of the original profile model and the eight approximate profile models is shown in Figure 4. When the Reynolds number is between 0.001 and 10, the values of eh/em of all models are the same in all models (meaning the equivalent hydraulic aperture eh is a constant). So, the fluid state at this time should be within the laminar flow. When the Reynolds number is over 10, the value of eh/em of all models begins to decrease, indicating that the relationship between flow rate Q and pressure gradient ∇P deviates from the cubic law. The fluid flow has entered a nonlaminar state (nonlinear flow). As the Reynolds number increases, the value of eh/em decreases gradually, suggesting that as the Reynolds number increases, the characteristics of nonlinearity in the fracture seepage increase. This is consistent with the research results of Zhang et al. [48]. Qin et al. [49] found that the modified cubic law can be well adopted in the test when the hydraulic pressure is 0.1 MPa, but R2 of the fitting lines decreases with the increase of hydraulic pressure. The influence of nonlinear flow in the crack increases with the increase of hydraulic pressure, which narrows the effective water passage and then reduces the actual flow capacity of the crack. Meanwhile, as the number of decomposition levels increases, the detailed profiles at different scales are gradually removed, and the value of eh/em in the Darcy flow state gradually increases. The value of eh/em at the 7th and 8th decomposition levels is approaching 1, which is close to the seepage results of the smooth and parallel plate model.Figure 4
Variation of the relative hydraulic apertureeh/em with Reynolds number under different decomposition levels in the approximate profile models (Ai+ and Ai− represent the calculation results of model Ai in the forward and reverse directions, respectively).Although the difference between the forward and reverse seepages cannot be clearly demonstrated in Figure4, the relationship between relative effective aperture eh/emC and the Reynolds number in two directions under the same Reynolds number can be obtained through the difference between the values of eh/em in two directions. As is shown in Figure 5, when the Reynolds number varies from 0.001 to 10, the value of eh/emC is basically zero, suggesting that the seepage property of the fracture is only related to the geometry of the fracture in the laminar (linear) flow, no matter how the roughness of the fracture changes. Besides, the flow direction of the fluid in the fracture will not make a difference in the seepage characteristics. When Re > 10, the fluid state of the fluid in the fracture enters the nonlinear stage, and the values of eh/emC of the approximate models of decomposition levels ∼4 are not zero and increase with the increase of the Reynolds number. Meanwhile, the values of eh/emC of the approximate models of decomposition levels 5∼8 are almost zero and do not change with the Reynolds number. In conclusion, the existence of nonlinear seepage does not necessarily result in the discrepancy between the two directions’ hydraulic characteristics. However, the discrepancy between the hydraulic characteristics of the two directions must be due to nonlinear seepage.Figure 5
Variation of the difference of relative hydraulic apertureeh/emC with Reynolds number under different decomposition levels in the approximate profile models.


### 4.2. Evolution of Non-Darcy Coefficient

It can be seen from Section4.1 that as the Reynolds number increases, the relationship between the flow rate Q and the pressure gradient ∇P is no longer linear, and equation (5) is not applicable. In describing the nonlinear relationship between pressure gradient and flow rate, the Forchheimer equation [50] is the most widely accepted:(6)−∇P=AQ+BQ2,where A and B are the coefficients of the linear and nonlinear terms. Cao et al. [11] pointed out that the nonlinearity of fluid flow caused by inertia force had a remarkable influence on hydraulic behavior with high hydraulic pressure. When the Reynolds number is small, the inertia effect can be neglected (the nonlinear term BQ2 can be neglected). In this case, equation (6) evolves into equation (5). From the changing trend of eh/em with the Reynolds number, it is appropriate to adopt equation (6) to describe the relationship between pressure gradient ∇P and flow rate Q. So we add Reynolds numbers of 300 and 750 among 10, 50, 100, 500, and 1000 to increase the accuracy of the fitted coefficients A and B. It can be drawn from Figure 6 that the Forchheimer equation can well reflect the nonlinear relationship between ∇P and Q. The seepage coefficients in the forward direction (A+, B+) and reverse direction (A−, B−) under different decomposition levels and their regression coefficients are listed in Table 1. As is shown in Figure 6, with the increase of decomposition level, the nonlinear characteristic between ∇P and Q gradually weakens and finally tends to be stable. The pressure gradient ∇P is obviously different in two directions when the original profile model to the third approximate profile model shares the same flow rate boundary, while ∇P in two directions in the 4th∼8th approximate profile models are basically the same.Figure 6
Analysis of nonlinear seepage results.
(a)(b)Table 1
Statistical table of non-Darcy fitted parameters at different scales.
A+B+R2β+A−B−R2β−(×109 kg·s−1·m−5)(×1012 kg·m−8)(m−1)(×109 kg·s−1·m−5)(×1012 kg·m−8)(m−1)S0.5391.8250.9999201.8890.5591.6680.9999184.502A10.5501.3740.9998155.9080.5761.2470.9998141.403A20.4660.9170.9997119.0030.4760.8130.9997105.455A30.3690.6110.999791.1690.3840.5340.999679.607A40.2430.2530.999549.1330.2420.2560.999649.628A50.1800.1590.999636.2600.1800.1590.999636.244A60.1610.1360.999633.3870.1610.1360.999633.313A70.1580.1340.999633.2660.1580.1340.999633.214A80.1570.1330.999633.1100.1570.1320.999633.075For single fractures, the coefficients can be obtained by the two following equations [51, 52]:(7a)A=μkAh=12μweh3,(7b)B=βρAh2=βρw2eh2,where k is the hydraulic conductivity of rock fracture, β is the non-Darcy coefficient, and Ah is the area of equivalent cross section. The linear term coefficient A represents the fracture’s inherent permeability, which is only related to the equivalent hydraulic aperture eh. It is clearly shown in Figure 5 that the equivalent hydraulic aperture of the forward and reverse directions is almost the same in the Darcy flow state. Therefore, the coefficients of A in two directions should also be the same. Although there are differences between A+ and A−, which are fitted by equation (6), the primary source of nonlinear seepage comes from the nonlinear term. The difference of seepage in two directions is also derived from the nonlinear term. Therefore, the difference of coefficient A is not discussed and studied here; coefficient B is mainly focused instead, as opposed to porous media where both coefficients A and B should be simultaneously parameterized [51, 53]. The difference between forward and reverse seepages can only result from the non-Darcy coefficient β, so the non-Darcy coefficient β can be used to quantitatively characterize the discrepancy of seepage in two directions in the rough fracture. According to the fitted B+ and B−, the non-Darcy coefficients β+ and β−, which are the non-Darcy coefficients in two directions, can be calculated, respectively, by equation (7b). The results can be found in Table 1.The variation of the non-Darcy coefficientsβ+ and β− with the decomposition level is plotted in Figure 7(a). As the decomposition level increases, β+ and β− decrease gradually first and then tend to be stable when the decomposition level is over 4. The non-Darcy coefficient β+ is significantly larger than β− when the decomposition level is from 0 to 3. βC is defined as the difference between β+ and β−, and the variation of βC with decomposition level is plotted in Figure 7(b). With the increase of the decomposition level, βC is slightly decreased. βC of decomposition levels 4 and 5∼8 are −0.43 and almost zero, respectively, while that of decomposition levels 1∼3 is over 10, suggesting that the seepage direction no longer affects the fracture’s hydraulic characteristics in the 4th approximate profile (primary roughness) models, and the secondary roughness may be the geometric basis for the difference between seepages in two directions.Figure 7
The correlation between Non-Darcy coefficientsβ in the forward and reverse directions with the decomposition level: (a) non-Darcy coefficient β; (b) difference of the non-Darcy coefficient βC.
(a)(b)


## 4.1. Evolution of Equivalent Hydraulic Aperture

The simplest and most well-known model for describing fluid flow in a single rock fracture is the smooth and parallel plate model, which simplifies the rock fracture into two smooth parallel plates with a gap between the upper and lower surfaces. When the flow rate is low, the state of the fluid flow in the smooth and parallel plate model is laminar. In the incompressible Newtonian fluid flow, the quantitative relationship (equation (5)) between the flow rate Q and the cubic of the aperture e can be obtained under the constant pressure gradient, which is the famous cubic law:(5)Q=−we312μ∇P.Because the natural fracture surface is far from the smooth plane, equation (5) needs to be revised before the cubic law is applied to rough fractures, where the aperture will be replaced by the equivalent hydraulic aperture eh. The equivalent hydraulic aperture eh is acquired based on the Darcy flow state’s back-calculation, reflecting the rough fracture’s discharge capacity in the Darcy flow state. The cubic law is used to calculate the equivalent hydraulic aperture eh under all selected Reynolds number, and it can also be used to quantitatively describe the discharge capacity of the rough fracture under different Reynolds numbers [26, 28, 45]. In order to reflect the effect of fracture roughness and nonlinear seepage on the water conductivity of fracture, the equivalent hydraulic aperture is normalized (the ratio of the equivalent hydraulic aperture eh to the average mechanical aperture em). Li and Jiang [46], Qian et al. [47], and Zhang et al. [48] have established the relationship between em/eh, roughness, and Reynolds number Re. Besides, under the same roughness, em/eh will increase with the increase of Re. Cao et al. [11] also proved by experiments that increasing hydraulic pressure decreased hydraulic conductivity and em/eh discrepancy. Briggs et al. [45] define eh/em as the relative effective aperture, and the average mechanical aperture em in this paper is 0.5 mm.When the Reynolds number varies from 0.001 to 1000, the trend of the relative effective apertureeh/em of the original profile model and the eight approximate profile models is shown in Figure 4. When the Reynolds number is between 0.001 and 10, the values of eh/em of all models are the same in all models (meaning the equivalent hydraulic aperture eh is a constant). So, the fluid state at this time should be within the laminar flow. When the Reynolds number is over 10, the value of eh/em of all models begins to decrease, indicating that the relationship between flow rate Q and pressure gradient ∇P deviates from the cubic law. The fluid flow has entered a nonlaminar state (nonlinear flow). As the Reynolds number increases, the value of eh/em decreases gradually, suggesting that as the Reynolds number increases, the characteristics of nonlinearity in the fracture seepage increase. This is consistent with the research results of Zhang et al. [48]. Qin et al. [49] found that the modified cubic law can be well adopted in the test when the hydraulic pressure is 0.1 MPa, but R2 of the fitting lines decreases with the increase of hydraulic pressure. The influence of nonlinear flow in the crack increases with the increase of hydraulic pressure, which narrows the effective water passage and then reduces the actual flow capacity of the crack. Meanwhile, as the number of decomposition levels increases, the detailed profiles at different scales are gradually removed, and the value of eh/em in the Darcy flow state gradually increases. The value of eh/em at the 7th and 8th decomposition levels is approaching 1, which is close to the seepage results of the smooth and parallel plate model.Figure 4
Variation of the relative hydraulic apertureeh/em with Reynolds number under different decomposition levels in the approximate profile models (Ai+ and Ai− represent the calculation results of model Ai in the forward and reverse directions, respectively).Although the difference between the forward and reverse seepages cannot be clearly demonstrated in Figure4, the relationship between relative effective aperture eh/emC and the Reynolds number in two directions under the same Reynolds number can be obtained through the difference between the values of eh/em in two directions. As is shown in Figure 5, when the Reynolds number varies from 0.001 to 10, the value of eh/emC is basically zero, suggesting that the seepage property of the fracture is only related to the geometry of the fracture in the laminar (linear) flow, no matter how the roughness of the fracture changes. Besides, the flow direction of the fluid in the fracture will not make a difference in the seepage characteristics. When Re > 10, the fluid state of the fluid in the fracture enters the nonlinear stage, and the values of eh/emC of the approximate models of decomposition levels ∼4 are not zero and increase with the increase of the Reynolds number. Meanwhile, the values of eh/emC of the approximate models of decomposition levels 5∼8 are almost zero and do not change with the Reynolds number. In conclusion, the existence of nonlinear seepage does not necessarily result in the discrepancy between the two directions’ hydraulic characteristics. However, the discrepancy between the hydraulic characteristics of the two directions must be due to nonlinear seepage.Figure 5
Variation of the difference of relative hydraulic apertureeh/emC with Reynolds number under different decomposition levels in the approximate profile models.


## 4.2. Evolution of Non-Darcy Coefficient

It can be seen from Section4.1 that as the Reynolds number increases, the relationship between the flow rate Q and the pressure gradient ∇P is no longer linear, and equation (5) is not applicable. In describing the nonlinear relationship between pressure gradient and flow rate, the Forchheimer equation [50] is the most widely accepted:(6)−∇P=AQ+BQ2,where A and B are the coefficients of the linear and nonlinear terms. Cao et al. [11] pointed out that the nonlinearity of fluid flow caused by inertia force had a remarkable influence on hydraulic behavior with high hydraulic pressure. When the Reynolds number is small, the inertia effect can be neglected (the nonlinear term BQ2 can be neglected). In this case, equation (6) evolves into equation (5). From the changing trend of eh/em with the Reynolds number, it is appropriate to adopt equation (6) to describe the relationship between pressure gradient ∇P and flow rate Q. So we add Reynolds numbers of 300 and 750 among 10, 50, 100, 500, and 1000 to increase the accuracy of the fitted coefficients A and B. It can be drawn from Figure 6 that the Forchheimer equation can well reflect the nonlinear relationship between ∇P and Q. The seepage coefficients in the forward direction (A+, B+) and reverse direction (A−, B−) under different decomposition levels and their regression coefficients are listed in Table 1. As is shown in Figure 6, with the increase of decomposition level, the nonlinear characteristic between ∇P and Q gradually weakens and finally tends to be stable. The pressure gradient ∇P is obviously different in two directions when the original profile model to the third approximate profile model shares the same flow rate boundary, while ∇P in two directions in the 4th∼8th approximate profile models are basically the same.Figure 6
Analysis of nonlinear seepage results.
(a)(b)Table 1
Statistical table of non-Darcy fitted parameters at different scales.
A+B+R2β+A−B−R2β−(×109 kg·s−1·m−5)(×1012 kg·m−8)(m−1)(×109 kg·s−1·m−5)(×1012 kg·m−8)(m−1)S0.5391.8250.9999201.8890.5591.6680.9999184.502A10.5501.3740.9998155.9080.5761.2470.9998141.403A20.4660.9170.9997119.0030.4760.8130.9997105.455A30.3690.6110.999791.1690.3840.5340.999679.607A40.2430.2530.999549.1330.2420.2560.999649.628A50.1800.1590.999636.2600.1800.1590.999636.244A60.1610.1360.999633.3870.1610.1360.999633.313A70.1580.1340.999633.2660.1580.1340.999633.214A80.1570.1330.999633.1100.1570.1320.999633.075For single fractures, the coefficients can be obtained by the two following equations [51, 52]:(7a)A=μkAh=12μweh3,(7b)B=βρAh2=βρw2eh2,where k is the hydraulic conductivity of rock fracture, β is the non-Darcy coefficient, and Ah is the area of equivalent cross section. The linear term coefficient A represents the fracture’s inherent permeability, which is only related to the equivalent hydraulic aperture eh. It is clearly shown in Figure 5 that the equivalent hydraulic aperture of the forward and reverse directions is almost the same in the Darcy flow state. Therefore, the coefficients of A in two directions should also be the same. Although there are differences between A+ and A−, which are fitted by equation (6), the primary source of nonlinear seepage comes from the nonlinear term. The difference of seepage in two directions is also derived from the nonlinear term. Therefore, the difference of coefficient A is not discussed and studied here; coefficient B is mainly focused instead, as opposed to porous media where both coefficients A and B should be simultaneously parameterized [51, 53]. The difference between forward and reverse seepages can only result from the non-Darcy coefficient β, so the non-Darcy coefficient β can be used to quantitatively characterize the discrepancy of seepage in two directions in the rough fracture. According to the fitted B+ and B−, the non-Darcy coefficients β+ and β−, which are the non-Darcy coefficients in two directions, can be calculated, respectively, by equation (7b). The results can be found in Table 1.The variation of the non-Darcy coefficientsβ+ and β− with the decomposition level is plotted in Figure 7(a). As the decomposition level increases, β+ and β− decrease gradually first and then tend to be stable when the decomposition level is over 4. The non-Darcy coefficient β+ is significantly larger than β− when the decomposition level is from 0 to 3. βC is defined as the difference between β+ and β−, and the variation of βC with decomposition level is plotted in Figure 7(b). With the increase of the decomposition level, βC is slightly decreased. βC of decomposition levels 4 and 5∼8 are −0.43 and almost zero, respectively, while that of decomposition levels 1∼3 is over 10, suggesting that the seepage direction no longer affects the fracture’s hydraulic characteristics in the 4th approximate profile (primary roughness) models, and the secondary roughness may be the geometric basis for the difference between seepages in two directions.Figure 7
The correlation between Non-Darcy coefficientsβ in the forward and reverse directions with the decomposition level: (a) non-Darcy coefficient β; (b) difference of the non-Darcy coefficient βC.
(a)(b)


## 5. Mechanism of the Discrepancy between the Forward and Reverse Seepages

### 5.1. Discrepancy in Eddy Currents between Seepages in Two Directions

In this paper, the mechanical aperture is only 0.5 mm, while the model is 100 mm long, which results in the large aspect ratio of the overall model. Therefore, the eddy current distribution characteristics of the entire model cannot be displayed well. A local representative area is taken out in this section to study the influence of flow direction on the seepage field, revealing the discrepancy between the forward and reverse seepages. TheX-direction velocity contour and flow line distribution in the range of X = [51 mm, 53 mm] of the original model are selected, and the results of forward and reverse seepage are shown in Figure 8 when the Reynolds number is at 1, 10, 100, and 1000, respectively. When Re > 10 and only the forward (or reverse) seepage direction is analyzed, as the Reynolds number increases, the eddy currents are about to emerge and are gradually enlarged at the location where the profile fluctuates, which is consistent with the observations of Lee et al. [8] and Lee et al. [37] using the micro-PIV (micro particle image velocity) system in the laboratory. Zou et al. [27] and Briggs et al. [45] also obtained similar conclusions through numerical simulation. When Re = 1 and 10, the seepage in two directions was analyzed. There is no eddy current, and the distribution of velocity is the same, except that the velocity direction of the forward seepage is opposite to that of the reverse seepage. When Re = 100, there are eddy currents with different shapes and sizes in the forward and reverse seepage. The eddy current can be regarded as a “circulation zone” or “dead water zone,” which will reduce the flow channel inevitably. As the eddy current’s locations and sizes vary in two directions, the flow line and the velocity contour in two directions differ. When Re = 1000, the eddy zone becomes larger compared with that under Re = 100. Furthermore, where there is no eddy zone at Re = 100, a new eddy zone is created at Re = 1000. The expansion of the original eddy zone and the creation of a new eddy zone further compress the effective seepage channel. This has also led to a further increase in nonlinear seepage. Due to the inconsistent distribution of eddy zones, the difference between forward and reverse seepages has also increased. Cao et al. [11] pointed out that the discrepancy and em/eh ratio were influenced by recirculation zones, infilling materials, roughness, contact areas, and inertia force. A two-dimensional numerical model showed that the vortices were the important cause of damage to the radiation flow [54].Figure 8
X-direction velocity contour and flow line of the original profile model under different Reynolds numbers in the forward and reverse directions (the contour represents the velocity distribution, m/s). (a) Re = 1. (b) Re = 10. (c) Re = 100. (d) Re = 1000.
(a)(b)(c)(d)In conclusion, the “dead water zone” generated by the eddy current significantly reduces the flow channel, which makes the effective flow channel narrow and then makes the distribution of the velocity and flow line of theX-direction more complex. So, the eddy current’s existence has a decisive effect on the existence and development of nonlinear seepage [27, 28, 45]. The difference between the eddy current distribution shapes and locations in two directions is the main reason for the discrepancy between the nonlinear seepages in two directions. The eddy current’s location and region are related to the direction of the seepage and the Reynolds number. With the increase of the Reynolds number, the eddy current region becomes larger, and even new eddy currents are generated. It is indicated that the difference between non-Darcy seepages in two directions will be enlarged with the increase of the number.


### 5.2. Effect of the Detailed Profile on Eddy Current Generated by Seepage in Two Directions

To further study the non-Darcy seepage characteristics of two directions in the approximate profile model at different decomposition levels, the velocity and flow line distribution in the local region ofX = [51 mm, 53 mm], as shown in Figure 9, is selected. When only the forward (reverse) seepage is analyzed, with the increase of the decomposition level, the region of the eddy current is gradually narrowed [26] compared with the original profile model (Figure 8(d)). At the fourth decomposition level, there is not any eddy. Therefore, the secondary roughness plays a leading role in the generation of the eddy current, which is following the results of Zou et al. [27], Wang et al. [28], and Dou et al. [7]. When the forward and reverse seepage fields are compared and analyzed, with the increase of the decomposition level, the eddy current regions are both narrowed, and the flow line is parallel to the fracture profile. At the 8th decomposition level, the velocity contour and flow line are similar to those of the smooth and parallel plate models. Therefore, with the increase of the decomposition level, the eddy current shape and region decrease and disappear at last in two directions, resulting in the decrease of the discrepancy between non-Darcy seepages in two directions. Moreover, this discrepancy will vanish when the eddy current disappears.Figure 9
X-direction velocity and flow line of two directions in the approximate models when Re = 1000 (the contour represents the velocity distribution, m/s). (a) Level 1. (b) Level 2. (c) Level 3. (d) Level 4. (e) Level 5. (f) Level 6. (g) Level 7. (h) Level 8.
(a)(b)(c)(d)(e)(f)(g)(h)In summary, as the decomposition level increases, the detailed information of the fracture profile containing is less, resulting in the shape and region of the eddy current decreasing or even disappearing. It is revealed that the detailed profile is decisive in the generation and development of the eddy current. From the results of the 1∼4th decomposition levels, it can be seen that the local concave and convex area generated by the detailed profile in the fracture provides a potential place for the generation of eddy current in the nonlinear seepage. Simultaneously, with the decrease of the range of the eddy current, the difference between non-Darcy seepages in two directions decreases, which further indicates that the discrepancy between the eddy currents in two directions is the primary cause of the discrepancy between non-Darcy seepages in two directions. Besides, the secondary roughness (D1 + D2 + D3 + D4) is related to the existence and development of eddy currents. Therefore, it can be concluded that the secondary roughness is the geometric basis for the discrepancy between non-Darcy seepages in two directions.


## 5.1. Discrepancy in Eddy Currents between Seepages in Two Directions

In this paper, the mechanical aperture is only 0.5 mm, while the model is 100 mm long, which results in the large aspect ratio of the overall model. Therefore, the eddy current distribution characteristics of the entire model cannot be displayed well. A local representative area is taken out in this section to study the influence of flow direction on the seepage field, revealing the discrepancy between the forward and reverse seepages. TheX-direction velocity contour and flow line distribution in the range of X = [51 mm, 53 mm] of the original model are selected, and the results of forward and reverse seepage are shown in Figure 8 when the Reynolds number is at 1, 10, 100, and 1000, respectively. When Re > 10 and only the forward (or reverse) seepage direction is analyzed, as the Reynolds number increases, the eddy currents are about to emerge and are gradually enlarged at the location where the profile fluctuates, which is consistent with the observations of Lee et al. [8] and Lee et al. [37] using the micro-PIV (micro particle image velocity) system in the laboratory. Zou et al. [27] and Briggs et al. [45] also obtained similar conclusions through numerical simulation. When Re = 1 and 10, the seepage in two directions was analyzed. There is no eddy current, and the distribution of velocity is the same, except that the velocity direction of the forward seepage is opposite to that of the reverse seepage. When Re = 100, there are eddy currents with different shapes and sizes in the forward and reverse seepage. The eddy current can be regarded as a “circulation zone” or “dead water zone,” which will reduce the flow channel inevitably. As the eddy current’s locations and sizes vary in two directions, the flow line and the velocity contour in two directions differ. When Re = 1000, the eddy zone becomes larger compared with that under Re = 100. Furthermore, where there is no eddy zone at Re = 100, a new eddy zone is created at Re = 1000. The expansion of the original eddy zone and the creation of a new eddy zone further compress the effective seepage channel. This has also led to a further increase in nonlinear seepage. Due to the inconsistent distribution of eddy zones, the difference between forward and reverse seepages has also increased. Cao et al. [11] pointed out that the discrepancy and em/eh ratio were influenced by recirculation zones, infilling materials, roughness, contact areas, and inertia force. A two-dimensional numerical model showed that the vortices were the important cause of damage to the radiation flow [54].Figure 8
X-direction velocity contour and flow line of the original profile model under different Reynolds numbers in the forward and reverse directions (the contour represents the velocity distribution, m/s). (a) Re = 1. (b) Re = 10. (c) Re = 100. (d) Re = 1000.
(a)(b)(c)(d)In conclusion, the “dead water zone” generated by the eddy current significantly reduces the flow channel, which makes the effective flow channel narrow and then makes the distribution of the velocity and flow line of theX-direction more complex. So, the eddy current’s existence has a decisive effect on the existence and development of nonlinear seepage [27, 28, 45]. The difference between the eddy current distribution shapes and locations in two directions is the main reason for the discrepancy between the nonlinear seepages in two directions. The eddy current’s location and region are related to the direction of the seepage and the Reynolds number. With the increase of the Reynolds number, the eddy current region becomes larger, and even new eddy currents are generated. It is indicated that the difference between non-Darcy seepages in two directions will be enlarged with the increase of the number.


## 5.2. Effect of the Detailed Profile on Eddy Current Generated by Seepage in Two Directions

To further study the non-Darcy seepage characteristics of two directions in the approximate profile model at different decomposition levels, the velocity and flow line distribution in the local region ofX = [51 mm, 53 mm], as shown in Figure 9, is selected. When only the forward (reverse) seepage is analyzed, with the increase of the decomposition level, the region of the eddy current is gradually narrowed [26] compared with the original profile model (Figure 8(d)). At the fourth decomposition level, there is not any eddy. Therefore, the secondary roughness plays a leading role in the generation of the eddy current, which is following the results of Zou et al. [27], Wang et al. [28], and Dou et al. [7]. When the forward and reverse seepage fields are compared and analyzed, with the increase of the decomposition level, the eddy current regions are both narrowed, and the flow line is parallel to the fracture profile. At the 8th decomposition level, the velocity contour and flow line are similar to those of the smooth and parallel plate models. Therefore, with the increase of the decomposition level, the eddy current shape and region decrease and disappear at last in two directions, resulting in the decrease of the discrepancy between non-Darcy seepages in two directions. Moreover, this discrepancy will vanish when the eddy current disappears.Figure 9
X-direction velocity and flow line of two directions in the approximate models when Re = 1000 (the contour represents the velocity distribution, m/s). (a) Level 1. (b) Level 2. (c) Level 3. (d) Level 4. (e) Level 5. (f) Level 6. (g) Level 7. (h) Level 8.
(a)(b)(c)(d)(e)(f)(g)(h)In summary, as the decomposition level increases, the detailed information of the fracture profile containing is less, resulting in the shape and region of the eddy current decreasing or even disappearing. It is revealed that the detailed profile is decisive in the generation and development of the eddy current. From the results of the 1∼4th decomposition levels, it can be seen that the local concave and convex area generated by the detailed profile in the fracture provides a potential place for the generation of eddy current in the nonlinear seepage. Simultaneously, with the decrease of the range of the eddy current, the difference between non-Darcy seepages in two directions decreases, which further indicates that the discrepancy between the eddy currents in two directions is the primary cause of the discrepancy between non-Darcy seepages in two directions. Besides, the secondary roughness (D1 + D2 + D3 + D4) is related to the existence and development of eddy currents. Therefore, it can be concluded that the secondary roughness is the geometric basis for the discrepancy between non-Darcy seepages in two directions.


## 6. Discussion

### 6.1. Directional Roughness and Non-Darcy Coefficient

We know that when a rock sample with a single rough fracture is subjected to a forward and reverse shear experiment, the forward and reverse shear strength is inconsistent, and this is related to the inconsistency of the forward and reverse roughness. Zhang et al. [55] pointed out that surface roughness plays a significant role in affecting fracture flow in the review paper. Hence, it is necessary to consider the effect of surface morphology or roughness on fluid flow, which can uncover underlying flow mechanisms inside fractures and advance our understanding of the fracture flow.It can be seen from the previous chapter that, in the nonlinear seepage stage, the forward and reverse seepage characteristics are inconsistent, and whether this is related to the forward and reverse roughness will be briefly discussed in this section. Belem et al. [56], Zhang et al. [57], Grasselli et al. [58], Grasselli [59], and Tatone and Grasselli [60] have proposed methods and formulas for describing directional roughness. In this paper, the directional roughness presented by Tatone and Grasselli [60] is selected as an index to measure the inconsistency of forward and reverse roughness.Grasselli et al. [58] proposed a method and parameter for calculating directional surface roughness based on the 3D scanning point cloud. Tatone and Grasselli [60] based on the method of Grasselli et al. [58] proposed the directional roughness measurement of 2D fracture profile. The method used to calculate the roughness of the 2D fracture profile direction is as follows: 1. Use the least-squares method to establish the best-fitting line and set it as the X-axis. 2. Connect adjacent points to form a series of short straight line segments. 3. Select the analysis direction and calculate each line segment’s inclination angle. For the 2D fracture profile, only two shear directions (forward and reverse) need to be considered. A length score Lθ∗ that is steeper than the gradually higher angle threshold θ∗ can be obtained based on each line segment’s inclination angle. The relationship between the two can be fitted by equation (8). The detailed theory and method can be found in the paper by Tatone and Grasselli [60]:(8)Lθ∗=L0θmax∗−θ∗θmax∗C,where L0 is the length ratio where the angle threshold is 0°, θmax∗ represents the maximum apparent inclination in the analysis direction, and C is a dimensionless fitting parameter, which characterizes the distribution’s shape. By integrating the area under the curve between Lθ∗ and the angle threshold θ∗, the roughness L0θmax∗/C+1 in the analysis direction can be obtained by the following equation:(9)L0∫0θmax∗θmax∗−θ∗θmax∗Cdθ∗=L0θmax∗C+1.Figure10 shows the forward and reverse θ∗ and Lθ∗ scatter points calculated by the Tatone and Grasselli [60] method of the original profile and 8 approximate profiles. The dimensionless parameter C can be determined using equation (8) for the curve fitting of the scatter point. The forward and reverse L0 and θmax∗ can also be obtained through calculation and statistics, and then the directional roughness of the forward and reverse L0θmax∗/C+1 can be determined and are listed in Table 2. Figures 11(a) and 11(b) are, respectively, plotted as curves of the calculated forward and reverse direction roughness L0θmax∗/C+1 and its difference value L0θmax∗/C+1C with the decomposition level. Figure 11 shows that the overall forward roughness L0θmax∗/C+1 is smaller than the reverse roughness. As the number of decomposition levels increases, the difference between the forward roughness and reverse roughness L0θmax∗/C+1C decreases. When it is decomposed to the fourth level, the roughness in the forward and reverse directions have tended to zero, and the subsequent fifth to eighth levels also tend to zero. This is contrary to the non-Darcy coefficient βC curve with the decomposition levels in Figure 7(b). The difference of forward and reverse non-Darcy coefficients βC and the difference of directional roughness L0θmax∗/C+1C under different decomposition levels corresponded. Figure 12 shows a good negative linear relationship between the difference of directional roughness L0θmax∗/C+1C and the difference of non-Darcy coefficient βC between forward and reverse directions. It can be seen from this section that the difference of directional roughness between forward and reverse directions can reflect the difference between forward and reverse non-Darcy flows (the difference of non-Darcy coefficients). However, it should also be seen that there are few samples used in this article, and the exact quantitative relationship between the difference of directional roughness and the difference of non-Darcy flow needs further research.Figure 10
Angle probability distribution and fitting curve of rough fracture profile under different decomposition levels. (a) Original. (b) Level 1. (c) Level 2. (d) Level 3. (e) Level 4. (f) Level 5. (g) Level 6. (h) Level 7. (i) Level 8.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Table 2
Statistical table of relevant parameters of direction roughnessL0θmax∗/C+1 under different decomposition levels.
ForwardReverseL0θmax∗/C+1 (o)L0θmax∗CR2L0θmax∗CR2ForwardReverseS0.45174.5061.9490.9980.54971.0281.9330.98811.39613.297A10.43574.6281.9530.9980.56558.5051.5490.99710.99512.970A20.44671.8461.9760.9980.55451.7931.2960.99710.76812.496A30.46867.5921.9910.9990.53256.0961.5840.99510.57511.547A40.53446.5591.6310.9980.46648.8091.4250.9979.4499.379A50.51543.3633.2330.9930.48534.3272.0860.9985.2765.395A60.52811.0211.0330.9930.47216.8431.7840.9942.8632.855A70.29010.5000.7970.9880.7106.7501.7760.9581.6951.726A80.3885.3180.9720.9120.6123.0610.7710.9681.0461.058Figure 11
The relationship between the roughness in the forward and reverse directions (a)L0θmax∗/C+1 (b) and its difference L0θmax∗/C+1C with the decomposition level.
(a)(b)Figure 12
Relationship between the difference of non-Darcy coefficientβC and the difference of directional roughness L0θmax∗/C+1C.


### 6.2. Single Rough Fracture in Fracture Network

The single fracture is the basic unit of the fracture network, and the seepage characteristics of the single fracture will also affect the seepage characteristics of the fracture network. The research on the difference between the forward and reverse seepages in single rough fracture can lay the foundation for future research on the forward and reverse seepage characteristics of the rough fracture network. When it is considered that the single fractures in the fracture network obey the cubic law, the seepage characteristics of the fracture network can be described by the anisotropic Darcy’s law, as shown in the following equation:(10)uj=−KjiJi,where uj is the velocity in direction j, Ji is the hydraulic gradient in direction i, and Kji is the equivalent permeability coefficient tensor. As shown in this paper’s research results, in the case of Darcy seepage in rough single fractures, there is no difference in the seepage characteristics of the forward and reverse seepage directions, and both strictly obey the cubic law. In the Darcy seepage state, the permeability coefficient K of the fracture network composed of rough single fractures in the relative seepage direction (such as 0° and 180°, 90° and 270°) is also consistent, which conforms to the characteristics of a tensor. Thus, the equivalent permeability coefficient tensor Kij can express the permeability characteristics of the seepage in the Darcy seepage state.With the increase of hydraulic gradient, not only the single fracture but also the fracture network will enter non-Darcy seepage. For the fracture network, Liu et al. [42], Liu et al. [61], and Yin et al. [62] point out that when the hydraulic gradient is greater than a certain value, the ratio of flow velocity to hydraulic gradient is no longer constant, which no longer satisfies Darcy’s law. The seepage flow in the fracture network will also enter the state of non-Darcy seepage. Similar to the use of Forchheimer law (equation (6)) to describe the non-Darcy flow in a single rough fracture, the non-Darcy flow in a fracture network is also described using the Forchheimer law (equation (11)) of anisotropic seepage [63]:(11)−Ji=Aijuj+Bijuuj=Kij∗Kuj+βijguuj.In the equation,Aij and Bij are, respectively, the coefficient tensors of the viscous force term and the inertial force term, which can also be expressed as a matrix related to the permeability coefficient tensor Kij and the non-Darcy coefficient tensor βij,u is the modulus of flow velocity, Kij∗/K is the inverse matrix of the equivalent permeability coefficient matrix, and βij is the equivalent non-Darcy coefficient tensor.Since the forward and reverse seepage characteristics of a single smooth fracture will not differ between Darcy and non-Darcy seepage conditions, for a fracture network composed of smooth single fractures, the fracture network’s seepage characteristics under forward and reverse seepage will not be different. In other words, the permeability coefficientK of the smooth fracture network will remain consistent in any state under the relative seepage direction (forward and reverse direction). For smooth fracture networks, the equivalent non-Darcy coefficient β in the relative seepage direction obtained by equation (11) will also be consistent. The equivalent non-Darcy coefficient β conforms to the tensor’s symmetry, which can theoretically constitute the equivalent non-Darcy coefficient tensor βij [64]. However, this paper’s research shows that when Re > 10, the fracture seepage of a single rough fracture enters the non-Darcy seepage stage, and then the difference of forward and reverse seepage occurs. The quantitative description by th Forchheimer law (equation (9)) shows that, in the state of non-Darcy seepage, the permeability coefficient k has nothing to do with the direction of seepage, and the non-Darcy coefficient β is related to the direction of seepage. Therefore, for the fracture network composed of single rough fractures, the difference of the forward and reverse flow direction of the single rough fracture will also lead to the difference in the rough fracture network’s relative flow direction. If equation (11) is used, the equivalent non-Darcy coefficient β in the relative seepage direction will no longer be consistent and will no longer conform to the tensor’s symmetry. So, for rough fracture networks, the applicability of equation (11) will be limited. This will also raise new problems and challenges for describing the rough fracture network’s non-Darcy seepage characteristics accurately and quantitatively.In the actual natural environment and engineering environment, the single fractures that make up the fracture network in the rock mass are rough. In the non-Darcy seepage state of these rough single fractures, the forward and reverse seepage directions’ seepage characteristics will be different. How much do these differences affect the difference in the direction of seepage flow in the fracture network? In turn, will the difference in the flow direction of the rough fracture network for seepage flow and solute movement have much impact on the project? These issues need to be addressed by an in-depth research in the follow-up work.


## 6.1. Directional Roughness and Non-Darcy Coefficient

We know that when a rock sample with a single rough fracture is subjected to a forward and reverse shear experiment, the forward and reverse shear strength is inconsistent, and this is related to the inconsistency of the forward and reverse roughness. Zhang et al. [55] pointed out that surface roughness plays a significant role in affecting fracture flow in the review paper. Hence, it is necessary to consider the effect of surface morphology or roughness on fluid flow, which can uncover underlying flow mechanisms inside fractures and advance our understanding of the fracture flow.It can be seen from the previous chapter that, in the nonlinear seepage stage, the forward and reverse seepage characteristics are inconsistent, and whether this is related to the forward and reverse roughness will be briefly discussed in this section. Belem et al. [56], Zhang et al. [57], Grasselli et al. [58], Grasselli [59], and Tatone and Grasselli [60] have proposed methods and formulas for describing directional roughness. In this paper, the directional roughness presented by Tatone and Grasselli [60] is selected as an index to measure the inconsistency of forward and reverse roughness.Grasselli et al. [58] proposed a method and parameter for calculating directional surface roughness based on the 3D scanning point cloud. Tatone and Grasselli [60] based on the method of Grasselli et al. [58] proposed the directional roughness measurement of 2D fracture profile. The method used to calculate the roughness of the 2D fracture profile direction is as follows: 1. Use the least-squares method to establish the best-fitting line and set it as the X-axis. 2. Connect adjacent points to form a series of short straight line segments. 3. Select the analysis direction and calculate each line segment’s inclination angle. For the 2D fracture profile, only two shear directions (forward and reverse) need to be considered. A length score Lθ∗ that is steeper than the gradually higher angle threshold θ∗ can be obtained based on each line segment’s inclination angle. The relationship between the two can be fitted by equation (8). The detailed theory and method can be found in the paper by Tatone and Grasselli [60]:(8)Lθ∗=L0θmax∗−θ∗θmax∗C,where L0 is the length ratio where the angle threshold is 0°, θmax∗ represents the maximum apparent inclination in the analysis direction, and C is a dimensionless fitting parameter, which characterizes the distribution’s shape. By integrating the area under the curve between Lθ∗ and the angle threshold θ∗, the roughness L0θmax∗/C+1 in the analysis direction can be obtained by the following equation:(9)L0∫0θmax∗θmax∗−θ∗θmax∗Cdθ∗=L0θmax∗C+1.Figure10 shows the forward and reverse θ∗ and Lθ∗ scatter points calculated by the Tatone and Grasselli [60] method of the original profile and 8 approximate profiles. The dimensionless parameter C can be determined using equation (8) for the curve fitting of the scatter point. The forward and reverse L0 and θmax∗ can also be obtained through calculation and statistics, and then the directional roughness of the forward and reverse L0θmax∗/C+1 can be determined and are listed in Table 2. Figures 11(a) and 11(b) are, respectively, plotted as curves of the calculated forward and reverse direction roughness L0θmax∗/C+1 and its difference value L0θmax∗/C+1C with the decomposition level. Figure 11 shows that the overall forward roughness L0θmax∗/C+1 is smaller than the reverse roughness. As the number of decomposition levels increases, the difference between the forward roughness and reverse roughness L0θmax∗/C+1C decreases. When it is decomposed to the fourth level, the roughness in the forward and reverse directions have tended to zero, and the subsequent fifth to eighth levels also tend to zero. This is contrary to the non-Darcy coefficient βC curve with the decomposition levels in Figure 7(b). The difference of forward and reverse non-Darcy coefficients βC and the difference of directional roughness L0θmax∗/C+1C under different decomposition levels corresponded. Figure 12 shows a good negative linear relationship between the difference of directional roughness L0θmax∗/C+1C and the difference of non-Darcy coefficient βC between forward and reverse directions. It can be seen from this section that the difference of directional roughness between forward and reverse directions can reflect the difference between forward and reverse non-Darcy flows (the difference of non-Darcy coefficients). However, it should also be seen that there are few samples used in this article, and the exact quantitative relationship between the difference of directional roughness and the difference of non-Darcy flow needs further research.Figure 10
Angle probability distribution and fitting curve of rough fracture profile under different decomposition levels. (a) Original. (b) Level 1. (c) Level 2. (d) Level 3. (e) Level 4. (f) Level 5. (g) Level 6. (h) Level 7. (i) Level 8.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Table 2
Statistical table of relevant parameters of direction roughnessL0θmax∗/C+1 under different decomposition levels.
ForwardReverseL0θmax∗/C+1 (o)L0θmax∗CR2L0θmax∗CR2ForwardReverseS0.45174.5061.9490.9980.54971.0281.9330.98811.39613.297A10.43574.6281.9530.9980.56558.5051.5490.99710.99512.970A20.44671.8461.9760.9980.55451.7931.2960.99710.76812.496A30.46867.5921.9910.9990.53256.0961.5840.99510.57511.547A40.53446.5591.6310.9980.46648.8091.4250.9979.4499.379A50.51543.3633.2330.9930.48534.3272.0860.9985.2765.395A60.52811.0211.0330.9930.47216.8431.7840.9942.8632.855A70.29010.5000.7970.9880.7106.7501.7760.9581.6951.726A80.3885.3180.9720.9120.6123.0610.7710.9681.0461.058Figure 11
The relationship between the roughness in the forward and reverse directions (a)L0θmax∗/C+1 (b) and its difference L0θmax∗/C+1C with the decomposition level.
(a)(b)Figure 12
Relationship between the difference of non-Darcy coefficientβC and the difference of directional roughness L0θmax∗/C+1C.


## 6.2. Single Rough Fracture in Fracture Network

The single fracture is the basic unit of the fracture network, and the seepage characteristics of the single fracture will also affect the seepage characteristics of the fracture network. The research on the difference between the forward and reverse seepages in single rough fracture can lay the foundation for future research on the forward and reverse seepage characteristics of the rough fracture network. When it is considered that the single fractures in the fracture network obey the cubic law, the seepage characteristics of the fracture network can be described by the anisotropic Darcy’s law, as shown in the following equation:(10)uj=−KjiJi,where uj is the velocity in direction j, Ji is the hydraulic gradient in direction i, and Kji is the equivalent permeability coefficient tensor. As shown in this paper’s research results, in the case of Darcy seepage in rough single fractures, there is no difference in the seepage characteristics of the forward and reverse seepage directions, and both strictly obey the cubic law. In the Darcy seepage state, the permeability coefficient K of the fracture network composed of rough single fractures in the relative seepage direction (such as 0° and 180°, 90° and 270°) is also consistent, which conforms to the characteristics of a tensor. Thus, the equivalent permeability coefficient tensor Kij can express the permeability characteristics of the seepage in the Darcy seepage state.With the increase of hydraulic gradient, not only the single fracture but also the fracture network will enter non-Darcy seepage. For the fracture network, Liu et al. [42], Liu et al. [61], and Yin et al. [62] point out that when the hydraulic gradient is greater than a certain value, the ratio of flow velocity to hydraulic gradient is no longer constant, which no longer satisfies Darcy’s law. The seepage flow in the fracture network will also enter the state of non-Darcy seepage. Similar to the use of Forchheimer law (equation (6)) to describe the non-Darcy flow in a single rough fracture, the non-Darcy flow in a fracture network is also described using the Forchheimer law (equation (11)) of anisotropic seepage [63]:(11)−Ji=Aijuj+Bijuuj=Kij∗Kuj+βijguuj.In the equation,Aij and Bij are, respectively, the coefficient tensors of the viscous force term and the inertial force term, which can also be expressed as a matrix related to the permeability coefficient tensor Kij and the non-Darcy coefficient tensor βij,u is the modulus of flow velocity, Kij∗/K is the inverse matrix of the equivalent permeability coefficient matrix, and βij is the equivalent non-Darcy coefficient tensor.Since the forward and reverse seepage characteristics of a single smooth fracture will not differ between Darcy and non-Darcy seepage conditions, for a fracture network composed of smooth single fractures, the fracture network’s seepage characteristics under forward and reverse seepage will not be different. In other words, the permeability coefficientK of the smooth fracture network will remain consistent in any state under the relative seepage direction (forward and reverse direction). For smooth fracture networks, the equivalent non-Darcy coefficient β in the relative seepage direction obtained by equation (11) will also be consistent. The equivalent non-Darcy coefficient β conforms to the tensor’s symmetry, which can theoretically constitute the equivalent non-Darcy coefficient tensor βij [64]. However, this paper’s research shows that when Re > 10, the fracture seepage of a single rough fracture enters the non-Darcy seepage stage, and then the difference of forward and reverse seepage occurs. The quantitative description by th Forchheimer law (equation (9)) shows that, in the state of non-Darcy seepage, the permeability coefficient k has nothing to do with the direction of seepage, and the non-Darcy coefficient β is related to the direction of seepage. Therefore, for the fracture network composed of single rough fractures, the difference of the forward and reverse flow direction of the single rough fracture will also lead to the difference in the rough fracture network’s relative flow direction. If equation (11) is used, the equivalent non-Darcy coefficient β in the relative seepage direction will no longer be consistent and will no longer conform to the tensor’s symmetry. So, for rough fracture networks, the applicability of equation (11) will be limited. This will also raise new problems and challenges for describing the rough fracture network’s non-Darcy seepage characteristics accurately and quantitatively.In the actual natural environment and engineering environment, the single fractures that make up the fracture network in the rock mass are rough. In the non-Darcy seepage state of these rough single fractures, the forward and reverse seepage directions’ seepage characteristics will be different. How much do these differences affect the difference in the direction of seepage flow in the fracture network? In turn, will the difference in the flow direction of the rough fracture network for seepage flow and solute movement have much impact on the project? These issues need to be addressed by an in-depth research in the follow-up work.


## 7. Summary and Conclusion

This paper’s primary purpose is to investigate the discrepancy between forward and reverse seepage characteristics of single rough fractures. The wavelet transform technique is used to decompose the fracture profile by multiple scales (multidecomposition levels). The approximate models of the detailed profiles with different frequencies (different decomposition levels) are obtained. The original profile model and eight approximate models are then numerically simulated with the Reynolds number varying from 0.001 to 1000 utilizing FVM (Finite Volume Method). Through the above research, the following conclusions can be drawn:(1)
The necessary condition for the various hydraulic characteristics of two directions in single rough fractures is nonlinear seepage. When Re < 10, the correlation between the pressure gradient and flow rate in this paper strictly satisfies the cubic law in all fracture models, and there is no difference between seepages in two directions. When Re > 10, the correlation between the pressure gradient and flow rate in this paper no longer strictly satisfies the cubic law in all fracture models. There is a discrepancy between the forward and reverse seepages only from the original fracture model to the 4th approximate fracture model. Still, there is no difference between seepages in two directions from the 5th to 8th approximate fracture model.(2)
The discrepancy between the forward and reverse hydraulic seepage characteristics of the single rough fracture is the different shapes and distribution of the eddy currents generated by the forward and reverse nonlinear seepage. In the nonlinear seepage state, from the original fracture model to the 3rd approximate fracture model, eddy currents’ shape and location are different in the seepage of two directions. The region of eddy currents in both directions gradually narrows with the increase of the decomposition level. So, the discrepancy of nonlinear seepage characteristics in two directions is also gradually reduced.(3)
The secondary roughness provides space for the generation and development of eddy currents. Simultaneously, it provides a geometric basis for the discrepancy of hydraulic characteristics in two directions. In this paper, the 4th approximate profile is taken as the primary roughness. The sum of the first four detailed profiles (D1 + D2 + D3 + D4) is regarded as the secondary roughness. The difference between the effective apertures in the forward and reverse seepages, the difference between non-Darcy coefficients of two directions, and the comparison of local eddy currents show that the 4th approximate profile is the boundary point of whether there will be a discrepancy in the forward and reverse hydraulic characteristics.In this paper, the difference between the forward and reverse nonlinear seepage flows under the single rough fracture’s multiscale roughness was analyzed in detail. We adopt the fracture model whose upper and lower fracture profiles are the same, without considering different upper and lower fracture profiles in nature. Only a 2D fracture model is discussed in this paper. More complex 3D fracture models still need to be further studied. We will gradually carry out the above-unfinished researches in the succeeding stage.



---

*Source: 6648522-2021-03-22.xml*