The Criteria for Transition of Fluid to Nonlinear Flow for Fractured Rocks: The Role of Fracture Intersection and Aperture

: Understanding the ﬂuid pattern is of special signiﬁcance for estimating the hydraulic conductivity of fractured rock masses. The nonlinearity of ﬂuid ﬂow in discrete fracture networks (DFNs) originates from inertial effects and is enhanced by complex geometric topologies, which produces additional viscous friction and is subject to inertia effects, consequently transitioning the ﬂuid to the nonlinear ﬂow regime. Therefore, it is important to obtain the critical conditions for the transition of a ﬂuid from laminar to turbulent ﬂow. To investigate the role of fracture aperture and fracture intersection on the onset of the transition of a ﬂuid to nonlinear ﬂow in fractured rocks, the ﬂuid dynamic computation was performed by solving Navier–Stokes (N–S) equations in DFN models. The results show that the ﬂow ﬂux initially linearly correlates with the hydraulic gradient ( J ) and the permeability of DFNs initially remains constant. As the hydraulic gradient increases, the ﬂow ﬂux presents a strong nonlinear relationship with the hydraulic gradient, and the permeability decreases dramatically. In particular, signiﬁcant inertial effects appear earlier with a large fracture aperture or a dense fracture intersection. A critical hydraulic gradient ( J c ) is proposed to judge the onset of nonlinear ﬂow. The mathematical expression of J c and Forchheimer coefﬁcients A and B involving the fracture aperture and fracture intersection density is established through a multiple regression algorithm. Finally, the reliability of the predictive model was veriﬁed by comparing the results of the prediction and ﬂuid dynamic computation of a series of DFN models with well-known geometric distributions. The consistency of the ﬁtted equations and a correlation coefﬁcient greater than 0.9 between them indicate that the predictive model proposed in this study is reliable.


Introduction
The hydraulic properties of fractured rock have attracted significant concern in many geophysical processes and rock engineering applications, including water resource management [1][2][3], sustainable urban drainage [4][5][6], contaminant pollution control [7,8], hazardous waste isolation [9][10][11], and underground resource excavation [12][13][14][15].Since the permeability of fractures is usually several magnitudes higher than that of the rock matrix, the fracture is recognized as the predominant pathway of fluid flow, and the role of the matrix can be neglected in fluid flow analyses [16,17].Natural rock fractures control the fluid flow and mass transport in fractured rock and display complex hydraulic characteristics that come from the internal morphology of the fractures (i.e., the surface roughness and contact within individual fractures) and their arrangement to form anisotropic networks (i.e., the geometry of the fracture skeleton) [18,19].Due to the concealment of the underground rock mass, and the complexity of the spatial distribution and surface morphology of fractures, the minor structures of fractures are usually difficult to completely obtain through existing methods such as boreholes or exposure [20].Presently, the methodologies employed for generating Discrete Fracture Networks (DFNs) are commonly classified into three distinct categories: deterministic, stochastic, and hybrid.The constraints associated with acquiring geologic data pose significant challenges in attaining comprehensive information about fractures.Therefore, it is crucial to generate fractured rock models using stochastic methods based on limited geological information and to evaluate the linear and nonlinear flow characteristics of fluid flow in fracture networks through fluid dynamic computation [21][22][23].
Conventionally, the description of Newtonian fluid flow in fractured rock is given by the N-S equations, which are composed of a set of coupled nonlinear partial derivatives of varying orders [24].However, solving the N-S equation is difficult due to its complex nonlinear terms, and the existence and uniqueness of closed-form solutions of the complete N-S equations in three dimensions have not yet been proven [25].The simplified form that is derived from the N-S equation, such as the linear Stokes equation and the Reynolds lubrication equation, assumes that the inertia effects can be neglected and the hydraulic gradient is linearly proportional to the cube of fracture aperture, which is valid for laminar flow [26,27].In Darcy's flow, the linear correlation between the pressure drops and flow flux comes from the diminishing inertia effects (i.e., kinetic energy) that can only be satisfied for low flow rates.With the increase in flux, inertial terms cannot be ignored considering that the viscous forces gradually become stronger and the increase in hydraulic gradient is faster than flow flux, which is referred to as nonlinear fluid flow [28,29].Fluid flow behavior and the onset of nonlinear flow in rock fractures have been investigated empirically [30], experimentally [31], and numerically [32,33].Many efforts have been devoted to analyzing geometric parameters that affect the seepage regime of fractured rock masses, such as the joint roughness coefficient, tortuosity, fractal dimension, and contact [34,35].
The above studies represent a significant step forward in deepening our qualitative understanding of the nonlinear fluid flow and hydromechanical behavior of a single fracture [36][37][38] or intersection fracture [39].To date, the nonlinear fluid flow in fractured rock is still unclear.Therefore, substantial research has concentrated on nonlinear fluid flow behavior in fractured rock when the inertial effects become considerable [40,41].By investigating the role of geometric characteristics of the DFN model, including fracture length, fracture orientation, fracture surface roughness, equivalent hydraulic aperture, and connectivity in the transition of a fluid from laminar to turbulent flow, researchers found that the nonlinearity of fluid flow in the fractured rock becomes increasingly significant with increasing Reynolds number (Re) [42].Among them, fracture connectivity and aperture play a first-order effect in the transition of fluid behavior, and the roughness of the fracture surface would enhance the nonlinear fluid flow early [43].However, the role of fracture connectivity and aperture in the hydraulic properties of DFN models has not been quantitatively estimated, and the critical threshold of nonlinear fluid flow has not been systematically investigated.On the other hand, Re may not apply to DFNs, due to the flow in each fracture having a different Re, and understanding the Re in each fracture is challenging.In contrast, the hydraulic gradient (J), defined as the ratio of the hydraulic head difference to the DFN side length, is convenient for determining and describing the macroscopic flow characteristics.Therefore, it is essential to propose an expression to estimate the critical hydraulic gradient (J c ) for the onset of nonlinear flow in DFNs.
The fracture aperture and the connectivity of the fracture network are two crucial factors that control the permeability of the fractured rock [44,45].The utilization of numerical simulation has the benefits of cost-effectiveness and ease of execution on a wide scale in comparison to field trials and laboratory experiments.To investigate the role of the two factors mentioned above, a series of DFN models with segment lengths ranging from 30 mm to 60 mm with a gradient of 10 mm were constructed, and then fracture apertures between 0.5 mm and 1.5 mm with a gradient of 0.25 mm were assigned to each model.In total, fluid dynamics analysis was performed in 20 models (4 types of fractures segment length × 5 types of fractures aperture) using the laminar flow module of COMSOL Multi- By comparing the permeability obtained through the computation of the N-S equations and cubic law, the critical hydraulic gradient was proposed.Then, the expressions of J c and Forchheimer's coefficients were established through a multiple regression algorithm.Finally, a series of DFNs with well-known geometric distributions were constructed to verify the accuracy of the prediction model.

Governing Equations for Fluid Flow
The Navier-Stokes equations, which are equations of motion describing the conservation of momentum in viscous incompressible fluids, are used as the controlling equations in fluid dynamics computation.Due to the higher-order nonlinear terms in the NS equation, it is particularly difficult to solve it in complex fracture networks.Therefore, its simplified forms, the cubic law and Forchheimer's equation, are proposed for use as the governing equations for linear and nonlinear flows of fluids, respectively.

Navier-Stokes Equations
For incompressible Newtonian viscous fluids at the microscopic scale, the fluid flow is governed by the well-known N-S equations [46], which are written as where ρ is the fluid density, U is the velocity vector, P is the hydraulic pressure, µ is the viscous coefficient, and f is the body force.In Equation ( 1), the terms ρ ∂U ∂t represent force components caused by the rate of momentum change, and the terms ρU • ∇U represent force components caused by convective acceleration and are the only source of nonlinearity that is strongly dependent on the geometric variation of fractures.These two terms together represent the force from the inertial effect, and µ∇ 2 U represents the viscous force.

Cubic Law
The N-S equation is composed of a set of coupled nonlinear partial derivatives of varying orders [47], which makes it difficult to compute fluid flow in models with a large number of meshes.To overcome the drawback of difficulty in solving the N-S equation caused by nonlinear terms, empirical governing equations, e.g., the Stokes equation, Reynolds equation, and cubic law, have been proposed.Among them, the cubic law is the easiest to solve regardless of the complexity of the geometric model.In the cubic law, the fluid flux is proportional to the pressure gradient and the cube of the equivalent hydraulic aperture [48]: where Q is the flow flux; w is the width of the model, which is usually assumed to be 1 m in a two-dimensional model; h is the hydraulic aperture; and hydraulic gradient ∇P = dP/dx.Even though the cubic law was proposed based on some premise, multiple methods have been successfully proposed to theoretically prove the validity of the cubic law, such as flow experiments on fractured rock under triaxial conditions [49] and numerical simulations [50].

Forchheimer's Equation
The cubic law was simplified from N-S equations and was used to describe the laminar flow in a parallel-plate model at low velocity.For non-Dracy fluid flow, the cube law will bring nonnegligible deviations.Therefore, it needs to be described using other governing equations if the fluid is in the nonlinear flow regime.
The most extensively adopted mathematical description of the nonlinear flow relationships between the flow flux and pressure gradient in fractured and porous media is Forchheimer's equation, which is a zero-intercept quadratic equation to macroscopically describe nonlinear fluid flow and is described by where A and B are two coefficients correlated with the model geometries and represent the pressure drop due to linear and nonlinear effects, respectively.If the nonlinear effect is very small compared with the linear effect, the nonlinear term B Q 2 can be negligible, which indicates that Equation (3) can be seen as an improvement to Darcy's law by adding the nonlinear term.Hydraulic gradient J, a dimensionless parameter to characterize the hydraulic pressure in the flow system and defined by J = ∇P/(ρg), can be written as Equation (3): where AQ represents the term of linear fluid flow and BQ 2 represents the term of nonlinear fluid flow.It is notable that if the fluid rate is small enough, the nonlinear term can be neglected compared with the linear term, that is, Comparing Equations ( 2) and ( 5), it is found that Equation (4) can be regarded as adding a nonlinear term to Equation (2).Moreover, substituting Equation (2) into Equation (5) yields Equation (6) shows that the fracture aperture is the factor that determines the linear coefficient A.

Geometric Modeling and Fluid Dynamic Computation
Figure 1 shows the model with fracture segment lengths (L seg ) from 30 mm to 60 mm with a gradient of 10 mm, and this interval covers most of the distribution of fracture segments in fractured rock [32,33].These models have a size of 200 mm × 200 mm.Conventionally, the mechanical aperture of fractures and faults is usually in the range of 0.1 to 10 mm.The number of meshes exponentially increases with decreasing fracture aperture in a fixed discrete fracture network, which would greatly increase the difficulty of geometric discretization and fluid dynamic computation.In addition, fractures with a small aperture typically have small Re and flow rates, and the nonlinear characteristic would diminish with a small aperture.Considering the above, a fracture aperture from 0.5 mm to 1.5 mm with a gradient of 0.25 mm was selected and assigned to each fracture in the DFN models.In total, 20 DFN models (four types of fracture distance × five types of fracture aperture) were quantitatively analyzed.The fluid dynamic analyses were performed based on the FEM (finite element method) by directly solving the N-S equations in COMSOL.Achieving a balance between the amount and quality of meshes is a complex undertaking, and unstructured mesh dissection algorithms are capable of addressing this drawback.Skewness was employed as a The fluid dynamic analyses were performed based on the FEM (finite element method) by directly solving the N-S equations in COMSOL.Achieving a balance between the amount and quality of meshes is a complex undertaking, and unstructured mesh dissection algorithms are capable of addressing this drawback.Skewness was employed as a metric to assess the quality of the mesh, and the skewness values of the DFN model in relation to this study exceed 0.8 after geometry discretization.Therefore, the geometric model discretization is carried out using an unstructured grid, and then the N-S equations are solved at the grid nodes.The fluid enters from the left boundary and flows out from the right boundary, and the other two boundaries are impermeable.The hydraulic gradient applied to the boundary ranges from 10 −8 to 10 0 , which covers most of the possible transitions of a fluid in the DFN model from linear to nonlinear flow.

The Results of Fluid Dynamic Computation
Figure 2 shows the effect of the hydraulic gradient on the flow rate under different fracture apertures and fracture segment lengths.The flow rate increases quadratically with increasing pressure gradient.The data in Figure 2 were fitted using Forchheimer's equation, and the quadratic coefficient B and first-order coefficient A of Forchheimer's equation are listed in Table 1.The adequate acceptable consistency between the data and fitting curve demonstrates that the relationship between the flow rate and hydraulic gradient obtained using the N-S equation can be properly described using Forchheimer's law.With the increase in the equivalent hydraulic aperture, the flow area increases gradually, and the flow rate also continuously increases.Simultaneously, the curvature of the curve continues to increase, which indicates that the nonlinear relationship between the flow rate and hydraulic gradient strengthens.In addition, when the fracture intersection is dense, it is easier for the nonlinear fluid flow to occur in a lower hydraulic gradient, because fracture intersections can transform the direction of fluid flow, thus making it easier for the turbulent flow to appear.

Effect of h on A and B under Different Lseg
Figure 3 shows that with the increase in fracture aperture, coefficients A and B decrease rapidly first and then slowly, following a negative exponential trend.Equations ( 7) and ( 8) present the expressions of A and B established using a multiple regression algorithm based on the results of fluid dynamic computation:

Effect of h on A and B under Different L seg
Figure 3 shows that with the increase in fracture aperture, coefficients A and B decrease rapidly first and then slowly, following a negative exponential trend.Equations ( 7) and ( 8) present the expressions of A and B established using a multiple regression algorithm based on the results of fluid dynamic computation: where the unit of parameter It is noted that Equations ( 7) and ( 8) are just a quantitative relationship, and the dimensions are not considered.With the change in fracture aperture from 0.5 mm to 1.5 mm and the variance in fracture segment length from 30 mm to 60 mm, the variation in A and B is within an order of magnitude, and the value of B is two orders of magnitude higher than A. Equation (7) shows that the fracture aperture is the main factor affecting the linear coefficient A, which is consistent with Equation ( 6).The nonlinear coefficient correlates with the distribution of fractures, geometry of fracture intersections, fracture surface roughness, etc.This study focuses on investigating the influence of fracture aperture and fracture segment length, so Equation ( 8) has a similar form of Equation ( 7).The values of A and B calculated using the fluid dynamic computation and prediction model (Equations ( 7) and ( 8)) are listed in Table 1.Taking the values of fluid dynamic computation and prediction values as horizontal and vertical coordinates, respectively, Figure 4 compares the derivation between the fluid dynamic computation and prediction model for A and B. The data in Figure 4 are fitted by y = x, and the correlation coefficients are both greater than 0.99, verifying that the predicted model is reliable.
the distribution of fractures, geometry of fracture intersections, fracture surface roughness, etc.This study focuses on investigating the influence of fracture aperture and fracture segment length, so Equation ( 8) has a similar form of Equation (7).The values of A and B calculated using the fluid dynamic computation and prediction model (Equations ( 7) and ( 8)) are listed in Table 1.Taking the values of fluid dynamic computation and prediction values as horizontal and vertical coordinates, respectively, Figure 4 compares the derivation between the fluid dynamic computation and prediction model for A and B. The data in Figure 4 are fitted by y = x, and the correlation coefficients are both greater than 0.99, verifying that the predicted model is reliable.

Effect of Hydraulic Gradient J on Permeability k
The equivalent permeability of the DFN model can be obtained by back-calculating Darcy's law, which is written as where S is the cross-sectional area of the outlet, m 2 .When a small hydraulic gradient is applied to the DFN models, permeability k is constant, which indicates that the fluid flow conforms to laminar flow, BQ 2 = 0; the permeability is not correlated with the hydraulic gradient; and the flow rate can be calculated using the cubic law.As the hydraulic gradient J increases, a nonlinear fluid flow appears, the flow rate does not linearly increase with the increase in hydraulic gradient, and the permeability k decreases gradually.Forchheimer's equation needs to be considered to compute the fluid flow through DFN models.Figure 5 shows the effect of the hydraulic gradient J in the range from 10 −8 to 10 0 on the permeability k with different fracture apertures and fracture segment lengths.Taking the fracture aperture equal to 1.5 mm in Figure 5a as an example, when J < 10 −3 , k is a constant and the fluid conforms to a laminar flow; when J > 10 −1 , k dramatically decreases with increasing hydraulic gradient and a nonlinear fluid flow is predominant; when 10 −3 < J < 10 −1 , k slowly decreases with increasing hydraulic gradient and a nonlinear fluid flow begins to appear.

Effect of Hydraulic Gradient J on Permeability k
The equivalent permeability of the DFN model can be obtained by back-calculating Darcy's law, which is written as where S is the cross-sectional area of the outlet, m 2 .
When a small hydraulic gradient is applied to the DFN models, permeability k is constant, which indicates that the fluid flow conforms to laminar flow, 2 0 BQ = ; the permeability is not correlated with the hydraulic gradient; and the flow rate can be calculated using the cubic law.As the hydraulic gradient J increases, a nonlinear fluid flow appears, the flow rate does not linearly increase with the increase in hydraulic gradient, and the permeability k decreases gradually.Forchheimer's equation needs to be considered to compute the fluid flow through DFN models.Figure 5 shows the effect of the hydraulic gradient J in the range from 10 −8 to 10 0 on the permeability k with different fracture apertures and fracture segment lengths.Taking the fracture aperture equal to 1.5 mm in Figure 5a as an example, when , k slowly decreases with increasing hydraulic gradient and a nonlinear fluid flow begins to appear.
Figure 6 compares the permeability evolution versus hydraulic gradient under different fracture apertures and segment lengths.The permeability evolution is quantified by the ratio of permeability under arbitrary flow regimes and permeability in a laminar flow.In a small hydraulic gradient, permeability remains constant in each DFN model for various hydraulic apertures and segment lengths.As the hydraulic gradient increases to 10 −3 , the permeability of the model gradually decreases at different rates due to differences in fracture aperture and segment length of each model.With a large fracture aperture, the permeability obviously drops under the same degree of hydraulic gradient increase.The evolution of the permeability with fracture segment lengths of the DFN models is contrary to the fracture aperture.That is, the longer the fracture segment length is, the weaker the drop in model permeability with increasing hydraulic gradient.Figure 6 compares the permeability evolution versus hydraulic gradient under different fracture apertures and segment lengths.The permeability evolution is quantified by the ratio of permeability under arbitrary flow regimes and permeability in a laminar flow.In a small hydraulic gradient, permeability remains constant in each DFN model for various hydraulic apertures and segment lengths.As the hydraulic gradient increases to 10 −3 , the permeability of the model gradually decreases at different rates due to differences in fracture aperture and segment length of each model.With a large fracture aperture, the permeability obviously drops under the same degree of hydraulic gradient increase.The evolution of the permeability with fracture segment lengths of the DFN models is contrary to the fracture aperture.That is, the longer the fracture segment length is, the weaker the drop in model permeability with increasing hydraulic gradient.

Effect of Hydraulic Gradient J on Permeability Deviation δ
To quantitatively evaluate the effect of hydraulic gradient J on permeability evolution, permeability deviation δ is introduced, which is defined as Figure 7 shows the effect of the hydraulic gradient on permeability deviation under different fracture apertures and segment lengths.The permeability deviation δ increases with increasing hydraulic gradient.When the hydraulic gradient J is below 10 −3 , the various fracture apertures and segment lengths have no effect on the permeability deviation and δ = 0, which agrees with the conclusion obtained in the previous section.With J increasing to larger than 10 −3 , the fluid flow is in a nonlinear regime, a deviation in permeability appears, and δ increases with the increase in J.
A nonlinear flow requires solving N-S equations to obtain accurate fluid flow characteristics, which is a heavy task for numerical computation.To balance the computational accuracy and workload, a critical hydraulic gradient J c is essential as the criterion for judging whether the fluid regime transitions from laminar flow to turbulence. Figure 7 shows that J c is roughly in the range of 10 −3 to 10 −1 and decreases with increasing fracture aperture or decreasing fracture segment length.

A Prediction Model for J c
Conventionally, parameter E is applied to quantify the nonlinearity contribution in the fluid flow through DFNs, which is defined by Water 2023, 15, 4110 10 of 16 different fracture apertures and segment lengths.The permeability deviation δ increases with increasing hydraulic gradient.When the hydraulic gradient J is below 10 −3 , the various fracture apertures and segment lengths have no effect on the permeability deviation and 0 δ = , which agrees with the conclusion obtained in the previous section.With J increasing to larger than 10 −3 , the fluid flow is in a nonlinear regime, a deviation in permeability appears, and δ increases with the increase in J .A nonlinear flow requires solving N-S equations to obtain accurate fluid flow characteristics, which is a heavy task for numerical computation.To balance the computational accuracy and workload, a critical hydraulic gradient c J is essential as the criterion for judging whether the fluid regime transitions from laminar flow to turbulence. Figure 7 shows that c J is roughly in the range of 10 −3 to 10 −1 and decreases with increasing fracture aperture or decreasing fracture segment length.

A Prediction Model for c J
Conventionally, parameter E is applied to quantify the nonlinearity contribution in the fluid flow through DFNs, which is defined by An obvious drawback of this parameter is that it cannot explain the derivation of permeability obtained using the cubic law and the N-S equation.To determine whether the nonlinear flow needs to be considered, in other words, whether the deviation in permeability calculated using the cubic law and the N-S equation is considerable, the hydraulic gradient corresponding to δ equal to 10% is regarded as the threshold for nonlinear fluid flow and cannot be neglected.
Figure 8 summarizes the relationship between fracture aperture, segment length, and critical hydraulic gradient J c .With the increases in fracture aperture, the critical hydraulic gradient J c decreases following a power pattern.In contrast, the critical hydraulic gradient increased linearly as the fracture segment length increased.The prediction model of hydraulic gradient J c involved with fracture aperture and fracture segment length is constructed using multiple regression algorithms based on the results of dynamic analyses: The values of J c calculated through numerical simulation and prediction using Equation ( 14) are listed in Table 1.In addition, the simulation values and prediction values are taken as horizontal and vertical coordinates, respectively.The data in Table 1 are plotted in Figure 9 and fitted with y = x.The correlation coefficients are greater than 0.98, indicating that the prediction in Equation ( 12) is reliable.
of hydraulic gradient c J involved with fracture aperture and fracture segment length is constructed using multiple regression algorithms based on the results of dynamic analyses: The values of c J calculated through numerical simulation and prediction using Equation ( 14) are listed in Table 1.In addition, the simulation values and prediction values are taken as horizontal and vertical coordinates, respectively.The data in Table 1 are plotted in Figure 9 and fitted with y = x.The correlation coefficients are greater than 0.98 indicating that the prediction in Equation ( 12) is reliable.5.The Application of Prediction Formulas for A, B, and c J

Generation of DFNs
The prediction equations (Equations ( 7), (8), and ( 12)) are obtained through a specific geometric distribution.To verify its rationality, a series of discrete fracture networks with

The Application of Prediction Formulas for A, B, and J c 5.1. Generation of DFNs
The prediction equations (Equations ( 7), (8), and ( 12)) are obtained through a specific geometric distribution.To verify its rationality, a series of discrete fracture networks with a more realistic and well-known geometric distribution is generated.The geometric model was first generated using the Monte Carlo method, and then fluid dynamic computation was performed in COMSOL by directly computing the Navier-Stokes equations.
The power law is widely used to describe the distribution of fracture length, which is defined by [51] n where l is the fracture length in the range of [l min , l max ], n(l) • dl is the number of fracture lengths within [l, l + dl], D is the power exponent varying between 2.0 and 4.5, and α is the proportionality factor of the density.Geological surveys show that the orientation of fractures usually forms clusters in one or more directions, which can be described using the Fisher distribution [47].The deviation angle from the mean orientation θ is given by a cumulative probability density function with a Fisher constant κ: The fracture centers are assumed to follow a Poisson distribution [48].The locations of fracture centers are generated by a random process using the recursive algorithm that adopts the decimal part of numbers using the following equation: where R i is a random number in the range [0, 1], int(x) is the integer part of the number x, and an initial value R 0 is generated from the multiplicative congruence algorithm.If the generation space is along two coordinate ranges x g1 < x g2 and y g1 < y g2 in a Cartesian coordinate system, the midpoint coordinates (x i , y i ) of each fracture can be generated by In Figure 10, DFN models of five different fracture segment lengths are generated using the Monte Carlo method, and five types of apertures are assigned to these models.To diminish the effect of the randomness of fracture distribution, each DFN model is generated 10 times using 10 random numbers.In total, the fluid flow is performed in 250 models (5 types of L seg × 5 types of apertures × 10 times).
Water 2023, 15, x FOR PEER REVIEW 13 of 17 where i R is a random number in the range [0, 1], int(x) is the integer part of the number x, and an initial value 0 R is generated from the multiplicative congruence algorithm.If the generation space is along two coordinate ranges xg1 < xg2 and yg1 < yg2 in a Cartesian coordinate system, the midpoint coordinates (xi, yi) of each fracture can be generated by ( ) ( ) In Figure 10, DFN models of five different fracture segment lengths are generated using the Monte Carlo method, and five types of apertures are assigned to these models.To diminish the effect of the randomness of fracture distribution, each DFN model is generated 10 times using 10 random numbers.In total, the fluid flow is performed in 250 models (5 types of Lseg × 5 types of apertures × 10 times).

Comparison of Results of Fluid Dynamic Computation and Prediction Model
Figure 11 shows the results of fluid dynamic computation and the prediction model of A, B, and J c .Consistent with the results shown in Figures 3 and 8, the changes in A, B, and J c follow a similar trend with various fracture apertures under different fracture segment lengths.The averaged values of the fluid dynamic computation are close to the predicted results, and the similar fitting coefficients indicate that the A, B, and J c of the DFN model can be properly estimated using the prediction model.In addition, the maximum and minimum of the prediction values A, B, and J c are within one magnitude of the prediction values.With the increase in fracture aperture and/or the decrease in fracture segment length, the discreteness of the fluid dynamic computation values becomes weaker, and the predicted values are closer to the fluid dynamic computation values.Overall, the prediction model is more accurate for the A, B, and J c estimation of discrete fracture networks with smaller fracture segment lengths and larger fracture apertures.

Conclusions
This study established a series of DFN models with different fracture apertures and segment lengths, and a fluid flow was performed by directly computing N-S equations in COMSOL.The influence of the hydraulic gradient on the fluid flow characteristics of fractured rock was analyzed, and the permeability was calculated using Darcy's law for comparison.Based on the study of applied hydraulic gradients between 10 −8 and 10 0 on the DFN models, the following main conclusions can be obtained:

Conclusions
This study established a series of DFN models with different fracture apertures and segment lengths, and a fluid flow was performed by directly computing N-S equations in COMSOL.The influence of the hydraulic gradient on the fluid flow characteristics of fractured rock was analyzed, and the permeability was calculated using Darcy's law for comparison.Based on the study of applied hydraulic gradients between 10 −8 and 10 0 on the DFN models, the following main conclusions can be obtained: (1) When a lower hydraulic gradient is applied to the DFN model, the fluid flow conforms to laminar flow.The flow rate Q is linearly correlated with the hydraulic gradient J, which ensures that the permeability k remains constant.The permeability of fractured rocks can be accurately estimated by the cubic law.(2) With the increase in hydraulic gradient, the flow flux variance nonlinearly changes with the evolution of the hydraulic gradient.The permeability k sharply decreases with the increase in the hydraulic gradient, and the cubic law is not suitable for the computation of fluid flow through DFN models because the permeability deviation caused by the pressure loss would be underestimated.(3) With the increase in fracture aperture or the decrease in segment length, the fluid flow transition from the linear regime to the nonlinear regime in DFN models would occur at a lower J c .Therefore, a critical hydraulic gradient J c is proposed, below which the cubic law can be used to compute fluid flow through DFN models and above which Forchheimer's law should be applied.(4) Mathematical expressions for prediction coefficients A and B and critical hydraulic gradient J c are proposed using a multiple regression algorithm.To verify the reliability of the proposed expressions, a set of discrete fracture network modes with well-known geometric distributions are generated.The coefficients A, B, and critical hydraulic gradient J c of these models are calculated through fluid dynamic computation and prediction models, and a correlation coefficient R 2 larger than 0.9 shows that the prediction model can provide a reasonable prediction of A, B, and J c .
This work focuses on the influences of fracture intersections and the aperture on the onset of nonlinear fluid flow in DFN models.It should be noted that other factors, such as fracture surface roughness and fracture number, can also strongly affect the flow nonlinearity in DFNs.More extensive studies are required to investigate the flow behavior of fractured rocks considering these factors and to derive more generous expressions for critical hydraulic gradients J c and Forchheimer's coefficients A and B. Furthermore, the conclusion is derived from the two-dimensional discrete fracture network (DFN) model, which may not align with the actual scenario.This aspect presents an intriguing area for exploration in future research endeavors.

6 (
Chongqing, China) to understand the fluid flow mechanisms in fractured rock.

Figure 3 .
Figure 3. Forchheimer coefficients as a function of fracture aperture and fracture segment lengths: (a) linear term coefficient A; (b) nonlinear term coefficient B.

Figure 3 .
Figure 3. Forchheimer coefficients as a function of fracture aperture and fracture segment lengths: (a) linear term coefficient A; (b) nonlinear term coefficient B.

Water 2023 , 17 Figure 4 .
Figure 4.The comparison of fluid dynamic computation and prediction values of (a) linear term coefficient A and (b) nonlinear term coefficient B.

,,
k is a constant and the fluid conforms to a laminar flow; k dramatically decreases with increasing hydraulic gradient and a nonlinear fluid flow is predominant; when

Figure 4 .
Figure 4.The comparison of fluid dynamic computation and prediction values of (a) linear term coefficient A and (b) nonlinear term coefficient B.

Figure 6 .
Figure 6.Effect of pressure gradient on the decay of permeability under different fracture segment lengths and different fracture apertures: (a) L seg = 30 mm; (b) L seg = 40 mm; (c) L seg = 50 mm; (d) L seg = 60 mm.

Figure 8 .
Figure 8. Evolution of the critical hydraulic gradient in terms of fracture aperture and fracture seg ment length.

Figure 8 . 17 Figure 9 .
Figure 8. Evolution of the critical hydraulic gradient in terms of fracture aperture and fracture segment length.Water 2023, 15, x FOR PEER REVIEW 12 of 17

Figure 9 .
Figure 9.The simulation and prediction values of the critical hydraulic gradient J c .

Figure 10 .
Figure 10.The model generated by the Monte Carlo method with different Lseg.

Figure 11
Figure 11 shows the results of fluid dynamic computation and the prediction model of A, B, and c J .Consistent with the results shown in Figures 3 and 8, the changes in A, B, and c J follow a similar trend with various fracture apertures under different fracture segment lengths.The averaged values of the fluid dynamic computation are close to the

Figure 10 .
Figure 10.The model generated by the Monte Carlo method with different L seg .

Water 2023 , 17 Figure 11 .
Figure 11.The fluid dynamic computation, prediction, and average fluid dynamic computation values of A, B, and c J of DFN models.

Figure 11 .
Figure 11.The fluid dynamic computation, prediction, and average fluid dynamic computation values of A, B, and J c of DFN models.

Table 1 .
Geometric parameters of the generated DFNs, fluid dynamic computation (N), and predicted results of A and B.