Numerical Simulation on Non-Darcy Flow in a Single Rock Fracture Domain Inverted by Digital Images

,


Introduction
Many important underground resources, such as groundwater, oil, gas, coalbed methane, and geothermal energy, are exploited in low-permeability reservoirs with highly developed fractures [1,2]. Hydraulic fracturing of oil-and gasrich strata and coal seams requires accurate control of the amount of fluid injected into the rock fractures. The deep burial of nuclear waste and CO 2 geological storage should reduce the fluid infiltration into the rock mass as much as possible [3]. In situ leaching (ISL), an alternative mining technology, requires the injection of leaching solution into the artificial fracture to dissolve target minerals in impervious host rocks [4]. Compared with intact low-permeability rocks, the fractures formed by rock failure greatly improve its permeability. The fractures formed by rock cracked change the stability and permeability of surrounding rocks, which are easy to induce geological disasters [5][6][7][8][9][10]. Obtain-ing the fracture flow characteristics is the premise of rock seepage control. Therefore, the fluid flow through fractures of the rock mass has always been the focus of engineering research [11].
The classical cubic law of linear laminar flow was developed in the early study of rock fracture flow using a smooth parallel plate model. It was used to evaluate the flow capacity of fractures [12,13]. However, due to the complexity of the fracture geometry and flow regime, the cubic law does not adequately describe fluid flow behavior in natural fractures, and non-Darcy flow may occur as a result of nonnegligible inertial losses. Previous experimental investigations reported that Darcy's law fails to predict pressure drops in fractures when inertial effects are relevant before the fully developed turbulence [14,15]. The rough structure of the fracture surface causes non-Darcy flow [16], and the exact solution must be obtained by solving the Navier-Stokes equation, which is difficult to obtain in engineering applications [17]. Therefore, scholars have tried to develop a characterization method for determining the fracture roughness and the reduction of the coarse structure of natural fractures to conduct non-Darcy flow behavior research. Su et al. [18] used artificial rough surface fractures to simulate natural rock mass fractures, but the actual rough fracture surface is much more complicated than the artificial rough surface. Barton and Choubey [19] used the joint roughness coefficient (JRC) to define the roughness and created 10 standard curves, which quantitatively reflect the fracture roughness. Ju et al. [20] used polymethyl methacrylate (PMMA) to make planar models of fractures with different roughness and used a high-speed camera to record the entire rough fracture water seepage process. Xie et al. [21] used a laser scanner system with a spacing grid of 0.1 mm to conduct a two-dimensional numerical simulation of single fractures during shear displacement with the Navier-Stokes equations. However, a fracture is a three-dimensional space with irregular scale and characteristics, so only the accurate inversion of a three-dimensional fracture domain with actual rough surfaces can objectively reflect the geometric characteristics of a natural fracture.
As an accurate spatial structure measurement and digital representation on the microscopic scale of materials, digital image processing technology has been widely used in the analysis of geotechnical microstructures. Angelin et al. [22] used image processing techniques of K-means algorithms and Watershed algorithms to analyze microscopic images for void identification in cement matrices. Based on digital image analysis, Thomas et al. [23] analyzed the distribution and size of macroporosity under different mixing conditions by computerized axial tomography, scanning electron microscopy, and a new developed methodology. Zhu et al. [24,25] utilized the digital image technique to numerically characterize the heterogeneity of the structural characteristics of coal and rock masses, and they studied the influence of heterogeneity on coal rock fractures and seepage using numerical analysis. In view of the fact that the roughness of the fracture surface results in different pixel values, in this paper, we attempt to restore the rough structure of the fracture surface by converting the pixels into a 0 to 1 data distribution (normalization).
There are many types of non-Darcy flow such as lowvelocity non-Darcy flow caused by boundary layer [26] and high-velocity non-Darcy flow caused by inertia force, of which the latter type is adopted in this manuscript. The Forchheimer equation [27][28][29][30] is commonly used to describe non-Darcy flow. Because the coefficients in the Forchheimer equation are closely related to the geometric characteristics of fractures [26,31], the variation in the coefficients is a necessary condition for describing non-Darcy flow. Many studies have found that surface roughness determines the nonlinear characteristics of fluid flow in natural rock fractures [32,33]. Xia et al. [34] observed that non-Darcy seepage in rough fractures depends on the void space and composite morphology of the fracture surface. Yin et al. [35][36][37] reported a large number of experimental studies on nonlinear flow characteristics of fractured rock samples. Zhang et al. [38] presented a hydromechanical framework for porous materials exhibiting two dominant porosity scales that accommodates transverse isotropy induced by distributed microfractures and non-Darcy flow through the nanometer-scale pore spaces. Therefore, it is necessary to quantitatively study the influence of the fracture's geometric characteristics on non-Darcy flow behavior.
In this work, the reverse model for rock fracture domain will be carried out according to the digital images of rock fracture surfaces, and the model will be implemented to numerical simulation to realize the visualization of non-Darcy flow in rock fracture. The simulation results are fitted according to the Forchheimer equation describing the nonlinear flow, and the effects of the geometric characteristics of the fracture domain on the non-Darcy flow are further analyzed. This paper provides a reference for the current methods of non-Darcy flow in rough rock fractures.

Inversion of the Fracture Domain Based on Digital Images
The digital image is the reflection of the objective object, which is composed of many pixels with a matrix arrangement, so the digital image can reflect the fluctuation height of the rock fracture surface. In gray level images or binarization images [39,40], the gray values are 0-1 and 0-255, respectively. The discrete function of the gray value or chromaticity of the corresponding digital image can accurately reflect the distribution characteristics of the material surface, which is the basis of inversion of the fracture domain based on digital image technology. The process of inverting the fracture domain based on the digital image is shown in Figure 1. Fractured rocks form fracture domains of different scales, which are actually surrounded by rough rock surfaces. As shown in Figure 2, a single fracture domain is formed between the upper and lower surfaces. Therefore, we used CCD camera to get the images of rock fracture surfaces, which are 24-bit true color pictures with pixel size 1000 × 600. In addition, the images are denoised to reduce the interference of external factors.

Geofluids
The images of the fracture surface have a wealth of information, which can well reflect the mesostructure of the surface through different grayscales and colors. Digital image is composed of rectangular image elements, also known as pixel points. The digital image of the surface of the rock sample can be stored as an image of different accuracy by image processing software. The so-called accuracy here is controlled by pixels. There are many kinds of color space such as RGB composed of three primary colors and HSI composed of three variables, i.e., hue, saturation, and illuminance [41,42]. To facilitate the transition of both kinds of color space, illuminance (I) is generally defined as the arithmetic mean of the three RGB components, i.e., I = ðR + G + BÞ/3 [43,44], which is also used in this paper. In the digital image of the rock fracture surface, the values of R, G, and B range from 0 to 255. After the image is normalized, the values of R, G, and B are converted to the I value, which has a 0-1 distribution. It shows the normalized result of the rock fracture surface images in Figure 3. The lowest pit on the rock fracture surface is considered to be in the horizontal plane, so the I value is always distributed between 0 and 1.
In this section, the fracture reconstruction method of Zhao et al. [45] is referenced and improved for inversion of the fracture domain. For the lower surface of the fracture, the digital image matrix function of the fracture surface is obtained and recorded as ½B, and then, ½B is normalized to obtain the normalization matrix ½ B. The minimum and maximum values in ½B are determined using MATLAB and are recorded as b min and b max , respectively. Thus, the height distribution of the lower surface of the fracture is as follows: For the upper surface of the fracture, the digital image matrix function of the fracture surface is obtained and recorded as ½T, and then, normalization processing is carried out to obtain the normalization matrix ½ T. Using the same method, the minimum t min and maximum t max in ½T are determined. Thus, the surface height distribution on the fracture is as follows: As shown in Figure 4, the height distribution of the upper and lower surfaces of the fracture is restored. By comparison, the results of the restored surfaces are consistent with the actual height distribution on the fracture surfaces. When further inversion of the fracture domain is needed, the relative positions of the upper surface and the lower surface need to be determined. As shown in Figure 5, we select the horizontal plane, where b min is located, as the reference plane I with height h 1 , so the height function of any point on the lower surface of the fracture will be Then, we select the horizontal plane, where t max is located, as the reference plane II with height h 2 ,s ot h e height function of any point on the surface of the same fracture will be By combining the top and bottom of the fracture specimen (see Figure 6), the rock fracture domain is obtained, and the spatial distribution function of the fracture aperture eðxÞ will be:

Model Implementation
3.1. Mathematical Model. The single-phase flow of incompressible fluid in rock fractures can be expressed by mass conservation and Navier-Stokes equations, which can be expressed as [12] ρ ∂u ∂t where u, ρ, ∇P, μ, and F represent the flow velocity vector (m/s), the fluid density (kg/m 3 ), the fluid pressure gradient along the flow direction (MPa/m), the viscosity coefficient (N·s/m 2 ), and the body force vector (N), respectively. Due to the increasingly obvious inertia effect as the flow rate and fracture roughness increase [28,46], the fluid usually exhibits non-Darcy flow behavior in rough rock fractures. In general, Forchheimer equation in porous media seepage theory is introduced to describe the non-Darcy flow [26], which is defined as where Q (m 3 /s) is the flow rate; A (kg·s -1 ·m -5 ) is the coefficient of the linear term; B (kg·m -8 ) is the coefficient of the nonlinear term. Both coefficients A and B rely on knowledge of the fluid properties and the geometric characteristics of rough fractures [29], which can be expressed as where β (m -1 ) is the non-Darcy coefficient, which varies with the geometric characteristics of the fractures. Reynolds number Re is a hydraulic parameter, which is a dimensionless number used to judge the flow state of viscous fluid. Its physical meaning is the ratio between the inertial force and the viscous force of the fluid, which can be expressed as In order to further explain the mechanism of non-Darcy flow, Zeng and Grigg [47] proposed a non-Darcy flow effect factor E, which is defined as Fracture aperture e(x) Height distribution

Geofluids
The non-Darcy flow effect factor E indicates the degree of non-Darcy flow, which is between 0 and 1. Combining equations (8), (9), (10), and (11), the Reynolds number Re can be rewritten as The critical Reynolds number Re c characterizes the onset of the flow transition from linear flow to nonlinear flow. In recent studies, the critical condition for the onset of nonlinear flow has been defined as the point at which the nonlinear pressure drop contributes 10% to the overall pressure drop [28], which is equal to E =0:1. Considering this condition, the critical Reynolds number for nonlinear fluid flow in rough fractures was suggested by Javadi et al. [29]: 3.2. Numerical Model. We consider a 3D single fracture domain 50 mm wide and 90 mm long, and the roughness varies with different fracture domains. The geometric profile and boundary conditions are shown in Figure 7. The fracture aperture distribution of each single fracture domain can be obtained by equation (5). In order to analyze the influence of fracture aperture on non-Darcy flow, the average value of fracture aperture distribution is selected as the relative description of fracture aperture size, in which of the model is 10 mm. We assume a single fluid flow in the fracture domain, the pressure at the inlet boundary on the left was 1.0 MPa, and the pressure at the outlet boundary on the right was 0.1 MPa. The other boundaries are set as no flow boundaries. The following water parameters are assumed in the simulations: temperature T =20°C; density ρ = 998:2 kg/m 3 ; dynamic viscosity μ =0:001 Pa · s.
In this work, the COMSOL Multiphysics code is used for the numerical model implementation based on the Finite Element Method (FEM). The model was divided into 133082 grids using free tetrahedral node. The numerical simulations are obtained by computational convergence of the stationary studies.
Using the rock fracture domain inversion method, we obtained 4 fracture domains with different roughness as shown in Figure 8. Table 1 lists the JRC range of some fracture contours of these 4 samples, which can be used to evaluate the overall roughness of the fracture domains.

Results and Discussion
4.1. Effects of the Different JRCs. Figure 9 shows the velocity distribution in the fracture domains with different JRCs. Overall, the undulating structure of the fracture surface makes the velocity distribution very uneven. On the upper surface of the fracture, each protruding position is relatively low, while the concave position has a relatively high velocity. There is a low-speed boundary layer at the entrance of each sample, which surrounds the high-speed mainstream area. As the JRC increases, the effect of the boundary layer becomes more significant, resulting in an uneven velocity distribution at the entrance. In Figure 9, because the fracture surface of Sample 1 is relatively smooth, the velocity fluctuation is not significant, and the maximum velocity is up to 0.6 m/s. As the surface roughness increases, the fluctuation in the fracture surface increases gradually, and the maximum velocity decreases gradually. The maximum velocity of Sample 4 is only 0.4 m/s. The essence of non-Darcy flow is that the growth of the flow and the pressure gradient no longer satisfies a linear relationship. In order to investigate this nonlinear flow behavior, different water pressures were set at the inlet boundary in the simulations. The relationship between the hydraulic gradient and the flow rate is shown in Figure 10. −∇P represents the macroscopic pressure gradient along the flow direction, which is equal to the pressure drop between the inlet and outlet divided by the fracture length l. Based on the relationship between the hydraulic gradient and the flow rate, the relationship between each pressure gradient and flow rate deviates from the linear relationship. Because the influence of the inertial force becomes more significant as the flow rate increases, the degree of deviation increases. When the flow rate is the same, as the JRC increases, the hydraulic gradient increases, the slope between the pressure gradient and the flow rate becomes steeper, and the deviation from the linear relationship increases. When the 5 Geofluids pressure gradient is the same, the flow rates of the four samples are significantly different. The sample with the lowest roughness has the highest flow rate, and the roughness of the other samples increases gradually and the flow rate decreases.
As the fracture roughness increases, the circuitous degree of the fracture flow path increases, and the inertial force of the fluid flow increases. The nonlinear pressure gradient loss caused by the inertial force accounts for most of the total    Figure 11. Overall, there are always relatively low velocity boundary layers within the fracture domains, which surround the high-speed mainstream area. Observing the internal section, there are significant differences in the velocity at different positions in the fracture aperture. For the position where the fracture aperture is small, the velocity increases rapidly, while for the position where the fracture aperture is large, the velocity decreases accordingly. For the same pressure gradient, the maximum velocity is 0.3 m/s when the fracture aperture is 6 mm. The maximum velocity increases as the fracture aperture increases, and the velocity reaches 0.5 m/s when the fracture aperture is 12 mm. The streamline color represents the velocity distribution, and the change in the streamline color reflects the increase of the overall velocity in the fracture domain with increasing fracture aperture. Based on the flow direction, as the fracture aperture increases, the streamline tends to become smoother. As shown in Figure 12, based on the Forchheimer fitting curves of the pressure gradient and the flow rate, the relationship between the pressure gradient and the flow rate is still nonlinear regardless of the change in the fracture aperture. The determination coefficients R 2 of fitting curves are 0.9993, 0.9997, 0.9998, and 0.9996, respectively, indicating that the fitting effect is satisfactory. As can be seen from the diagram, when the flow rate is the same, i.e., the x value of each regression equation is the same, the y value, i.e., the required pressure gradient, decreases with increasing fracture  7 Geofluids aperture e. When the pressure gradient is the same, the flow rate increases with increasing fracture aperture e.
The flow path of a single fracture with the same JRC becomes relatively relaxed and smooth with increasing fracture aperture, which means that the corresponding tortuous degree decreases and the inertial force becomes weaker in the fluid flow. The nonlinear pressure gradient loss caused by the inertial force accounts for part of the total pressure gradient loss and increases the proportion of the fluid kinetic energy. In addition, due to the difference in the fracture aper-ture e, the primary term coefficient A and the quadratic term coefficient B of each regression equation are also different. The variation in these coefficients will be further analyzed in the next section.

Variation of the Coefficients of the Forchheimer Equation.
According to the physical meaning of the Forchheimer equation (equation (7)), A (kg · s −1 ·m −5 ) is the coefficient of the linear term, which represents the energy losses due to viscous dissipation mechanisms; and B (kg · m −8 ) is the coefficient of  8 Geofluids the nonlinear term, which describes the energy losses arising from the inertial effects [26]. Both coefficients A and B rely on knowledge of the fluid properties and the geometric characteristics of rough fractures [14,15]. In order to further analyze the change in the coefficients of the Forchheimer equation, we carried out numerical simulation for all four samples when the fracture apertures were adjusted to 6 mm, 8 mm, 10 mm, and 12 mm, respectively. The relationship between the pressure gradient and the flow rate was extracted, and the results of coefficients A and B were obtained by polynomial fitting (Table 2). Figure 13 shows the relationship between coefficient A of the Forchheimer equation (equation (8)) and the fracture aperture e. For fractures with the same JRC, coefficient A decreases with increasing the fracture aperture e, indicating that the viscous effect has been gradually weakened. For the same fracture aperture, coefficient A also increases with increasing JRC, which also means that the viscous force is increasing. This is consistent with the observation of Xiong et al. [48]. Figure 14 shows the relationship between coefficient B of the Forchheimer equation (equation (7)) and the fracture aperture e. From the point of view of the changing trend, the JRC and the fracture aperture e are closely related to coefficient B. For the same JRC, coefficient B decreases with increasing fracture aperture e, indicating that the fluid inertial force weakens. For the same fracture aperture e, as the JRC increases, coefficient B also increases, which indicates that the inertial force of the fluid is increasing. This observation is consistent with the observation of Chen et al. [49].
Variation of the coefficients A and B shows a much similar pattern. For the same JRC, both coefficients A and B decrease with increasing fracture aperture e. For the same fracture aperture e, both coefficients A and B increase with increasing JRC. Both coefficients A and B of rough samples experience a decrease in 2 orders of magnitude while those coefficients of smooth samples decrease in 1 order of magnitude as the fracture aperture increases from 6 to 14 mm. Table 2 produced from the calculations using equation (13), the distribution of the critical Reynolds number   9 Geofluids Re c is obtained. As shown in Figure 15, the critical Reynolds number Re c increases with increasing fracture aperture e, and rougher samples have smaller critical Reynolds numbers Re c . Combined with the previous analysis, rougher fractures and smaller fracture apertures will make the flow paths more tortuous. Then, the proportion of the nonlinear pressure gradient loss increases. This easily leads to nonlinear flow, resulting in a smaller critical Reynolds number Re c .

Variation of the Critical Reynolds Number. Based on the data in
The critical Reynolds number (CRN) model in equation (14) provides a simple method of clear physical significance to quantify Re c for fluid flow through rock fractures. This model is useful for numerical simulation of fluid flow in fractured networks, in which a decision can be flexibly made to include the nonlinear effect [50].
where λ and m are regression coefficients. As illustrated in Figure 15, simulation data are fitted well the results of the CRN model, which manifest that the numerical simulation can suit for fluid flow in the fracture domain. In further analysis, the effects of fracture aperture and roughness on the coefficients A and B should be related to the physical meaning of the fluid flow process. In fact, the small fracture aperture and the great rough extent are generally accompanied with more tortuous flow paths, which results in significant inertial effects. This statement is similar to many previous studies. For example, Javadi et al. [29] investigated the role of shear processes on nonlinear flow through rough-walled rock fractures, showing that the coefficients A and B experience 4 and 7 orders of magnitude reduction with increasing shear displacement, respectively, mainly as a result of shear dilation (or equivalently, the increase in fracture aperture) of the fractures. Xiong et al. [51] developed a numerical procedure about nonlinear flow in threedimension discrete fracture networks (DFN), showing that both the linear coefficient A and the nonlinear coefficient B of the Forchheimer law decrease with increasing percolation density, but increase with increasing JRC. The experimental results conducted by Ni et al. [15] on a seepage apparatus have shown that the monomial coefficient and the quadratic coefficient decrease with the increase of the fracture aperture, and with the increase of joint roughness coefficient, the non-Darcy influence coefficient of rough fracture increases.

Conclusions
This paper presents a model to investigate non-Darcy flow in single rock fracture domain inverted by digital images, which is further used in flow numerical simulation. In addition, we further discuss the influence of the geometric characteristics of rock fracture domain on the non-Darcy flow. The main conclusions are as follows.
(1) The rough structure of the fracture surface produces different pixel values in the image. The rough undulating structure of the rock fracture surface can be accurately reduced using digital image processing technology, and then, it can be combined with the actual measurement data to determine the relative position of the upper and lower surfaces of the fracture. This method can be used for inversion of rock fracture domains (2) The JRC and fracture aperture significantly influence the fracture fluid flow. For the same conditions, as the JRC increases, the tortuous degree of the fracture flow path increases, the pressure gradient loss caused by the inertial force increases, and the proportion of the kinetic energy of the transformed fluid decreases.
For the same conditions, as the fracture aperture increases, the fluid flow path becomes relaxed, the   10 Geofluids pressure gradient loss caused by the inertial force decreases, and the fluid velocity increases (3) The geometric characteristics of fracture domain have obvious influence on the coefficients of Forchheimer equation. Both coefficients A and B decrease with increasing fracture aperture and increase with increasing JRC. Meanwhile, as the JRC increases or the fracture aperture decreases, the tortuous degree of the fracture seepage path increases, leading to an increase in the proportion of nonlinear pressure gradient loss caused by the inertial force. The critical Reynolds number Re c decreases accordingly, indicating that the nonlinear flow is more likely to occur at this time (4) The conclusion of non-Darcy flow in rock fracture in this paper is consistent with the previous studies, which verifies the feasibility of understanding non-Darcy flow in rock fracture domains by inversion based on digital images Although the inversion model for the fracture domain presented in this paper provides some insights into the investigation of non-Darcy flow behaviors, compared with CT image method, the inversion method of the fracture domain has limitations in accuracy. Thereby, the inversion method of the fracture domain by digital image needs further improvement.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

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