Numerical Calculation of Fracture Seepage in Rough Rock and Analysis of Local Pressure Drop

The seepage performance of a rock mass mainly depends on the rock fractures developed in it. Numerical calculation method is a common method to study the permeability properties of fractures. Seepage in rock fractures is affected by various factors such as fracture aperture, roughness, and filling, among which aperture and roughness are the two most widely influenced factors. The Navier-Stokes (NS) equation can be solved directly for the seepage flow in rock fractures with good accuracy, but there are problems of large computational volume and slow solution speed. In this paper, the fracture aperture space data is substituted into the local cubic law as an aperture function to form a numerical calculation method for seepage in rough rock fractures, namely, the aperture function method (AFM). Comparing with the physical seepage experiments of rock fractures, the calculation results of AFM will produce a small amount of error under the low Reynolds number condition, but it can greatly improve the calculation efficiency. The high efficiency of calculation makes it possible to apply AFM to the calculation of largescale 3D rough fracture network models. The pressure drop of fluid in the fracture has viscous pressure drop (VPD) and local pressure drop (LPD). VPD can be calculated using the AFM. After analyzing the results of solving the NS equation for fracture seepage, it is concluded that the LPD includes the pressure drop caused by area crowding in the recirculation zone (RZ), kinetic energy loss in the RZ, kinetic energy loss in the vortices, and other reasons.


Introduction
Rock fractures are widely distributed in subsurface rock masses, and for fluid transport, the permeability of rock fractures is significantly greater than that of the rock matrix [1]. The hydraulic properties of rock fractures play an important role in assessing the performance of subsurface engineering, such as geothermal energy development, enhanced recovery, nuclear waste disposal, and groundwater pollution control [2][3][4][5][6]. Fluid flow in rock fractures can be calculated using the Navier-Stokes (NS) equation. The NS equation is used to solve the fracture flow with high accuracy, but the fracture aperture dimension is much smaller compared to the dimension in the extension direction, which makes the modeling difficult and the calculation workload is large. Zhou et al. [7] used the NS equation to solve a real rock fracture seepage model with a size of 150 mm × 120 mm, and the number of model elements exceeded 10 6 to ensure the computational accuracy.
Due to the complexity of the NS equation solution, a concise fracture seepage equation, the cubic law, can be derived by neglecting the nonlinear terms in fluid flow. The same conclusion was obtained for the seepage experiments using two smooth flat plates [8]. The natural rock fracture surface is rough, which restricts the flow of fluid in the fracture, and using the cubic law to calculate the fracture flow would overestimate the permeability of the fracture [9,10]. Through theoretical derivation and experimental research, some scholars proposed that under the condition of low Reynolds number (Re < <1), the effect of inertial flow of fluid can be neglected, and the local cubic law (Reynolds equation) can be used to calculate the fracture seepage flow, and its calculation results do not deviate much [11][12][13][14]. Zimmerman and Bodvarsson [12] found by single fracture seepage experiment that when Re is in the range of 1~10, there is a weak inertia effect zone within the fluid in the fracture, and the fracture seepage flow gradually deviates from the linear relationship with the hydraulic gradient but becomes a nonlinear relationship. By correcting the geometric mean aperture, arithmetic mean aperture, surface roughness factor, and curvature factor, the local cubic law can still be used to evaluate the fluid flow in the fracture [11,14,15]. Zhu et al. [16] concluded by theoretical derivation and experimental study that the cubic law still has good applicability in the laminar flow range (Re ≤ 2300) when the flow of the fracture is long enough.
Javadi et al. [17] and Liu et al. [18] found that there is a deviation between the calculated results of linear and nonlinear flow in the fracture, but the difference is not too big at low Reynolds number conditions. Javadi et al. [17] considered that there are two types of pressure drop during fluid flow in the fracture: one is viscous pressure drop (VPD) due to fluid viscosity; the other is the local pressure drop (LPD) due to the sudden change of aperture. For laminar flow in a smooth flat plate, the VPD can be calculated using the cubic law. For rough fractures, corrections for fracture wall roughness and fracture curvature are required.
In this paper, the aperture function method (AFM) based on the local cubic law is introduced, which utilizes the fracture aperture function as the parameter of the numerical model, where one fracture wall surface is used as the geometric model, and the rough fracture seepage is calculated based on the local cubic law. The results of the AFM calculations are compared with the physical experiments and the results of the numerical model of rough fracture seepage based on the solution of the NS equation. Multiple pressure boundary conditions are applied to the inlet of the fracture models, and the rough fracture seepage flow is solved based on the NS equation, and the formation of the recirculation zone (RZ) in the fracture flow and the reasons for the influence on the seepage flow are analyzed in detail. The rotational motion of fluid flow is called vortex, and the presence of vortices has a significant effect on the motion of the fluid due to the vortices. This paper presents the mathematical definition and identification of vortices in fluids and identifies the vortices and their distribution in rough fracture flows. A comparative analysis of vortices in fracture flow and RZ is carried out, and the effect of vortices on fracture seepage flow is studied.

Methodology
2.1. Aperture Function Method (AFM) 2.1.1. Aperture Interpolation Method. Rough rock fracture surface morphology is complex and can be obtained by using high precision laser scanning technology [19][20][21], and the obtained data can well reflect the real fracture surface morphology. For the obtained fracture surface data, let the fracture surface undulation height be z and use x and y to represent the fracture plane coordinates; then, the fracture surface geometry can be represented by the function z = gðx, yÞ in the Cartesian coordinate system. The 3D scanned fracture surface is processed to form a list of values that can be expressed by the function z = gðx, yÞ. For subsequent calculations, the x-direction is set to move in the shear direction. The upper fracture surface function is g U ðx, yÞ, and the lower fracture surface function is g L ðx, yÞ. To calculate the aperture between the two fracture walls, the following equation can be used [22][23][24].
where eðx, yÞ is the upper and lower fracture wall aperture functions, u s is the shear direction displacement, and u v is the normal displacement. When eðx, yÞ does not correspond exactly to the spatial coordinates ðx, yÞ of g U ðx, yÞ (or g L ðx, yÞ), it is necessary to use the interpolation method to find the aperture value at the desired spatial location. In Figure 1(a), the g i,j point on the geometric slit surface corresponds to the e i,j ′ point in the aperture, but there is no aperture value at this point, so the interpolation method is used to obtain the aperture value at the e i,j ′ point.
The fracture aperture is a function of the spatial coordinates ðx, yÞ, so each of the four points can be expressed as ðx 1 , y 1 , eðx 1 , y 1 ÞÞ, ðx 1 , y 2 , eðx 2 , y 1 ÞÞ, ðx 2 , y 1 , eðx 2 , y 1 ÞÞ, and ðx 2 , y 2 , eðx 2 , y 2 ÞÞ (Figure 1(b)). If we want to get the aperture value eðx, yÞ at a certain position in the middle of the four points, we can use linear interpolation to estimate it: First, interpolate x to find eðx, y 1 Þ and eðx, y 2 Þ: Then, use the same idea to find eðx, yÞ: The above equation gives the value of the aperture at any position within the spatial distribution of the aperture. Each point g i,j ðx, yÞ on the fracture wall grid can also be found corresponding to the desired aperture value.
When discretizing the fracture numerical model, each fracture wall grid point corresponds to a tiny geometric grid cell. Each grid cell can find the corresponding fracture aperture data, and within this cell, the aperture value is a given value; then, the geometric grid cell and the aperture value can form a tiny virtual flat flow model. For each tiny virtual flat plate model, the fracture flow calculation can be performed using the cubic law, and the overall view is an 2 Geofluids application of the local cubic law. Since the fracture surface and fracture aperture values of each cell are a list of values that vary with space, a rough fracture model with two walls is formed. A triangular grid Figure 2(a) or a rectangular grid Figure 2(b) can be used when the fracture data are discrete.
Substituting the aperture function eðx, yÞ into Equation (4) yields a clearer expression of the equation [27].

∂ ∂x
where h is the hydraulic head, μ is the dynamic viscosity, ρ is the fluid density, and g is the gravitational acceleration, and eðx, yÞ is the aperture function. Figure 2 shows the meaning of Equation (5) graphically by dividing the rough fracture surface into a series of connected small parallel plates. In this paper, the calculation method of using the aperture function to characterize the aperture value and solving the rough fracture flow based on the local cubic law is referred to as the aperture function method (AFM).

Recirculation Zone in Fracture Seepage
Flow. Lee et al. [25], Zhou et al. [28], and Zhou et al. [7] found through physical and numerical experiments that in the steady laminar flow state, the fluid in the fracture will generate a recirculation zone (RZ) when flowing through the abrupt change of aperture. The fluid in the RZ also consumes part of the flow energy, and because the fluid in the RZ occurs back on its own, it does not produce effective flow at the out-let of the fracture. The size of the RZ is nonlinearly related to the fluid flow rate and the fracture roughness procedure, and its magnitude cannot be calculated using a simple equation. In addition, Lee et al. [29] found experimentally that under specific conditions, fluid can also produce slip phenomenon at the fracture wall, which also affects the overall seepage flow rate of the fracture.
In order to further verify whether a fluid RZ is generated in the fracture model in this paper and to analyze whether the RZ effect is a linear or nonlinear problem, with other geometric and aperture conditions unchanged, a variety of pressure boundary conditions are set, and two methods, NS solution and AFM solution, are used to calculate and compare, respectively, ( Table 1). The JRC1 fracture model inlets were applied with pressures ranging from 10 Pa to 6000 Pa (−dh/dl = 0:01~5:1) for a total of 52 submodels. JRC3 fracture model inlets were applied with pressures ranging from 10 Pa to 890 Pa (−dh/dl = 0:01~0:76) for a total of 36 submodels. JRC6 fracture model inlets were applied 25 Pa~650 Pa pressure (−dh/dl = 0:02~0:55), respectively, for a total of 26 sub-models. JRC9 fracture model inlets were applied with 10 Pa~1000 Pa pressure (-dh/dl =0.01~0.85), respectively, for a total of 34 submodels. The outlet pressure of all fracture models is 0 Pa. The fluid density in the model is 998.2 kg/m 3 , and the dynamic viscosity coefficient is 0.0010093 Pa.s. The calculation software was carried out using COMSOL Multiphysics finite element method.

Identification of Vortices in Seepage
Flow. The phenomenon of vortices is prevalent in the fluid flow process. Intuitively, vortices represent the rotational motion of the fluid; that is, where there is fluid rotation, there are vortices, and conversely, where there are vortices, there is fluid rotation [30]. In order to give a quantitative mathematical definition of vortex, a parameter Ω is introduced, representing the ratio of the size of the rotating part of the vortex to the size of the total vortex, calculated as [31] x y z e i+1,j+1 x y e (x 1 , y 1 ) e (x 2 , y 1 ) where k k F represents the Frobenius parametrization of the matrix, A is the symmetric part of the fluid velocity gradient tensor, and B is the antisymmetric tensor (vorticity). According to its physical meaning, it is obvious that Ω takes values in the range 0 ≤ Ω ≤ 1, which can be interpreted as the concentration of vorticity, and more specifically, Ω represents the rigidity of fluid motion, and when Ω = 1, it represents the fluid doing rigid body rotation. When Ω > 0:5 represents that the antisymmetric tensor B dominates over the symmetric tensor A. Therefore, Ω slightly larger than 0.5 can be used as the criterion for vortex identification. Liu et al. [31] proposes to use Ω = 0:52 to determine the boundary of vortex.

Analysis of AFM Calculation Results.
In order to verify the rationality of the AFM, this paper compares the physical experimental results of rough fracture fluid flow completed by Zhu et al. [16,32]. The typical definition of rock fracture roughness curve (JRC) in rock mechanics was adopted for the experiment, and the 1st, 3rd, 6th, and 9th of these curves were selected and combined with a smooth flat plate to form four fracture models with different roughness, each with a fracture length of 100 mm and a minimum aperture value of 0.51 mm. The four models were numbered as JRC1, JRC3, JRC6, and JRC9, respectively. Multiple experiments were conducted for each artificial fracture sample, and the pressure gradient and flow rate were averaged.
Based on the AFM, the numerical model of the four fractures was completed and calculated using the finite element

Geofluids
The error values in Table 2 are calculated as where ε represents the flow error, q s represents the experimentally measured flow rate per unit width of the fracture, and q c represents q p calculated according to the average aperture and the cube law or q n calculated using the AFM. Table 2 shows the comparison of different calculation methods with experimental results. It can be seen that the results calculated by the average aperture and the overall cubic law deviate very much from the experimental measured results. When the fracture roughness increases, the error of the overall cubic law calculation results is even larger. The results calculated by the numerical model using the AFM are in good agreement.

Comparison of the Results of AFM and NS Method.
The results of the flow calculated by the AFM are compared with the results of the NS equation solution. We plotted the flow rate variation curves with hydraulic gradient derived from the two calculation methods and plotted the error curves of the two calculation methods (Figure 4).
The error curve shows that the AFM calculation result overestimates the permeability of the fracture, causing its calculated flow rate to be large, which is consistent with the results of Yeo et al. [9] and Bauget and Fourar's study [10]. The errors of NS and AFM calculations are small at low hydraulic gradient and become larger as the hydraulic gradient gradually increases, which is consistent with the results of Oron and Berkowitz [11] and Zimmerman et al. [13]. The increase of fracture roughness also increases the flow error (the overall error of JRC9 model is the largest).
The flow error variation also shows a nonlinear variation pattern, rather than a simple linear increase.  Table 2 shows that the results calculated by the AFM under the low hydraulic gradient conditions are in good agreement, but there is still a certain degree of error. The theory of fracture seepage calculation used in the AFM is the local cubic law, which makes a linear relationship between flow and pressure drop (corresponding to the hydraulic gradient) because it ignores the effect of fluid inertia forces. In contrast, Javadi et al. [17] considers the pressure drop during the flow of fluid in a fracture as VPD and LPD. VPD is due to fluid viscosity and can be calculated using the cubic law.
LPD is the pressure drop (or energy loss) that occurs when the fluid flows through a part of the fracture where the aperture changes abruptly. When using the local cubic law-based AFM for seepage calculations, the full LPD cannot be calculated, which is the main reason for the deviation of the AFM calculation results. Zhu et al. [16] and Liu et al. [18] also found the phenomenon of LPD in their studies. The RZ in the fracture flow is one of the important reasons for the LPD. When using the NS solution method for fracture seepage calculations, the extent of the RZ can be determined in the fracture model using streamlines. In this paper, we take JRC3 and JRC6 models as examples to discuss the evolution of RZ and the influence of RZs on fracture seepage.  (Figures 5(c) and 6(c)). With the increase of the applied hydraulic gradient, the influence of fluid inertial force gradually increases, and small RZs appear in the region of large changes in the fracture wall (Figures 5(f) and 6(f)), and the range of RZs is small and has little influence on the mainstream zone. The further increase of hydraulic gradient, the influence of fluid inertia force is further revealed, not only to expand the range of RZ but also the emergence of a number of RZs of different sizes, resulting in the narrowing of the mainstream zone, bending   Figures 7 and 8 shows the variation of the fluid velocity profile with hydraulic gradients on a cut line in the JRC3 and JRC6 fractures. When the hydraulic gradient is small, the velocity distribution on the intercept line has a parabolic shape; when the hydraulic gradient continues to increase, the velocity profile has a bimodal shape due to the creation of the RZ, reflecting the formation of the RZ.
The most classical mathematical description of nonlinear flow of a fluid in a fracture is Forchheimer's law [13,33] where ∇p is the hydraulic gradient ∇p = ρg ðdh/dlÞ, v is the average flow velocity of the fluid under hydraulic gradient conditions v = q/e h , and k is the intrinsic permeability of the fracture k = e h 2 /12. β½m −1 is the inertial resistance coefficient of the fluid, e h is the equivalent hydraulic aperture of the fracture, μ is the viscosity coefficient of the fluid, and ρ is the fluid density. It should be noted that k is the permeability when the hydraulic gradient (−dh/dl) is tiny and the effect of fluid inertial forces can be neglected. When the hydraulic gradient increases, e h changes when there is the creation of a RZ. The kðe h Þ and the coefficient β can be obtained from the applied hydraulic gradient and the flow fitting (Equation (9)). The fitted equations of the Forchheimer equation for the JRC3 and JRC6 models are labeled in Figure 9, respec-tively, and the fitted curves are perfect, reflecting the phenomenon that the nonlinear flow in the fracture becomes progressively more pronounced with the increase of the hydraulic gradient.
Reynolds number is a quantified parameter to compare the proportional relationship between inertial and viscous forces in the fluid, and for the fluid passing through the fracture, the Reynolds number is calculated as [13]: where v is the average flow velocity of the fluid in the fracture, μ is the viscous coefficient of the fluid, ρ is the fluid density, and D is the characteristic length in the flow system; here, the average aperture of the fracture is taken. Q is the fracture flow rate, and w is the width of the fracutre sample perpendicular to the direction of the hydraulic gradient. We calculated the Reynolds number for each hydraulic gradient condition in the JRC3 and JRC6 fracture NS seepage model, respectively, and the variation of the Reynolds number is represented as a curve (blue solid line in Figure 9). The obtained Reynolds number can only indicate the trend of increasing fluid inertia force with the hydraulic gradient, but it does not tell the kind of conditions in which significant nonlinear flow occurs.
In order to quantify the nonlinear flow, Javadi et al. [34] and Zhou et al. [33] proposed to use the critical Reynolds number (Re c ) to make a judgment. Re c indicates the beginning of the transition from fluid flow to non-Darcy flow and can be defined as the Reynolds number when the percentage of nonlinear pressure drop (βρv 2 ) to total pressure drop (μv/k + βρv 2 ) reaches the critical point α. Re c can be obtained from the following equation.
where α = 5% and e h , β, and k are obtained from the Forchheimer's law fitting equations for the two models JRC3 and JRC6. Using Equation (11), we obtained that Re c = 109:27 (Figure 9(a)) for JRC3, corresponding to a hydraulic gradient −dh/dl = 0:30 (Table 3) and Re c = 74:03 (Figure 9(b)) for JRC6, corresponding to a hydraulic gradient −dh/dl = 0:16 (Table 3). Referring to the meaning of the critical Reynolds number, it can be considered that the hydraulic gradient corresponding to the critical Reynolds number, theoretically the fluid inertia force for the overall pressure drop is equal to 5%; combined with the Equation (10) to calculate the variation of Reynolds number, the smaller the hydraulic gradient, the smaller the Reynolds number, the impact of the inertia force is also smaller.

Geofluids
RZ or only a small RZ is formed, and the linear flow is dominant at this time. And when greater than the critical Reynolds number (Re c ), the formation of obvious RZ, and with the increase of hydraulic gradient and increase, the nonlinear flow in the fracture flow is more significant.

The Effect of Recirculation Zone.
In order to further study the influence of RZ on fracture seepage, we analyze the area occupied by the RZ and the kinetic energy of fluid in the RZ, respectively. The extent of the RZ is difficult to be directly defined, but there must be a velocity vector in 8 Geofluids the RZ that is opposite to the main seepage direction. The main seepage direction of all models in this paper is the x -direction, i.e., the range of velocity u < 0 in the x-direction is part of the RZ. According to the interrelationship between the velocity vector and the flow line in Figure 10, it is found that the range of u < 0 can be approximated to half of the area of the RZ. The area of u < 0 under a certain hydraulic gradient is calculated and taken as two times its size as the area of the RZ, although there is a certain deviation from the real RZ area, but it can reflect the RZ area change law. The fluid in the RZ for circular flow, its internal fluid kinetic energy is calculated using a similar method. Calculate the kinetic energy of the fluid with velocity u < 0 in the x -direction and then multiply it by 2 as the total kinetic energy of the fluid in the RZ. Then, analyze the law of kinetic energy change in the RZ. The kinetic energy of the fluid is calculated as follows: where E is kinetic energy, Ω is the area of RZs, ρ is fluid density, U is fluid flow velocity (scalar), and ΔV is element volume.
The main purpose of this paper is to analyze the proportion of the RZ area to the total fracture area and the change of the proportion of the kinetic energy of the RZs to the total kinetic energy of the fluid. The influence of the RZs is analyzed according to the proportional change curve. The calculated results of solving the NS equation for four models, including JRC1, JRC3, JRC6, and JRC9, are counted to show the area of the RZs, the ratio of kinetic energy to total area, and total kinetic energy under each hydraulic gradient, and the change curves are plotted.
The JRC3 model is used as an example for illustration ( Figure 11(a)). The hydraulic gradient of JRC3 fracture model ranges from 0.009 to 2.55, and the ratio of the RZ area to the total area is 0.06% to 12.02%. With the increase of the hydraulic gradient, the percentage change of the RZs area is very obvious. The proportion of the kinetic energy of the RZ area to the total kinetic energy is 0.0002%~0.21%. The kinetic energy percentage is much lower compared to the area, and under the condition of small hydraulic gradient, the effect of kinetic energy in the RZs can even be ignored. The remaining three models JRC1, JRC6, and JRC9 have similar characteristics (Figures 11(b)-11(d)).
Therefore, it can be considered that the main effect of the RZs on the rough fracture seepage is that the fluid occurs back in the RZ by itself and does not produce an effective flow at the fracture outlet. The fluid flow in the RZs is slow, and the loss of kinetic energy of the fluid in the fracture is not significant. However, the RZs squeeze the seepage channel of the fracture, resulting in the reduction of the actual seepage capacity of the fracture. The area occupied by the RZ increases with the hydraulic gradient and shows a nonlinear variation.
In addition, the area of RZs and the proportion of kinetic energy do not increase infinitely with the increase of 9 Geofluids hydraulic gradient, but gradually reach a peak, and then stabilize or slightly decrease. And the rougher the fracture is, the easier it is to reach the peak. In terms of the peak area share of the RZ, the four models JRC1, JRC3, JRC6, and JRC9 reach the peak at the hydraulic gradient of 5.1, 1.68, 0.76 (presumed), and 0.57, respectively.
3.4. Analysis of the Effect of Vortices. The RZ is circled according to the shape of the streamline, and its shape is influenced by the number of streamlines and the calculation method of the flow line. The RZ is intuitively seen as a backflow or rotation of the fluid, which is similar to the physical concept of a vortex. However, the vortex determined according to a strict mathematical definition differs significantly from the RZ (Figure 12). Vortices are generated near the fracture wall, inside the RZ, at the transition between the reflux and main flow zones, or even within the main flow zone. That is, when there is a velocity difference within the fluid, vortices may be generated. Figure 12 shows that the range of vortices also includes part of the mainstream zone; that is, there is also local fluid rotation within the mainstream zone. This part of the vortex is different from the RZ; it will lose part of the fluid energy but at the same time can produce effective flow.
To further analyze the effect of vortices on the seepage in the fracture, the value of the vortex parameter Ω in the fluid within the fracture is calculated using Equation (6), and the kinetic energy of the fluid in the range Ω ≥ 0:52 is calculated    11 Geofluids using Equation (12), along with its percentage of the total kinetic energy. Then, the variation of total kinetic energy, intravortex kinetic energy, and the proportion of intravortex kinetic energy of the fluid is shown as curves using the hydraulic gradient as the horizontal coordinate. Figure 13 shows the curves for the four fracture models. It can be seen that (a) the intravortex kinetic energy increases nonlinearly with the hydraulic gradient. (b) The variation of intravortex kinetic energy ratio is much more complicated. When the fluid Reynolds number is smaller than the critical Reynolds number (Re c ), the intravortex kinetic energy ratio shows a decreasing trend. (c) The proportion of kinetic energy in the vortex first decreases, then rises, and then the curve changes to a smooth, indicating that the kinetic energy in the vortex is not completely lost energy, and the fluid in the vortex both rotational motion, but also with the mainstream of the flat motion. (d) Compared with the proportion of kinetic energy in the RZ, the proportion of kinetic energy in the vortex is significantly larger. In the JRC1 model, for example, the highest is 22.16%, and the lowest is 14.85%, so it is obvious that the proportion of kinetic energy in the vortex is much larger and has a greater impact.

Discussion
For the rough fracture seepage model, it is obvious that solving the NS equation directly can reflect the nonlinear flow state in the fracture more accurately; for example, the fit of Forchheimer equation is very high, and the corresponding calculation results are more accurate. However, its disadvantage is also very obvious; that is, the computational effort is very large. In comparison, the computational efficiency of AFM is much higher than the former.
Comparing the four models, the number of discrete elements of the model is much larger than that of AFM when solving the NS equation directly (Table 1). This is because AFM only needs to discrete one fracture wall surface. The former has 8.76 times, 4.90 times, 3.41 times, and 3.3 times more than the latter (Table 1).
Numerical calculations were performed using the same computer (CPU: Intel(R) Xeon(R) CPU E5-1620 v3, RAM: 16G); the difference in computational efficiency between the two is huge. The JRC1, JRC3, JRC6, and JRC9 models took 169415 s, 176712 s, 58665 s, and 14772 s, respectively, to solve the NS equations (Table 3), while the computation time using AFM was 19 s, 28 s, 8 s, and 9 s, respectively (Table 3).

Geofluids
The huge computational efficiency advantage of AFM can complete more model computation work in a short time. Because the calculation time-consuming of a single rough fracture model is very short, the seepage calculation of a complex three-dimensional spatial rough fracture network (e.g., a fracture network composed of thousands of rough fractures) can be completed using only an ordinary computer, which is one of our future works. In addition, the geometric model of AFM requires only one fracture wall surface, and the modeling process of numerical model is simpler, while the aperture data can be produced more quickly with the help of third-party software or programming.
From this point of view, the huge numerical calculation efficiency of AFM can offset the impact of its calculation error, and the error is not significant at low Reynolds number conditions. This will help researchers to carry out seepage calculations and research work on large scale rock fracture networks.

Conclusions
From the above analysis and description, we can get several conclusions as follows.
(1) Under the low Reynolds number condition, the numerical calculation work of seepage flow in rough fractures can be done by using the AFM, but the calculation results will produce some small errors (2) Due to the abrupt change of the fracture wall to produce recirculation zone (RZ), the main effect of RZ is that it crowds the main flow range, making the effective seepage channel of the fracture narrower. The area of the RZs is related to the hydraulic gradient and varies nonlinearly. In addition, the fluid motion inside the RZ is slow, and the loss of kinetic energy of the fluid is not large, accounting for a small proportion of the overall kinetic energy (3) In addition to the RZ in the rough fracture, vortices are formed, i.e., local fluid rotation. Vortex and RZ area are very different, the fluid in the vortex loss of kinetic energy, but can produce effective flow, and its changes are more complex (4) Using the local cubic law to solve the rough fracture seepage flow, the main source of error in the calculation results is the LPD. The LPD in fracture seepage includes the area crowding in the RZ, kinetic energy loss in the RZ, kinetic energy loss in the vortex, and pressure drop caused by other reasons (5) The main advantage of AFM is the efficiency of its numerical calculation. Based on this advantage, AFM can quickly complete the numerical calculation of a huge number of 3D rough fracture networks Strictly speaking, the fracture model in this paper is a single-wall rough fracture model, which has some shortcomings of its own; for example, the nonuniformity and anisotropy of the fracture are not sufficiently reflected. This is mainly due to the need to be consistent with the fracture physical model for comparative analysis. These will be improved gradually in our subsequent studies.

Conflicts of Interest
The authors declare that they have no conflicts of interest.