Seepage Characteristics Study of Single Rough Fracture Based on Numerical Simulation

Featured Application: The application of research results is mainly in the ﬁelds of rock fracture seepage, geothermal development and utilization, etc. Abstract: Fracture seepage is an important aspect of groundwater research, but due to the closure of fractures and the randomness of wall surface roughness, it is a challenge to carry out relevant research. Numerical simulation serves as a good way to solve this problem. As such, the water ﬂow in single fracture with different shapes and densities of roughness elements (various bulges/pits on fracture wall surfaces) on wall surface was simulated by Fluent software. The results show that, in wider rough fractures, the ﬂow rate mainly depends on fracture aperture, while, in narrow and close rough fracture medium, the surface roughness of fracture wall is the main factor of head loss of seepage; there is a negative power exponential relation between the hydraulic gradient index and the average fracture aperture, i.e., with increase of rough fracture aperture, both the relative roughness of fracture and the inﬂuence of hydraulic gradient decrease; in symmetrical-uncoupled rough fractures, there is a super-cubic relation between the discharge per unit width and average aperture; the rough fracture permeability coefﬁcient K is not a constant which is affected by the scale effect and the density of roughness elements. Results found provide further understanding of rough fracture seepage.


Introduction
The bedrock fracture seepage characteristics are related to many research fields, such as groundwater utilization, underground engineering, hydraulic engineering, geothermal development, etc. [1][2][3][4]. Single fracture is the most basic component unit of fracture media, so the flow characteristics of single fracture is the key point to study the seepage and solute transport of fracture media [5].
The early research on bedrock fracture seepage mostly relied on Darcy's law, which characterized pore seepage, but the application of Darcy's law had preconditions, including media homogeneity and isotropy, temperature and pressure. Studies had shown that Darcy's law was only valid within a certain hydraulic gradient; if this range was not reached or exceeded, the linear relation described by Darcy's law no longer existed, and Darcy's law was no longer applicable [6,7].
Snow put forward the LCL (local cubic law) through simulation experiment of water flow in smooth parallel plates and found that the discharge through fracture cross-section 2 of 16 was proportional to third power of the fracture aperture [8]. LCL was based on a smooth parallel plates model, and an important parameter, the roughness of fracture wall surface, was not taken into consideration. There is no ideal smooth fracture in the nature, and the existence of roughness elements on fracture wall surface will narrow down the passage of water ( Figure 1). Therefore, in practical applications, LCL often overestimated the seepage capacity of rough fracture [9,10]. The real fracture surface is generally rough and uneven, and the aperture also changes randomly at different positions in fracture, so the concepts of hydraulic fracture aperture e h , average fracture aperture e, and mechanical fracture aperture e m were introduced. Hydraulic fracture aperture e h was also called equivalent hydraulic fracture aperture, which was an inference value, i.e., the rough fracture was equivalent to a smooth parallel plates fracture with aperture e h , and its value could be calculated according to LCL after the flow of the rough fracture was measured; the average aperture e was the average value of the fracture apertures along the fracture; the mechanical fracture aperture e m referred to the maximum closing deformation value of the two surfaces of fracture. Appl. Sci. 2022, 12, x FOR PEER REVIEW 2 of 16 Snow put forward the LCL (local cubic law) through simulation experiment of water flow in smooth parallel plates and found that the discharge through fracture cross-section was proportional to third power of the fracture aperture [8]. LCL was based on a smooth parallel plates model, and an important parameter, the roughness of fracture wall surface, was not taken into consideration. There is no ideal smooth fracture in the nature, and the existence of roughness elements on fracture wall surface will narrow down the passage of water ( Figure 1). Therefore, in practical applications, LCL often overestimated the seepage capacity of rough fracture [9,10]. The real fracture surface is generally rough and uneven, and the aperture also changes randomly at different positions in fracture, so the concepts of hydraulic fracture aperture eh, average fracture aperture e , and mechanical fracture aperture em were introduced. Hydraulic fracture aperture eh was also called equivalent hydraulic fracture aperture, which was an inference value, i.e., the rough fracture was equivalent to a smooth parallel plates fracture with aperture eh, and its value could be calculated according to LCL after the flow of the rough fracture was measured; the average aperture e was the average value of the fracture apertures along the fracture; the mechanical fracture aperture em referred to the maximum closing deformation value of the two surfaces of fracture. Other researchers had performed research work before and after the formal establishment of LCL theory, and proposed and developed the relation among discharge per unit width q, fracture representative aperture ẽ (including hydraulic fracture aperture eh, average aperture e , and mechanical aperture em), and fracture roughness (including relative roughness and absolute roughness), which modified and supplemented LCL theory (see Section 2.1). Xu et al. found that, when the average fracture aperture e was used to replace the hydraulic fracture eh, there was a super-cubic or sub-cubic relation between the discharge per unit width q and average aperture e [11]. In short, the following challenges still exist for the study of fracture seepage: how to quantify the effect of rough fracture wall surface on seepage and the effect of relative roughness on seepage; how to evaluate the applicability of LCL in rough fractures and its characterization equations.
The experimental method can be used to study the above problems, but in practice, due to the temporal and spatial variability of fracture media, the closure of fractures and the randomness of wall surface roughness and groundwater movement in fractures being quite complex, the experimental method is generally not viable. In recent years, with the development of technology and the application of new technologies, some scholars began to apply numerical simulation methods to the research of hydraulic engineering, groundwater and fracture seepage [12][13][14][15][16]. The advantage of numerical simulation method is that it is easy to set parameters of fracture and seepage, such as the fracture aperture, the shapes and densities of roughness elements on fracture wall surface, and a variety of seepage conditions can be simulated without the high cost associated with the experimental method [17,18]. In this paper, Fluent numerical simulation software [10] was used to simulate the following seepage characteristics in rough single fracture under the conditions Other researchers had performed research work before and after the formal establishment of LCL theory, and proposed and developed the relation among discharge per unit width q, fracture representative aperture e (including hydraulic fracture aperture e h , average aperture e, and mechanical aperture e m ), and fracture roughness (including relative roughness and absolute roughness), which modified and supplemented LCL theory (see Section 2.1). Xu et al. found that, when the average fracture aperture e was used to replace the hydraulic fracture e h , there was a super-cubic or sub-cubic relation between the discharge per unit width q and average aperture e [11]. In short, the following challenges still exist for the study of fracture seepage: how to quantify the effect of rough fracture wall surface on seepage and the effect of relative roughness on seepage; how to evaluate the applicability of LCL in rough fractures and its characterization equations.
The experimental method can be used to study the above problems, but in practice, due to the temporal and spatial variability of fracture media, the closure of fractures and the randomness of wall surface roughness and groundwater movement in fractures being quite complex, the experimental method is generally not viable. In recent years, with the development of technology and the application of new technologies, some scholars began to apply numerical simulation methods to the research of hydraulic engineering, groundwater and fracture seepage [12][13][14][15][16]. The advantage of numerical simulation method is that it is easy to set parameters of fracture and seepage, such as the fracture aperture, the shapes and densities of roughness elements on fracture wall surface, and a variety of seepage conditions can be simulated without the high cost associated with the experimental method [17,18]. In this paper, Fluent numerical simulation software [10] was used to simulate the following seepage characteristics in rough single fracture under the conditions of different shapes, distribution densities, and fracture apertures: (1) the influence of fracture wall roughness on fracture seepage; (2) the effect of fracture relative roughness on seepage; (3) super-cubic and sub-cubic relations of seepage in rough fractures; and (4) the scale effect of permeability coefficient K in rough fracture.

Theoretical Background
The current research on fracture seepage is mainly based on the LCL: where: q is discharge per unit width, mm 2 /s; g is gravitational acceleration, m/s 2 ; e is fracture aperture, mm; µ is the kinematic viscosity coefficient of water, mm 2 /s; J is the hydraulic gradient, dimensionless. Hydraulic gradient refers to the ratio of head loss along the seepage path to the length of the seepage path; it can be understood as the mechanical energy lost by the water flow to overcome the friction resistance through the infiltration path of unit length; or the driving force that makes water flow at a certain velocity to overcome friction; or is the head drop per unit distance along the flow direction in the aquifer (the ratio of the water level difference at any two points to the distance between the two points). For the establishment, development, and improvement of LCL, many researchers had introduced absolute roughness ∆ and relative roughness δ, and established the relations among discharge per unit width q, fracture aperture e and fracture roughness, so that LCL could be applied to study the seepage in rough fractures.
Lomize [19]: Flow of turbulent : q = e gJe[2.6 + 5.1 log 10 Louis [20]: Flow of laminar : q = g(e) 3 12µ J 1 Flow of turbulent : q = 4e gJe log 10 [1.9( ∆ 2e ) Amadei and Illangasekare [21]: Su et al. [22]: The above formulas can be generalized as [11,23]: In the above equations: q is discharge per unit width, mm 2 /s; e is average fracture aperture, mm; e is the representative fracture aperture, mm; g is gravitational acceleration, m/s 2 ; µ is kinematic viscosity coefficient of water, mm 2 /s; J is hydraulic gradient, dimensionless; C is the unit conversion parameter to ensure that the dimensions on both sides of the equation are consistent; σ e is mean square deviation of fracture aperture, mm 2 ; ∆ is absolute roughness of fracture, i.e., the height of roughness element a in Figure 1, mm; δ is relative roughness, δ = ∆/e m , i.e., ratio of roughness element height ∆ to fracture aperture e m , dimensionless; n, m, η are indices of fracture aperture e, hydraulic gradient J Appl. Sci. 2022, 12, 7328 4 of 16 and relative roughness ∆, respectively, and their values obtained by different researchers are quite different; ε is a coefficient of δ. Figure 1 illustrates the effect of fracture wall roughness elements on seepage. It can be seen that the streamline (blue streamline) near the rough surface is significantly affected by the roughness elements (rough surface). The blue streamline is curved, and vortexes occur in front and behind the roughness elements. The streamline (yellow streamline) further away from the rough surface will also be affected by the roughness elements (rough surface), but the degree of influence has been weakened. The streamline (green streamline) furthest from the rough surface and near the smooth surface is hardly affected by the roughness elements (rough surface), and the flow line is almost smooth [10].

Numerical Simulation Model
The purpose of this study is to research the effect of rough wall surface of fracture on seepage flow, and to explore the seepage characteristics in symmetrical (but uncoupled) rough fractures [11]. The model with symmetrical-uncoupled rough wall surfaces was devised for this study ( Figure 2). In order to simplify the calculation, the axisymmetric centerline was the symmetrical boundary, and there is no exchange of mass and heat on the symmetrical boundary.
Δ is absolute roughness of fracture, i.e., the height of roughness element a in Figure  δ is relative roughness, δ = Δ/em, i.e., ratio of roughness element height Δ to fractu ture em, dimensionless; n, m, η are indices of fracture aperture e, hydraulic gradie relative roughness Δ, respectively, and their values obtained by different researc quite different; ε is a coefficient of δ. Figure 1 illustrates the effect of fracture wall roughness elements on seepag be seen that the streamline (blue streamline) near the rough surface is significa fected by the roughness elements (rough surface). The blue streamline is curved, a texes occur in front and behind the roughness elements. The streamline (yellow line) further away from the rough surface will also be affected by the roughness e (rough surface), but the degree of influence has been weakened. The streamline streamline) furthest from the rough surface and near the smooth surface is hardly by the roughness elements (rough surface), and the flow line is almost smooth [10 The purpose of this study is to research the effect of rough wall surface of fra seepage flow, and to explore the seepage characteristics in symmetrical (but unc rough fractures [11]. The model with symmetrical-uncoupled rough wall surfaces vised for this study ( Figure 2). In order to simplify the calculation, the axisymme terline was the symmetrical boundary, and there is no exchange of mass and hea symmetrical boundary. Fluent numerical simulation software was applied to simulate the seepage i fracture. The main variables included the distribution density and shape of rou elements, and the fracture aperture. The length of the 2D fracture model (seepag was 100 mm, and the height (fracture aperture) was 3, 4, 5, 6, and 7 mm. Regardle width in the Z direction, it was 1000 mm by default. The dimension a was the h roughness elements and a = 1 mm, b was the spacing between two adjacent ro elements, and the distribution density of roughness elements was defined as A = value of a was a fixed value of 1 mm, and the value of b (roughness elements spaci set to 4 mm, 5 mm, and 6 mm, respectively, so the value of A was 4, 5, and 6, respe Different shapes roughness elements (triangular, rectangular, and sinusoidal) w lected to simulate the influence of rough fracture wall roughness on the change field [10].

Numerical Simulation Model
Grid scale is a problem that should be paid attention to in the research. It is propriate to be too large or too small. When the grid size is too big, it will cause simulation accuracy. The minimum fracture aperture is 3 mm, the roughness height is 1 mm, and the minimum roughness element spacing is 4 mm, so the g should be much less than 1 mm. However, if the grid size is too small, calculati would increase dramatically. Thus, in this study, the grid size selected 0.2 × 0.2 m Fluent numerical simulation software was applied to simulate the seepage in rough fracture. The main variables included the distribution density and shape of roughness elements, and the fracture aperture. The length of the 2D fracture model (seepage path) was 100 mm, and the height (fracture aperture) was 3, 4, 5, 6, and 7 mm. Regardless of the width in the Z direction, it was 1000 mm by default. The dimension a was the height of roughness elements and a = 1 mm, b was the spacing between two adjacent roughness elements, and the distribution density of roughness elements was defined as A = b/a. The value of a was a fixed value of 1 mm, and the value of b (roughness elements spacing) was set to 4 mm, 5 mm, and 6 mm, respectively, so the value of A was 4, 5, and 6, respectively. Different shapes roughness elements (triangular, rectangular, and sinusoidal) were selected to simulate the influence of rough fracture wall roughness on the change of flow field [10].
Grid scale is a problem that should be paid attention to in the research. It is not appropriate to be too large or too small. When the grid size is too big, it will cause loss of simulation accuracy. The minimum fracture aperture is 3 mm, the roughness element height is 1 mm, and the minimum roughness element spacing is 4 mm, so the grid size should be much less than 1 mm. However, if the grid size is too small, calculation time would increase dramatically. Thus, in this study, the grid size selected 0.2 × 0.2 mm. The grid size was much smaller than 1 mm (roughness element height), so the influence of grid size on the flow field could be ignored. The model and grid division structure are shown in

Parameters of Numerical Simulation
The numerical simulation parameters assigned to computer software Fluent are listed in Table 1.

Parameters of Numerical Simulation
The numerical simulation parameters assigned to computer software Fluent are listed in Table 1.

Fracture Roughness, Discharge per Unit Width, and Hydraulic Gradient
In this paper, 2D model was used to simulate the seepage in rough single fracture. The simulation roughness elements' shapes analyzed were triangular, rectangular and sinusoidal. The height of roughness elements, i.e., absolute roughness, was 1 mm, the density of roughness elements A (A = b/a, where a = 1 mm, b = 4, 5, 6 mm) was 4, 5, and 6 respectively, and the fracture aperture was 3, 4, 5, 6, and 7 mm, respectively. The simulation results of rough fracture seepage under three roughness element densities, three roughness element shapes, and five fracture apertures are shown in Figures 6-8.
* All parameters are the same as the simulation parameters in reference [10]. Reprinted/adapted with permission from Ref.
[Simulation on the water flow affected by the shape and density of roughness elements in a single rough fracture]. Copyright year 2019, copyright owners' name Qing Zhang & Jiazhong Qian. ** All italics indicate variables whose parameter values can be set. Non italicized words are explanatory notes.

Fracture Roughness, Discharge per Unit Width, and Hydraulic Gradient
In this paper, 2D model was used to simulate the seepage in rough single fracture. The simulation roughness elements' shapes analyzed were triangular, rectangular and sinusoidal. The height of roughness elements, i.e., absolute roughness, was 1 mm, the density of roughness elements A (A = b/a, where a = 1 mm, b = 4, 5, 6 mm) was 4, 5, and 6 respectively, and the fracture aperture was 3, 4, 5, 6, and 7 mm, respectively. The simulation results of rough fracture seepage under three roughness element densities, three roughness element shapes, and five fracture apertures are shown in Figures 6-8.  Equation (8) shows that there is a power exponential function relation between discharge per unit width q and hydraulic gradient J as follows: where Kc is coefficient of hydraulic gradient; m is termed hydraulic gradient index (Table  S1).
According to the simulation results under the above different conditions, the fitting parameters Kc and m were obtained, as shown in Table 2. It can be seen from the simulation results (Figures 6-8) that the discharge per unit width q is mainly affected by the fracture aperture: a larger fracture aperture results in larger seepage flow. At the same time, the flow is also affected by the roughness element density and relative roughness. Table 2 shows that the variation range of m value is 0.586 Equation (8) shows that there is a power exponential function relation between discharge per unit width q and hydraulic gradient J as follows: where K c is coefficient of hydraulic gradient; m is termed hydraulic gradient index (Table S1). According to the simulation results under the above different conditions, the fitting parameters K c and m were obtained, as shown in Table 2. It can be seen from the simulation results (Figures 6-8) that the discharge per unit width q is mainly affected by the fracture aperture: a larger fracture aperture results in larger seepage flow. At the same time, the flow is also affected by the roughness element density and relative roughness. Table 2 shows that the variation range of m value is 0.586~0.934, all less than 1. When the shape and density of the roughness elements remain unchanged, a decrease in fracture aperture will cause an increase of m value, and it means that, in narrower fractures, the influence of hydraulic gradient is more obvious. When the shape, height of the roughness elements, and fracture aperture remain unchanged, a higher density of the roughness elements will also cause an increase of m value, which means that, in rougher fracture, the influence of hydraulic gradient is also more obvious. The above results indicate that in wider fractures the flow rate mainly depends on fracture aperture, while in narrow and close fracture medium (fracture aperture is very small), the surface roughness of the fracture wall becomes the main factor of head loss of seepage.

Fracture Aperture and Hydraulic Gradient
The average fracture aperture and the corresponding hydraulic gradient index m are plotted as shown in Figure 9.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 10 of 16 ~ 0.934, all less than 1. When the shape and density of the roughness elements remain unchanged, a decrease in fracture aperture will cause an increase of m value, and it means that, in narrower fractures, the influence of hydraulic gradient is more obvious. When the shape, height of the roughness elements, and fracture aperture remain unchanged, a higher density of the roughness elements will also cause an increase of m value, which means that, in rougher fracture, the influence of hydraulic gradient is also more obvious.
The above results indicate that in wider fractures the flow rate mainly depends on fracture aperture, while in narrow and close fracture medium (fracture aperture is very small), the surface roughness of the fracture wall becomes the main factor of head loss of seepage.

Fracture Aperture and Hydraulic Gradient
The average fracture aperture and the corresponding hydraulic gradient index m are plotted as shown in Figure 9. The results show that there is a negative power exponential relation between the hydraulic gradient index m and the average fracture aperture, i.e., as the fracture aperture increases, the influence of hydraulic gradient decreases. The simulation values under different fracture conditions, whether the shape or the density of the roughness elements changes, would have a similar behavior: with an increase fracture aperture, the relative roughness of fracture decreases, and the m value also decreases. A decrease of fracture relative roughness occurs, since the height of roughness elements remains unchanged The results show that there is a negative power exponential relation between the hydraulic gradient index m and the average fracture aperture, i.e., as the fracture aperture increases, the influence of hydraulic gradient decreases. The simulation values under different fracture conditions, whether the shape or the density of the roughness elements changes, would have a similar behavior: with an increase fracture aperture, the relative roughness of fracture decreases, and the m value also decreases. A decrease of fracture relative roughness occurs, since the height of roughness elements remains unchanged while the fracture aperture increases, which means the height (relative roughness) is decreasing relatively.

Super-Cubic Phenomenon of Seepage in Rough Fractures
According to Equation (8), there is a power function relation between the ratio of discharge per unit width to hydraulic gradient q/J m and the average fracture aperture e, which can be expressed as: where K e is the parameter; n is fracture aperture index. The relations between q/J m and average fracture aperture e were drawn according to the simulation results, as shown in Figure 10.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 11 of 16 while the fracture aperture increases, which means the height (relative roughness) is decreasing relatively.

Super-Cubic Phenomenon of Seepage in Rough Fractures
According to Equation (8), there is a power function relation between the ratio of discharge per unit width to hydraulic gradient q/J m and the average fracture aperture e , which can be expressed as: where Ke is the parameter; n is fracture aperture index.
The relations between q/J m and average fracture aperture e were drawn according to the simulation results, as shown in Figure 10.  Table 3.
In the three types of rough fractures, the variation range of coefficient Ke is within the same order of magnitude, and its maximum value is within three times of the minimum value. The curves of different A values (roughness element density) in Figure 10 are very close to one another or even coincide. Therefore, only the fracture aperture index n is discussed here. The ratio of discharge per unit width to hydraulic gradient q/J m and average  Table 3. In the three types of rough fractures, the variation range of coefficient K e is within the same order of magnitude, and its maximum value is within three times of the minimum value. The curves of different A values (roughness element density) in Figure 10 are very close to one another or even coincide. Therefore, only the fracture aperture index n is discussed here. The ratio of discharge per unit width to hydraulic gradient q/J m and average fracture aperture e obtained based on the simulation results were fitted through regression analysis. Coefficient of determination and fracture aperture index n corresponding to different shapes and densities of roughness elements are shown in Table 3.
In the fracture of triangular roughness elements, the fracture aperture index n values are 4.554-4.764; in the fracture of rectangular roughness elements, they are 4.090-4.214, and, in the fracture of sinusoidal roughness elements, they are 3.299-3.608. All are greater than 3, i.e., there is a super-cubic relation in symmetrical-uncoupled fractures between the discharge per unit width q and average aperture e. The simulation results also verify Xu's research conclusion [11].
When studying the fracture aperture index n, the change of average fracture aperture e affects the value of relative roughness (negative correlation, the increase/decrease of the former will cause the decrease/increase of the latter), and the change of roughness density A essentially causes the change of relative roughness (positive correlation, the former and the latter increase and decrease at the same time), so it is impossible to qualitatively describe the effect of the change of roughness density A on the fracture aperture index n by controlling variables. This is also verified by the relation between the ratio q/J m of discharge per unit width to hydraulic gradient and average fracture aperture e in Figure 10: although they increase together on the whole, the increasing trend is difficult to quantify.

Effect of Rough Element Shape and Density on Permeability Coefficient
In the study of groundwater dynamics, especially in the solution of flow equation, it is generally considered that the permeability coefficient K is a constant. However, when evaluating the scale effect of dispersion, it was found that the permeability coefficient K varied with the scale of interest, i.e., K affected by scale effect. The ideal smooth parallel plates fracture can be considered as uniform medium, and the fracture permeability coefficient has no scale effect. However, the surface of the actual fracture is rough and uneven. Due to the existence of wall roughness elements, the roughness and bending degree of the fracture surface are difficult to predict. The flow movement in the real single fracture has strong heterogeneity, so the scale effect of permeability coefficient K is also reasonable.
According to the theory of hydrodynamics, in the fully developed turbulent flow, the relation between the average velocity and the hydraulic gradient is as follows: where V is average velocity, K is permeability coefficient of fracture media, and J is hydraulic gradient. Hydraulic gradient J can be expressed as: where H 1 and H 2 are the corresponding piezometric head values at points 1 and 2, respectively, and L is the horizontal distance between points 1 and 2. Then, the expression of permeability coefficient K can be obtained from Equations (11) and (12): Assuming that the seepage flow is one-dimensional, Equation (13) can be used to obtain the permeability coefficient K. This assumption is reasonable in the ideal smooth parallel plates fracture, but considering the influence of the roughness of the fracture wall on the flow, many small vortices will occur in the flow after bypassing the roughness elements, and they will hinder the seepage movement. Therefore, the two-dimensional method should be used to study the seepage movement in rough fractures. Although Equation (13) may no longer be suitable for studying the scale effect of permeability coefficient K in rough fractures, it can still be used to study the correlation between permeability coefficient K and the horizontal distance L from its corresponding point to the inlet, so as to indirectly study the scale effect of permeability coefficient K. In this paper, the simulation was carried out under the condition that the inlet velocity was set to 0.1 m/s in the fractures with triangular, rectangular, and sinusoidal roughness element shapes and density of 4, 5, and 6, respectively, and the observation points on the symmetrical center line were selected. Under the conditions of different roughness element densities, the relation in each observation point calculated by Equation (13) between K and L was determined. The results are shown in Table 4. As shown in Table 4, the relationship between the permeability coefficient K and the corresponding horizontal distance L can be expressed as a linear function of K = α + βL. Here, α represents the intercept value, and its variation range is −0.521-−0.126; β represents the slope of L, which varies from 5.99-20.05. Parameters α and β are related to fracture aperture, fracture roughness, and hydraulic gradient. In order to illustrate the influence of roughness element density on the scale effect of permeability coefficient K, the simulation results are presented in Figure 11.
As shown in Figure 11, the K-L relation diagrams under the conditions of different densities and shapes of roughness elements are drawn respectively. It can be seen that the farther away from the inlet, the greater the permeability coefficient, which indirectly proves that the permeability coefficient K has a scale effect. The closer to the outlet, the faster the K value increases. The roughness elements with different shapes have similar variation trends. In the same shape of roughness element fracture medium, with the increase of roughness element density, the influence of roughness element density on the permeability coefficient K becomes greater. It can also be seen from Table 4 that, in the data corresponding to the triangle roughness element fracture, when the density value A is 6, the correlation coefficient is 0.417, when A is 5, the coefficient increases to 0.577, and when A is 4, the coefficient decreases to 0.404; in the rectangular roughness element, when A is 6, the correlation coefficient is 0.542, and when A is 5 and 4, the corresponding correlation coefficients are 0.562 and 0.458, respectively. In the sinusoidal roughness element fracture, the correlation coefficient corresponding to A of 6 is 0.623, and the correlation coefficients corresponding to A of 5 and 4 are 0.790 and 0.708, respectively. is 6, the correlation coefficient is 0.417, when A is 5, the coefficient increases to 0.577, and when A is 4, the coefficient decreases to 0.404; in the rectangular roughness element, when A is 6, the correlation coefficient is 0.542, and when A is 5 and 4, the corresponding correlation coefficients are 0.562 and 0.458, respectively. In the sinusoidal roughness element fracture, the correlation coefficient corresponding to A of 6 is 0.623, and the correlation coefficients corresponding to A of 5 and 4 are 0.790 and 0.708, respectively. With the increase of roughness elements distribution density A, the correlation coefficient between permeability coefficient K and corresponding horizontal distance L first increases and then decreases, which may be due to the coupling influence of fracture aperture, fracture roughness element height, and roughness element distribution density on the scale effect of permeability coefficient K. Under the condition that the fracture aperture and roughness element height remain unchanged, the influence of roughness element distribution density A on the correlation between K and L also increases first and then decreases. This is because the fracture wall tends to be smooth when the distribution of rough elements is very dense or sparse, while the influence of the distribution of roughness elements on the flow in the real case is between the two.

Conclusions
Based on the Fluent numerical simulation software and a 2D numerical model, this paper studies the rough fracture seepage characteristics under the conditions of different With the increase of roughness elements distribution density A, the correlation coefficient between permeability coefficient K and corresponding horizontal distance L first increases and then decreases, which may be due to the coupling influence of fracture aperture, fracture roughness element height, and roughness element distribution density on the scale effect of permeability coefficient K. Under the condition that the fracture aperture and roughness element height remain unchanged, the influence of roughness element distribution density A on the correlation between K and L also increases first and then decreases. This is because the fracture wall tends to be smooth when the distribution of rough elements is very dense or sparse, while the influence of the distribution of roughness elements on the flow in the real case is between the two.

Conclusions
Based on the Fluent numerical simulation software and a 2D numerical model, this paper studies the rough fracture seepage characteristics under the conditions of different shapes, densities, and fracture apertures of roughness elements. The main conclusions are as follows: 1.
In wider rough fractures, the flow rate mainly depends on fracture aperture, while in narrow and close rough fracture medium, the surface roughness of fracture wall becomes the main factor of head loss of seepage.

2.
In rough fracture, there is a negative power exponential relation between the hydraulic gradient index m and the average fracture aperture e, i.e., with the increase of fracture aperture e, the relative roughness of fracture and the influence of hydraulic gradient both decrease. 3.
In symmetrical-uncoupled rough fractures, there is a super-cubic relation between the discharge per unit width q and average aperture e.

4.
The value of rough fracture permeability coefficient K is not a constant, and it is affected by the scale effect (the horizontal distance from the measuring point to the fracture inlet) and the density of the roughness elements.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/app12147328/s1, Table S1: Relation between discharge per unit width q and hydraulic gradient J.