Numerical Modeling of Fluid Flowing Properties through Porous Media with Single Rough Fractures

Due to the existence of pores in the fracture surface, the permeability of pore-fracture media is more complicated, and its permeability law and mechanism are worthy of an in-depth study. The rough fracture surface of this paper is obtained by 3D laser scanning technology, taken a small part of the middle area to carry out the finite element simulation of single fracture seepage, and studied the distribution of flow velocity fields and the nonlinear characteristics by changing the inlet flow rate, fracture thickness, and size of fracture surface. Given that for the same fracture opening, the flow velocity decreases from the middle position to both sides. When the roughness increases locally, the flow velocity suddenly increases. The nonlinear coefficient decreases as the fracture thickness increases. A dual-media model which consists of porous matrix and a single rough surface was established to study the seepage characteristics through two media. Given that it takes less time to reach a steady state when the permeability coefficient of the matrix is greater than that of the fracture. After reaching the steady state, the flow direction of the fluid in the matrix is consistent with the fluid direction in the fracture, and both are perpendicular to the isosurface of pore water pressure.

Focusing on flow behaviors in single rough fractures, researchers have devoted substantial efforts to permeability in rock fractures by separating intact rock blocks with negligible matrix permeability [13][14][15][16][17]. The results indicate that the flow behaviors, flow capacity, and permeability evolution through single fractures are influenced by many factors such as fracture opening, roughness coefficient, and inlet water pressure [6,[18][19][20]. Given the difficulties in monitoring the evolutions of aperture thickness and contact spots between fracture surfaces in the fluid flowing process during the experimental tests, the numerical simulation has been widely adopted recently [21][22][23]. Fluid flowing behaviors through a series of three-dimensional models of rough fractures were simulated by solving the Navier-Stokes equations [24]. Xiong et al. [25] reported integrated laboratory experiments and 3D numerical modeling results of fluid flow through roughwalled rock fractures, which considered the impacts of shear displacement on the transmissivity of fractures. Xiong et al. [26] investigated the influences of normal loading and shearing on hydraulic properties in rough-walled rock fractures. The results demonstrate that the deformation of fracture with increasing normal stress and shear causes nonuniform changes in void space geometry and further influences fracture permeability. These studies generally illustrated the impacts of surface roughness and weak inertial effects on the flow behaviors through 3D rough rock fractures. It should be noted that most of existing publications about simulation studies rarely present the details of flow velocity fields because the modified models considered only simplified fracture surface topography and parameters of specific flow conditions.
Additionally, simplified 2D rough fractures formed by cutting through cross sections of real fractures have been extensively applied to model fracture flows ( [27][28][29][30][31][32][33][34]). However, these simplified 2D models are locally selected, which have significant limitations in considering mechanical effects and nonlinear flow behaviors of the whole 2D rough fracture surfaces. Xiong et al. [25] studied the nonlinear fluid flow behaviors in rough-walled rock fractures by simulating two 3D fracture models, with and without shear, which showed that the effective transmissivity is a function of local apertures with serious uncertainties even when Re (Reynolds numbers) is small. However, the fracture model established in this study has only one fracture opening and the conclusions are not universal. More importantly, in the laboratory tests and numerical simula-tions of abovementioned studies, the effects of matrix porosity are not considered, thus to produce errors for evaluating the overall flow capacity of fractured rocks. Until recently, even though the hydraulic properties of pore-fracture dual media have been studied by scholars [35,36], the nonlinear flow behaviors, evolution of permeability, and contribution of fracture and matrix have rarely been reported.
The purpose of this study is to numerically investigate the flow regimes, nonlinear flow behaviors, and permeability evolution processes of pore-fracture combination models. First, single rough fracture models with different fracture openings (0.1-0.5 mm) were established to investigate the microscopic fluid flow patterns. Then, porefracture combination models with different boundary conditions were established. The pore water pressure in the two models, Darcy's velocity in matrix and fractures, and contributions of matrix and fractures to the overall discharge capacity of the pore-fracture models were individually discussed.

Model Establishment
2.1. Geometrical Models. In this study, a cylindrical granite sample with the size of φ50 × 100 mm was first conducted directly Brazilian splitting to produce a tensile rough-walled fracture [37], as shown in Figures 1(a) and 1(b). The two clamps at the left/right side of the samples in Figure 1(a) can prevent the dislocation in the splitting process. After splitting, the rough fracture surface with a projection size of a length of 100 mm and a width of 50 mm was produced (Figure 1(b)). The fracture surface morphology of the fracture samples after high temperature exposure before the hydromechanical tests was measured using a high-resolution noncontact 3D laser scanning [37,38]. The scanned data, which were automatically recorded in a 0.1 mm internal in the fracture plane, and have a vertical (z-direction) precision of ±1 μm, were then used to digitalize the fracture surfaces, as shown in Figure 1(c). Limited by the challenging computational capacity of solving the NSE for the 3D fracture models with highly fine representation of surface roughness, only a part of the surface at the zone of x = ½25, 35 mm and y = ½35, 55 mm (with the size of 20 mm in length and 10 mm in width) was extracted to form the 3D models for a generic study (Figure 1(d)). This part was chosen in order to avoid the effect of sample edges as well as the extreme asperity heights, thus to avoid only single asperity contact spot appearing in the computational model.
Two three-dimensional models, with/without the matrix, were, respectively, generated for simulations. Therefore, the fluid flow characteristics through a single rough fracture model and the pore-fracture combination model can be clearly demonstrated by the visual models. The single rough fracture model was built by a numerical "lift" step that simply extruded the scanned surface upwards in the vertical (height) direction (z-direction) by 0.1 mm, 0.2 mm, 0.3 mm, 0.4 mm,   3 Geofluids and 0.5 mm, respectively, to produce various mechanical aperture thickness d of the fracture surface ( Figure 1(d)). The pore-fracture combination model was created numerically by Boolean operation between the single fracture models (taking d = 0:2 mm as the research object) with a length of 100 mm and a width of 50 mm using a cylindrical matrix (Figure 1(e)). This numerical step finally created a pore-fracture combination model with a 3D rough-walled fracture surface and matrix, which can be used to simulate thermally treated cylindrical rock samples with a single rough fracture surface in laboratory tests [37].

Fracture Characterization.
The fluid flow behaviors through fractures depend on the fracture geometry characteristics, commonly characterized by the midsurface and aperture geometries. The midsurface is often applied to roughly characterize the tortuous flow paths in the vertical direction (z-axis). The aperture space directly determines the transmissivity and flow-wetted areas [39].
To describe the morphology characteristics of the fracture surface, 3D spatial distribution of the roughness was analyzed [40,41]. Based on the statistical results of the rough rock fracture surfaces using the high-resolution noncontact 3D laser scanning, the 3D spatial distribution parameters, including the asperity height, slope angle, and aspect direction of the fracture surface mesh elements, were, respectively, calculated. The frequency distributions of the asperity height, slope angle, and aspect direction were then evaluated, as shown in Figures 2(a)-2(c). The asperity height frequency distribution histograms present a typical normal distribution, with the standard deviations of 0.9972 mm and the maximum of 2.8242 mm. The slope angle frequency distribution histograms exhibit a log-normal distribution, with the standard deviations of 12.3737°and the maximum of 81.0239°.
The aspect direction frequency distribution shows a polar plot, with the standard deviations of 105.9922°.

Modeling of the Fluid Flow Behaviors
3.1. Governing Equations and Numerical Methods. The governing equations of fluid flowing through single fractures are the well-known Navier-Stokes equations (NSE), which represents mass and momentum conservations. For the isothermal, steady-state, and incompressible single Newtonian fluid flow, the NSE can be written as follows [8]: where ρ (kg/m 3 ), u (m/s), P (Pa), μ (Pa·s), and t (s) denote the density of a fluid, the velocity vector, the pressure, the viscosity coefficient, and time, respectively. The NSE are a set of nonlinear partial differential equations coupled with velocity and pressure fields. In this study, the commercial finite element software of COMSOL Multiphysics was employed to solve them.

Initial and Boundary
Conditions. The basic assumptions used in the numerical simulation are as follows. Firstly, the rock matrix is a typical kind of homogeneous, isotropic linear elastic medium and its deformation belongs to the category of small deformation. Secondly, no crack propagation and dislocation occur in the rock mass in the process of fluid flow tests. Finally, the compressibility of the fluid and the thermal effect of the fluid flow process are ignored. The density and dynamic viscosity coefficient of the fluid are kept constant. The basic assumptions were applied to both single fractures and pore-fracture combination models.
The numerical simulation of fluid flowing through single rough-walled fracture models adopts five parameters, as listed in Table 1, in which d is the fracture aperture thickness, and here in this study, the d value of 0.1, 0.2, 0.3, 0.4, and 0.5 mm was, respectively, selected. Q is the volume velocity at the inlet, L is the length of the fracture surface along the x -axis, and w is the fracture width along the y-axis. Taking d = 0:5 mm as an example, the two boundaries at x = 0 mm and x = 20 mm were set as the water inlet and outlet boundaries, respectively ( Figure 3). The inlet boundary conditions were specified with different values of constant flow rates, and the outlet boundary condition was set as zero pressure, i.e., P = 0. The rest of top, bottom, and side boundaries in the models were all set as no fluid flow and nonslip walls, i.e., u = 0 m/s. Such boundary conditions are consistent with many experimental conditions [6,13,42]. To reduce the inlet boundary effects of constant flow rate (which often leads to a constant flow velocity), a laminar entrance inflow condition, provided by the COMSOL Multiphysics software, was adopted. Using this laminar entrance inflow condition can largely improve the convergence rates, especially in nonlinear flow cases when the Reynolds number Re is relatively high. With such settings, the x-axis direction is the principal direction of fluid flow.
The pore-fracture combination model generated 209805 tetrahedral elements using a free tetrahedral mesh function. In order to compare and analyze the fluid flow characteristics of the pore-fracture model under the combination of different permeability coefficients, the model has two working conditions (the permeability coefficient of the matrix is greater or less than that of the fracture). The geometric parameters adopted by the model of the two working conditions are consistent. The Darcy flow physical field was adopted for both working conditions, and the fracture flow module in the physical field of Darcy flow was adopted for the boundary conditions of the fluid flow of the fracture.
M1 corresponds to the case that the permeability coefficient of the matrix is smaller than that of the fracture. The pore structures in the matrix are relatively developed under high temperature exposure, but the permeability coefficient is still relatively small, which corresponds to the working conditions that the fractured granite sample is under a confining pressure of 10 MPa and a high temperature of 400 degrees [37]. The permeability coefficient K m1 of the rock matrix is 2:04 × 10 −17 m 2 , and the calculated permeability coefficient K f1 of the fracture is 1:695 × 10 −16 m 2 .
M2 corresponds to the case that the permeability coefficient of the matrix is greater than that of the fracture. The pore structures in the matrix are highly developed under high temperature exposure, and the permeability coefficient is relatively large, which corresponds to the working conditions of the fractured granite sample under a confining pressure of 10 MPa and a high temperature of 800 degrees [37]. The permeability coefficient of the rock matrix K m2 is 2:01 × 10 −16 m 2 and the calculated permeability coefficient K f2 of the fracture is 1:204 × 10 −16 m 2 .

Numerical Simulation of Fluid Flow in Single Rough Fractures
The simulated Darcy velocity field is shown in Figure 4 6 Geofluids When the fracture aperture thickness is small, the flow velocity after the steady state is relatively large. The average volume flow velocity through the entire fracture surface for d = 0:1 mm is 7:6778 × 10 −12 m 3 /s, larger than that of 7:4686 × 10 −12 m 3 /s for d = 0:5 mm. This is because the larger the fracture opening, the greater the cross-sectional area. According to the law of mass conservation, the average volume flow rate at the water outlet is small when the water inlet flow rate is constant. For the same fracture opening, the velocity is larger in the middle position but smaller closer to the nonslip wall. The velocity will occur abnormal phenomenon when the local roughness is large; a partial enlargement of this phenomenon is shown in Figure 4(f), which presents a low velocity area at the corner and the surrounding velocity will increase sharply (Figure 4(e)).
Based on the mean pressure values along the inlet boundary and zero pressure along the outlet boundary of the five fracture models, the relationships between pressure gradient −∇P and the measured volume flow rate Q were plotted in Figures 5(a)-5(e), respectively. In order to study the nonlinear fluid flow characteristics in the low velocity region, the fluid flow behaviors can be fitted using the Forchheimer law, in which the pressure gradient is a quadratic function of the flow rate, written as follows: where a and b are the model coefficients, representing pressure drop components caused by linear and nonlinear effects, respectively. The pressure gradient refers to the pressure change per unit length along the fluid flow direction. The pressure on the outlet boundary was set as zero. The specific data obtained in the numerical simulation are shown in Table 2.
As d increases, both linear and nonlinear coefficients decrease, as shown in Figure 6. The linear coefficient a decreases rapidly in the d range of 0.1-0.2 mm but slowly for d larger than 0.2 mm. As d increases from 0.1 mm to 0.3 mm, the nonlinear coefficient b decreases rapidly and then stabilizes gradually. It can be seen that the nonlinear characteristics are more obvious when the fracture opening is small. At the same inlet velocity, the smaller the fracture opening, the larger the volume velocity of the fluid passing through the fracture. Even if the flow channel is rough, the flow characteristics are mainly linear when the flow rate is small, but the order of magnitude of the nonlinear coefficient is still very large, as shown in Table 3.

Flow
Characteristics through the Pore-Fracture Combination Model 5.1. Evolution of Pore Water Pressure and Flow Lines. Figures 7(a)-7(f) present the evolution of pore water pressure when the fracture permeability coefficient is greater than that of the matrix. In M1, when t = 0 s, the pore water pressure isosurface is relatively dense in the initial state, which is consistent with the actual situation (Figure 7(a)). When the initial water pressure is added, the water pressure on the facing surface is equivalent. The contour surfaces present a U shape at the fracture position along the x-axis direction (Figure 8(a)). The farther away from the water inlet position, the deeper the concave surface. This is because the permeability coefficient of the matrix is smaller than that of the fracture surface, and the fluid first passes through the fracture. The pore water pressure in the fracture will be greater than that in the matrix on the same cross section. Both flow lines and flow directions in the model are randomly distributed for the water pressure has not transferred to the interior of the matrix (Figure 9(a)). For t = 0:5 s, the fluid gradually flows into the inner part of the model for a certain distance, and the spacing of equal surfaces of the pore water pressure gradually increases (Figure 7(b)). The opening of the U-shaped isosurfaces of    (Figure 7(f)). There will be a slight change in the isosurfaces of pore water pressure at the rough-walled fracture positions (Figure 8(b)), for tortuosity of the rough fractures varies the flow direction. The flow lines in the matrix are already basically parallel, indicating that the fluid in the fracture and matrix is no longer transferred to each other at a stable fluid flowing state. Figures 7(g)-7(l) present the evolution of pore water pressure for M2, in which the permeability coefficient of the matrix is larger than that of the fracture, corresponding to the thermally treated fractured granite samples after 800°C [37]. The flow lines and flow directions of the whole sample are randomly distributed at t = 0 s (Figure 9(h)), which is similar to that in M1. The difference is that the U-shaped isosurface is in the opposite direction along the x-axis (Figure 8(c)), indicating that the fluid first penetrates from the matrix rather than the fracture. The matrix denotes a greater contribution to the fluid flow; thus, it will take less time to reach a steady state. In this case, due to the fact that the pore water pressure in the matrix is always larger than that of the fracture at the same cross section, the fluid diffusion will be transferred from the matrix into the fracture before reaching the steady state (Figures 9(i) and 9(j)). For t = 5 s (Figure 9(k)), fluid flowing has become stable, and the contour surfaces are perpendicular to the x-axis and the streamlines in the matrix are basically parallel.
Variations of the pore water pressure over time, taking the inlet water pressure of P = 2 × 10 5 Pa as an example, are shown in Figures 10(a) and 10(b). The locations of the measuring points are at an identical y of 5 mm, z of 0 mm, while various x values of 1, 5, 10, 15, and 19 mm, respectively. When the permeability coefficient of the matrix is less than that of the fracture (M1), the farther the measuring point from the water inlet, the slower the increase rate of pore water pressure with time, but the time required to achieve the steady state is basically the same for various measuring points. However, when the permeability coefficient of the matrix is larger than that of the fracture (M2), it is quicker to attain the steady state than that for M1 due to the fact that the matrix has a larger cross-sectional area than that of the fracture. The water head difference at the proximal end of the water inlet is gradually reduced and gradually increased away from the water inlet. Whether the fluid flowing process reaches the steady state can also be judged according to the changes of pore water pressure. The increase extent of pore water pressure for M1 and M2 is different, which is caused by fluid diffusion between the matrix and fracture, but the peak values after reaching the steady state are basically the same, indicating that the water pressure after steady state is related to the fluid flowing path, rather than the pore structures or fracture geometries. Figures 11(a) and 11(b) show the variations of pore water pressure along the flow length at different times. The measured cut-off line is arranged at the center position of the model, with the starting point and end point coordinates of (0, 5, 0) and (20, 5, 0), respectively. The pore water pressure at stable states is proportional to the length of the flow path, while independent on the physical properties of the matrix or fracture. However, before the stable state, a nonlinear relationship between the pore water pressure and flow path can be observed.

Evolution of Darcy's Velocity in Matrix and Fractures.
In the initial state, in M1, at the same time, the flow velocity in the fracture is obviously larger than that in the matrix, and the flow velocity decreases gradually along the flow path (Figure 12(a)). However, in M2, the flow velocity through the matrix is larger (Figure 12(c)), and the velocity near the fracture boundary is larger, which is due to fluid diffusion between the fracture and the matrix. Then, after the fluid flowing process reaches its steady state, the flow velocity in the matrix and fracture remains stable, but there will be obvious changes near the interface between the rough fracture and matrix. The flow rate in the fracture varies due to changes in the fracture geometry (Figures 13(a) and 13(c)). Due to high flow velocity in some local areas, eddy current occurs at this location (Figures 13(b) and 13(d)). Figure 14 shows the evolution of Darcy's velocity of the outlet boundary along the direction of fracture aperture thickness in the middle position after the steady state. Taking the inlet water pressure of 0.2 MPa as an example, it can be seen that the flow velocity at the outlet is not the same in the matrix and fracture. In M1, the flow velocity in the fracture is obviously larger than that in the matrix before reaching the steady state, and the velocity in the matrix decreases gradually from the middle position to both sides (Figure 14(a)). After reaching the steady state, the change of flow velocity at the interface between the fracture and matrix does not show a gentle variation. A U-shaped flow velocity distribution can be observed along the fracture aperture thickness. In M2 (Figure 14(b)), the flow velocity in the fracture is obviously smaller than that in the matrix, thus to produce a larger overall volume flow velocity of the matrix and shorter time to reach the steady state. When fluid flowing reaches the steady state, the fluid flow behaviors in the fracture and matrix are independent except the interface positions. Figure 14 can reveal the mechanism of water leakage in the underground engineering. When a subaqueous tunnel is built in the stratum with a low permeability coefficient, the difficulty of constructing the tunnel is actually the same as that of the mountain tunnel with rich water. Of course, the premise is to make a good advance geological prediction.
The contribution of the matrix and fracture to the overall flow velocity of the model is various. Taking the inlet pressure of 0.2 MPa as an example, in M1, the volume flow rate through the matrix Q m1 , fracture Q f1 , and whole model Q r1 is 3:108 × 10 −12 m 3 /s, 1:529 × 10 −11 m 3 /s, and 1:840 × 10 −11 m 3 /s, respectively. Q f1 accounts for 83% of the total volume flow rate, which is 4.92 times larger than that of Q m1 . Although the permeability coefficient of fracture is larger than that of the matrix, the volume flow rate in the matrix cannot be ignored. In M2, the volume flow rate through the matrix Q m2 , fracture Q f2 , and whole model Q r2 is 1:504 × 10 −10 m 3 /s, 2:440 × 10 −12 m 3 /s, and 2:440 × 10 −12 m 3 /s, respectively. In this case, the volume flow rate through the matrix is up to 98.40%, which is 61.64 times larger than that through the fracture.
To analyze the variations in flow velocity along the flow length, a straight line along the fracture length direction in the middle position was chosen as the measuring line, and the change process of flow velocity was quantitatively assessed, as shown in Figure 15. For t = 0 s, the flow velocity at the entrance of the fracture is larger, but the flow velocity is zero in the fracture for fluid does not begin to flow. The peak flow velocity decreases gradually with time. In M1, the peak Darcy's velocity decreases and tends to move towards the exit of the fracture. There are obvious changes in Darcy's velocity where the fracture surface fluctuates greatly, which further explains the occurrence of vortex or reflux of flow in rough fractures. As the velocity is stable, the Darcy's       velocity changes approximately to a straight line along the fracture length. From the local magnification, the Darcy's velocity after the steady state changes more violently along the fracture width. This is mainly due to the fact that the measured line passes through both the matrix and the fracture, which results in a sharp change in the flow rate at the interface position, characterized by a sudden decrease in the flow rate in M1, but a sudden increase in M2.
In order to further discuss the fluid flow characteristics of rock with combined fracture and matrix, different inlet water pressures were, respectively, selected in numerical simulations to analyze the relations between pressure gradient and volume flow rate in the two models. The selected parameters and calculated results are listed in Table 4. Figures 16(a) and 16(b) show the relations between the pressure gradient and volume flow rate of the outlet in M1 and M2, respectively. Two typical flow characteristics can be identified from the numerical simulation results [37]. In M1, the fluid flowing through the fracture accounts for a larger proportion of the overall flow capacity of the fractured model. Due to the fact that the cross-sectional area of the matrix is much larger than that of the fracture, the volume flow rate through the matrix cannot be ignored. In M2, the fluid flowing through the fracture accounts for a small proportion of the overall flow capacity, especially when the pressure gradient is small.

Conclusions
In this paper, the fluid flow process through single rough fracture models and pore-fracture combination models was, respectively, numerically analyzed using the COM-SOL Multiphysics finite element software. The fluid flowing behaviors through both the rough fracture and matrix, evolution characteristics of the flow streamlines, and volume flow rates were all investigated. The main conclusions are as follows: (1) For fluid flow through single fractures, the flow velocity closer to the fracture boundary is relatively small. Due to violent changes of the local fracture surface roughness, nonflowing area or the situation of reflux happens. The smaller the fracture aperture thickness,   19 Geofluids the more obvious the nonlinear characteristics of fluid flowing through the fracture for an identical flow velocity (2) The fluid flowing behavior is closely related to the pore water pressure, and the overall flow direction is vertical to the pressure isosurfaces. In the initial state, the fluid prefers to flow towards the medium with a relatively large permeability coefficient. After the fluid flow attains the steady state, the isosurfaces of pore water pressure are planes perpendicular to the fracture length direction, which is independent on physical properties of models (3) Due to various permeabilities of the fracture and matrix, fluid diffusion occurs. For fluid flowing in a stable state, the flow velocity in the fracture and matrix is significantly different, and fluid flowing through the two media is independent except for the interface position. When the permeability coefficient of the fracture is greater than that of the matrix, the flow capacity through the matrix cannot be ignored because the cross-sectional area of the matrix is much larger than that of the fracture

Data Availability
All raw data are available from the corresponding author if needed.

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