Study on the REV Size of Fractured Rock in the Non-Darcy Flow Based on the Dual-Porosity Model

For the problem of whether the representative elementary volume (REV) obtained in the Darcy flow is also applicable to the case of the non-Darcy flow, the study on the REV size within the non-Darcy flow is proposed tentatively. The concept of the REV in the non-Darcy flow is based on the definition of the REV. According to the determination of the REV in the Darcy flow, the intrinsic permeability k and non-Darcy coefficient β are used simultaneously for the determination of the REV in the non-Darcy flow. The pore pressure cohesive element (PPCE) is developed with the subroutine in ABAQUS. Then the simulation method of the Darcy and non-Darcy flow in the fractured rock mass is built using the PPCE. The proposed plan is examined through the comparison with existing research results. It is validated that this technic is efficient and accurate in simulating the Darcy and non-Darcy flow in the fractured rock mass. Combined with fracture networks generated byMonte Carlo Simulation technique, the PPCE is applied to the study on the REV size. Both conditions of the Darcy and non-Darcy flow are simulated for comparison. The simulation results of this model show that the REV of the non-Darcy flow is inconsistent with the REV of the Darcy flow, and the REV of the non-Darcy flow is more significant than the REV of the Darcy flow.The intrinsic permeability k tensors obtained in the Darcy flow and the non-Darcy flow are basically the same.


Introduction
Hydraulic property of the fractured rock mass plays a vital role in on-site engineering, like the underground basement, oil storage, CO 2 storage, and giant dams.The numerical simulation method can provide essential references for such projects, among which the Finite Element Method (FEM) is the most sophisticated numerical method.In the simulation of flow in the porous media, the equal intrinsic permeability should be determined first.However, as for the fractured porous media, before determining the equivalent intrinsic permeability, there is another critical parameter, the representative elementary volume (REV) to be determined.The REV is defined as the minimum volume of the sampling domain, beyond which the intrinsic permeability of the sampling domain remains essentially constant [1].Therefore, the determination of the REV size is extremely important and meaningful for understanding the hydraulic behavior of the fractured rock mass.
In 1982, Long et al. [2] first studied the DFN (Discrete Fracture Network), in which the scale can be regarded as porous media, and two proposed qualifications: (i) the REV exists; (ii) the equivalent hydraulic property at the certain REV can be approximated by an intrinsic permeability tensor.Oda [3,4] researched the fractured rock's REV, anisotropy, and inherent permeability tensor using numerical simulation.Min et al. [5] investigated the REV by the Monte Carlo Simulation technique method for the generation of the DFN based on the results of a site characterization of the field and proposed the "CV" values (coefficient of variation), which is defined as the ratio of standard deviation over the mean value, and the prediction error can be adopted for the determination of the REV.Baghbanan and Jing [6] investigated the intrinsic permeability of fractured rock mass with the constant aperture when considering the correlation between the disturbed fracture aperture and the trace length.The influence of different aperture distributions on the REV is studied by comparing.Min et al. [7] and Baghbanan and Jing [8] examined the effects of stress on the intrinsic permeability and fluid flow in fractured rock mass.The stress ratio was mainly considered for the study, and the results showed that it became more difficult to establish an equivalent tensor and the REV with increasing stress ratios.It is also found that there is a significant difference between correlated and noncorrelated apertures and the fracture length distribution [8].The studies on the REV in three dimensions have also been reported [9,10].The volumetric joint count calculation is applied to the determination of the REV in three dimensions [9].Thus, in general, there have been many studies on the REV and the intrinsic permeability tensor of the fractured rock mass.However, almost all the reviews of the REV above rely on two fundamental assumptions: (i) the rock matrix is impermeable, and the fluid flow can only occur through fractures.Therefore, the DFN analysis is widely applied; (ii) the fluid flow through the fractures obeys the cubic law.It is also important to take the intrinsic permeability of the rock matrix into consideration [11].Besides, naturally, the Darcy flow can occur at low flow rates and the fluid flow would enter the non-Darcy flow zone with the highpressure gradient unaffectedly [10,[12][13][14][15]. Thus, considering the intrinsic permeability of the rock matrix and the non-Darcy flow in the fracture is necessary for the further study on the REV of the fractured rock mass.
In 1960, Barenblatt et al. [16] first proposed the dualporosity model for the simulation of the fluid flow in the fractured rock mass.In the dual-porosity model, the fluid can both flow in the rock matrix and fracture.In addition, there is a fluid exchange between them, which is more realistic.Many achievements have been made [11,[16][17][18][19]. Ren et al. [20] extended a unified pipe network method (UPNM), which is a dual-porosity model, and applied it to on-site engineering.The results show that the UPNM has excellent advantages in simulating free surface flows in fractured rock slopes and oil underground storage projects with the water curtain system.Excellent reflection of the pore pressure field has been found by adopting the dual-porosity model [11,21].Chen et al. [22,23] proposed the composite element method (CEM) based on the dual-porosity model.The CEM is adopted to investigate the REV of the fractured rock mass and numerical simulation of the on-site project.However, at present, the study on the dual-porosity model is mostly based on self-development, and it is quite complicated to realize this point, which makes it difficult to become popular.Therefore, it becomes more necessary for the realization of the dual-porosity model in mature commercial software, and it can exert advantages by simplifying its application in on-site engineering and related researches.
Moreover, many studies on the non-Darcy flow in the fracture have been carried out.Currently, it is considered that inertial force mainly causes the non-Darcy flow.Many experiments include it [10,24,25], and it is found that the fluid would present the "weak inertia regime" before totally entering the fully developed turbulence from the Darcy flow.Moreover, it is mentioned that the loss of initial inertial effect mainly results from fracture geometries.As for such inertial losses, the Forchheimer's equation has been proved by many experiments, which describe such flow behavior very well.
Besides, as for the study of flow behavior in the single fracture or intersecting fractures, many works have been done on the factors like equivalent hydraulic apertures, roughness, and pressure gradients [26].Only a few works focused on the multifractures, especially the non-Darcy flow behavior in the fracture network and its REV size.Zimmerman et al. [27] measured and computed the intrinsic permeability of a natural sandstone with fractures.When the Reynolds number is larger than 20, both results of the experiment and numerical simulation exhibit a Forchheimer-type regime, in which the non-Darcy flow pressure drop is quadratic in the flow rate.Many non-Darcy flow experiments in the fracture network have also been conducted [28,29].It has been reported that when the hydraulic gradient is large enough, the inertial effect becomes nonnegligible and the nonlinear characteristic of the fluid flow should be taken into consideration necessarily [28].What is more, the study on the Darcy and non-Darcy flow based on the dual-porosity model is rare.The 2D models, based on the dual-porosity model and mathematic analysis, are established by taking the fracture as interfaces within the Forchheimer flow while the flow in the surrounding matrix is such that Darcy's law is adequate [30][31][32], but those studies only set few fractures in porous media.Generally, to take both the dual-porosity model and the non-Darcy flow behavior in the fracture network into consideration, the REV size should be further studied.
Throughout the previous literature, on one hand, most studies on the REV of the 2D fractured rock rely on the discrete model, and it is assumed that the rock matrix is impermeable and the fluid flow only occurs in the fracture.Considering that the pressure in the matrix would affect the velocity of flow in the breach, the REV size of fractured rock based on the dual-porosity model is investigated.Besides, the simulation method of the Darcy and non-Darcy flow in the fractured rock mass is built, which is conducive to promotion.On the other hand, almost all the studies on the REV involve the fluid flow in the fracture obeying the Cubic law, and as for the nonlinear flow, most studies focus on the single fracture and the fracture network.The survey on the REV in the non-Darcy flow has not been conducted in previous work.Based on the review of the non-linear flow behavior, the REV size of fractured rock in the non-Darcy flow is tentatively investigated.
In this paper, firstly, based on the dual-porosity model and the mathematics model [30], a simulation method is proposed with the pore pressure cohesive element (PPCE) in ABAQUS, which can simulate both the Darcy flow in the matrix and the non-Darcy flow in the rock fracture.According to the interpretation of the REV in the Darcy flow, the definition of the REV in the non-Darcy flow is given, when the concept of the REV in the non-Darcy flow is proposed based on the description of the REV.Then, focusing on the problem whether the representative elementary volume (REV) obtained in the Darcy flow is also suitable for the case of the non-Darcy flow, the REV size in the non-Darcy flow with the stochastic fracture network generated by the Monte Carlo method is investigated.In this paper, the REV of fracture rock in the Darcy flow is compared with that in the non-Darcy flow, when paying attention to the difference of REV sizes between the Darcy and non-Darcy flow.

The Simulation Theory and Method of the Darcy and Non-Darcy Flow Based on the Dual-Porosity Model
2.1.Darcy's Law and Forchheimer's Law.Considering that the velocity of the fluid in the porous rock matrix is much lower than in the fracture, the fluid flow in the matrix always obeys Darcy's law.For the single-phase flow and incompressible flow in porous media, Darcy's law can be expressed as (in the absence of gravity) where is the intrinsic permeability of porous media; ∇ [ML −2 T −2 ] is the pressure gradient in the porous media.
When the flow velocity is low in the fracture, the relationship between the volumetric flow rate and the pressure gradient is linear, and it can be described as the Cubic law: where  [L 3 T −1 ] is the volumetric flow rate or discharge;  [L] is the fracture aperture;  [L] is the aperture of the idealized parallel fracture, and it is the equivalent hydraulic aperture  ℎ for the rough fracture;  [ML −1 T −1 ] is the dynamic viscosity of the fluid;   [L 2 ] is the intrinsic permeability defined as   =  2 /12;  ℎ is the cross-section area.When the slope of the flow rate and the pressure gradient no longer coincide, the flow would enter into the non-Darcy regime, and its constitutive relationship can be described as Forchheimer's law: where ∇ [ML −2 T −2 ] is the pressure gradient along the flow direction;  [ML −5 T −1 ] and  [ML −8 ] are the coefficients, which describe the energy losses of the flow caused by viscosity and inertia, respectively. and  in (3) are commonly written as follows: where  [L −1 ] is called the non-Darcy coefficient, whose value is related to the aperture and roughness of the fracture.
Reynolds number Re is defined as the ratio of inertial force to viscous force, while Re can be calculated as follows: where  [ML −3 ] is the density of fluid.2.2.The Pore Pressure Cohesive Element.The pore pressure cohesive element (PPCE, 6-node element) is developed from a cohesive element (4-node element), which can simulate the fluid flow in fractured rock [33].Chen et al. [34] and Carrier and Granet [35] study hydraulic fractures with the PPCE, indicating the PPCE can simulate the fluid flow in fractured rock precisely.As shown in Figure 1, in the two-dimensional case, nodes 3, 4, 5, and 6 have not only a freedom degree of 8 (the default pore pressure of ABAQUS is 8), but also the freedom degree of 1 and 2 (the freedom degree of the displacement), whereas nodes 1 and 2 only have a freedom degree of 8.
In the simulation, nodes 1, 3, and 5 (2, 4, and 6) have different pore pressure values, but six nodes have the same coordinates when the displacement is excluded, so that the entire unit looks like a line unit.The fluid flow in the PPCE includes the tangential flow within the fracture and the normal flow across the fracture, as shown in Figure 1, while Figure 2 shows the meshing examples for 2D intersecting fractures with the PPCE in ABAQUS.
In ABAQUS, the default constitutive relationship of tangential flow is described as (2), and the constitutive relationship of normal flow is described as follows: where  [M −1 L 2 T] is the fluid leak-off coefficients at the top and bottom element surfaces;   and   [ML −1 T −2 ] are the midface pore pressure and the pore pressure on the sides, respectively.
In the flow exchange of fractures and the matrix, the entire basin must follow the law of mass conservation.As shown in Figure 3,  1 and  2 are porous media and  is the fracture.Vector  is a vector normal to .In  1 ,  2 , for single-phase incompressible fluids, it can be expressed as follows: In the flow exchange of fractures and the matrix, the law of mass conservation where −/2 V , dn, V , being the tangential flow rate in the fracture; v 1 ⋅ n 1 , v 2 ⋅ n 2 represent the flow rate into the fissure and out of the fracture, respectively.

The Analysis Method Based on the PPCE for Simulation of
Non-Darcy in Fractured Rock.For the standard calculation in linear FEM, when taking the theory of flow for example, the relationship can be written as follows: where {} is the volumetric flow rate of the node, [] is hydraulic conductivity of the matrix, and {} is pore pressure of the node.
If the relationship between the flow rate and pressure is linear, then the intrinsic permeability matrix is a constant matrix.However, for the non-Darcy flow, the problem is complicated.The Forchheimer flow rate  can be obtained by solving (4): Then, From ( 12), it is obvious that the hydraulic conductivity  is no longer a constant one as pore pressure changes.Figure 4 abstractly describes the nonlinear relationship between the flow rate and pressure and the "slope" of a certain point on the curve is the intrinsic permeability matrix [].There are many ways of solving this kind of nonlinear problem, such as the iterative method, the incremental method, and the incremental iteration method in FEM.
As for the non-Darcy flow, the constitutive default of tangential flow can be adapted with a subroutine (nodes 1 and 2 in Figure 1), and then the case of the non-Darcy flow can be simulated in fractured media.For a line element with 2 nodes, (10) can be written as follows: where  1 and  2 are the flow rate of the node, L is the length of the line element, K is hydraulic conductivity, and  1 and  2 are pore pressure of the node.In this study, the primary incremental method is used for the calculation of the non-Darcy flow.As shown in Figure 5, specific steps are as follows.Each increment (Δ  or Δ  ) is small enough in the calculation.At the first increment, the quadratic term in (3) is ignored, and then  −1 and  −1 can be obtained for the calculation of the intrinsic permeability matrix in the next increase.Repeat the above steps until the non-Darcy flow is steady.At the end of each increment, there is an error justification, Δ  < .If the error is beyond the threshold, the increment should be smaller for precision.

Validation.
To validate the correctness and effectiveness of the PPCE and its subroutine, several cases with the Darcy flow and the Forchheimer flow are calculated to compare with previous studies.

Fractured Media within the Darcy Flow in the Fracture.
In this case, the steady state of fractured media is simulated.The dimensions and boundaries are shown in Figure 6.The The discretization and pressure fields of this field obtained by ABAQUS-PPCE are compared with DFM (the discrete fracture model) and FM-Meshfree (fracture mapping with mesh-free), which are proposed by Lamb [36], as shown in Figures 7 and 8. To be specific, Figure 7 shows the approaches of four kinds of algorithms for the same model.The mesh nodes should be set on the fracture for the FM-Meshfree model while nodes and element edges need to be placed on the fracture for ABAQUS-PPCE and DFM.It can tell from Figure 8 that the pore pressure contours obtained by four algorithms are consistent, and Figure 9 shows the profile of pore pressure for a vertical section taken through the center of the fractured media for comparison.It can be found that ABAQUS-PPCE's pore pressure curve and other curves are consistent with the result calculated in Lamb [36], which can validate that ABAQUS-PPCE has its effectiveness in fractured media within the Darcy flow in the fracture simulation.

Free Surface Flow in the Fractured Rock Dam (The Darcy Flow).
In this case, ABAQUS-PPCE is used to simulate one of the most common problems in geotechnical engineering: unconfined flow analysis.The dam model's dimensions and water elevations of the upstream and the downstream are shown in Figure 10(a).To evaluate the effects of the fracture on free surface flow, a fracture with the aperture  = 0.001 m and the length  = 2.5 m is embedded in the dam as shown in Figure 10(b).The results obtained by ABAQUS-PPCE and Ren et al. [20] are compared in Figures 11 and 12.When the experiment result calculated by ABAQUS-PPCE is slightly higher than others (about 0.4 m), it can be found that the curves of free surface flow coincide with each other.Thus, it is reasonable to conclude that ABAQUS-PPCE is also correct and efficient in solving unconfined flow problems.

Fractured Media within Forchheimer
Flow in the Fracture.Frih et al. [30,32] validated the correctness of the interface model by being compared with the standard model, both in the case of the Darcy flow and in the case of the Forchheimer flow.The dimensions and boundaries are shown in Figure 13, where   = 10 −9 m 2 ,  = 0.01 m (ignoring cubic law),  = 10000 m −1 , and   = 10 −12 m 2 .The steady flow results of pore pressure and velocity simulated by ABAQUS-PPCE are shown in Figure 14.The comparison of tangential velocities between ABAQUS-PPCE and the interface model is given in Figure 15.According to the velocity along the fracture in Figure 15 and the recalculation equation ( 6), it can be seen that Re is about 3000.Combined with previous studies, the flow in this fracture is turbulent.It is clear that both the Darcy flow and Forchheimer flow of ABAQUS-PPCE results are close to the interface model.Besides, Figure 15 shows that the flow velocity in the fracture is affected by the pore pressure in porous media, which indicates it has the role of a barrier.Obviously, this situation cannot be achieved by the discrete media model, so the dual-media model has its advantages.
Above all, by comparing ABAQUS-PPCE's results with other existing examples, it can be found that ABAQUS-PPCE is effective and accurate when simulating the Darcy flow and the non-Darcy flow (Forchheimer flow) within the fractured porous media.

The Stochastic Fracture Network and Darcy and Non-Darcy Flow Simulations of Fractured Rock
In this section, a fundamental question in the field of flow in fractured rock, REV, and equivalent hydraulic properties is researched.Firstly, fracture networks are generated by Monte Carlo Simulation technique.Then REV and equivalent hydraulic properties are studies in the cases of Darcy and non-Darcy flows.Different criteria are adopted to determine whether there is inconsistency in the REV in the case of the Darcy and non-Darcy flow. 1 shows the necessary parameters of the fracture network for this study, and the intrinsic permeability of the rock matrix is 3 × 10 −15 m 2 .The program based on the Monte Carlo Simulation technique [37] is used for the generation, and the basic parameters are input.In order to eliminate the boundary effects, ten fracture networks which are large enough are first generated in the program, and then smaller fracture networks are cut by the boundaries of the "analysis domain," as shown in Figure 16.From each of the ten large "analysis domains," fracture network models are extracted with various sizes from 1 m × 1 m to 10 m × 10 m for calculation.In total, 100 models are generated and simulated.To justify whether the calculated intrinsic permeability has a tensor quality at a specific REV size, the fracture network models are rotated at an interval of 15 ∘ in the clockwise direction for the calculation.Figure 16 only shows the angles of 0 ∘ , 30 ∘ , 60 ∘ , and 90 ∘ models and totally 72 models are generated for this purpose.

Generation of the Stochastic Fracture Network. Table
The generated fracture network is imported into ABAQUS for modeling.The element type of the fracture is COH2D4P, and the rock matrix is CPE6MP, which is a

The Darcy Flow.
When calculating the fractured rock flow, one assumption is that the fracture is the parallel smooth plane, which means the roughness of the breach is ignored.Figure 18 shows the general boundary conditions for the calculation of two types of constant pressure and two types of linear pressures in each direction.Then the flow rates in the direction and -direction can be calculated, and the complete components of the 2D intrinsic permeability tensors of the models,   ,   ,   , and   can be obtained using the following:  where   is the flow rate crossing the boundary along with axis,  is the cross-section of the model,   is the intrinsic permeability tensor, and /  is pressure gradient along the  direction.
The values of 1/√  of each direction would be plotted on a polar diagram.If the values of 1/√  can be fitted into an ellipse, then an intrinsic permeability tensor at a given side length scale may exist.The purpose of checking the possible tensor quality, such as an average intrinsic permeability matrix   , was calculated at this range using the following (take  = 0 for reference axes): where  is the number of rotations,   and   are the directional cosines, and    is the calculated intrinsic permeability in each rotated model.experiments [15,24,38,39].Since all the splits in this paper are assumed to be smooth parallel plates with an aperture of 85 um (Table 1), all single breaches in the fracture network have the same coefficients  and , and therefore they have the same Forchheimer's equation.As this simple fracture model has a fixed aperture, to obtain the quadratic coefficient  and the non-Darcy coefficient  accurately, the method of using numerical experiments to determine the parameters is adopted in this paper.In this paper, CFD calculation software FLUENT is used for nonlinear simulation of the smooth parallel plate model with length L = 100 mm and aperture= 85 m.Non-Darcy coefficient  is obtained when the aperture  = 85 m, which is used to determine the REV size in the case of the non-Darcy flow.In the study, a rectangular grid is utilized to discretize the model.The unit length is 0.1 mm and the aperture is 8.5 um.The governing equation of hydrodynamics is the Navier-Stokes equation.For this paper, the nonlinear partial differential equations of the mass conservation equation (continuity equation) and the momentum conservation equation (momentum equation) of isothermal and incompressible Newtonian fluids can be written as

The Non
where  [T] is time.Due to the assumption that effects of gravity are neglected, the physical term is ignored in (17).The turbulence model is calculated using a Realizable - model.For the non-Darcy flow of the fractured rock mass, consistent with the non-Darcy flow of a single fracture, the slope between the flow rate and the pressure gradient decreases as the pressure gradient increases.Therefore, to determine the REV and non-Darcy coefficient tensors in the non-Darcy flow of the fractured rock mass, different hydraulic boundaries need to be applied in the same model.According to the research of Chen et al. [15] and Liu et al. [28], five different pressure gradients are applied (/  = 1 × 10 5 , 2 × 10 5 , 3 × 10 5 , 4×10 5 , and 5×10 5 ).Although the five pressure gradients can ensure that the flow regime in the fracture network is the non-Darcy flow in general, due to the complexity of the fracture network, different single fractures in the fracture network may have different flow regimes, and especially the connectivity of some branch fractures is poor (such as the fracture in the red circle in Figure 17).Since a section is in the matrix of the rock mass, the flow regime in the fracture should be the Darcy state.However, many scholars [40][41][42][43][44] think that (3) could probably be applied over the entire range of the flow rate/velocity, including that it reduces linear Darcy's flow at a sufficiently low flow rate/velocity.Therefore, the flow regime of all fractures in the fracture network can be described by Forchheimer's equation.
Ignoring the "preferential flow" effect [45], totally 500 models are generated.Then ( 18) is adopted to fit the simulation data, and, therefore, the equivalent intrinsic permeability tensor   and non-Darcy coefficient tensor   can be obtained.

𝜕𝑃 𝜕𝑥
where   is the non-Darcy coefficient tensor;  is a modulus of the flow rate throughout the entire region.The process of solving a model is as follows.Take the model and the boundary shown in Figure 18 as an example.According to the boundary shown in Figure 18, similar to the expansion of Darcy's equation ( 14), ( 18) can be expanded into four linear equations, as shown in the following: where , the pressure gradient of -direction is  1 / = ( 1 −  2 )/; the pressure gradient of -direction is  1 / = 0; the first two equations of ( 19) can be obtained; when there is pressure difference in -direction, the flow rate of , the pressure gradient of -direction is  1 / = 0; the pressure gradient of -direction is  1 / = ( 1 −  2 )/; then the last 2 equations are obtained.
In the simulation of two-dimensional flow, the side length of the square is , then  = .Thus, four equations can be obtained under a pressure gradient, but there are 4 inherent permeability tensors  and 4 non-Darcy coefficient tensor .In this paper, the fluid flow under five pressure gradients was simulated, forming 45 linear equations.Solve the overdetermined system of 8 coefficients, and then solve the overdetermined equations by the least square method, and the tensor of intrinsic permeability  and non-Darcy coefficient  of the model of the determined size and the rotation angle can be calculated.Bachmat [46], Hsu and Cheng [47], and Wang et al. [48] have proved and applied the theory of the tensor of non-Darcy coefficient .For the determination of the equivalent intrinsic permeability tensor and the intrinsic permeability ellipse, the same method in the Darcy flow is carried out.The only difference is that, as for the non-Darcy flow, the non-Darcy coefficient tensor also needs to be determined.First, the values of 1/√  of each direction would be plotted on a polar diagram.Then ( 20) is used to calculate the average non-Darcy coefficient   (take the  = 0 for reference axes): where    is the calculated non-Darcy coefficient in each rotated model.

Determination of the REV
In this section, firstly, the numerical results of the Darcy flow about size effect on the intrinsic permeability are presented.Then the statistical results of the non-Darcy flow are presented with size effect on inherent permeability and the non-Darcy coefficient.The graphic presentation of results is shown more clearly, with direct results of intrinsic permeability components (  ,   ,   , and   ) and non-Darcy coefficient components (  ,   ,   , and   ), including their simulation results, mean, maximum, and minimum.The coefficient of variation (CV) is used for the determination of the REV.For comparison of the REV scale, the values of "CV" with the increase of the model size are also plotted on the graph.Finally, for the evaluation of the intrinsic permeability tensor and the non-Darcy coefficient tensor, calculated and fitted ellipses are drawn at a certain model size.

Results of the Darcy Flow.
Figure 20 shows the pore pressure field at the model side length of 9 m and the rotation angle of 0 ∘ in the Darcy flow, from which the pore pressure variation of the flow in fractured rock can be obviously seen.Calculated values of the minimum, maximum, and mean of the intrinsic permeability components of 10 fracture network realizations of different sizes in the Darcy flow are shown in Figure 21.The results show that the range of variation of each intrinsic permeability components,   ,   ,   , and   , becomes smaller with the increase of the model side length, which indicates the existence of the REV.The "CV" values of the diagonal components of the intrinsic permeability tensor   and   are shown in Figure 22, which shows the values of "CV" both larger than 30% at the side length of 1 m, and it means that the difference of each model is somewhat big.As the side length of 6 m is exceeded,   and   never change significantly.The "CV" values are both below 20%.In addition to that, the "CV" values are both below 10% with the side length of 8 m.Besides, the values of   and   are almost close to zero with the increase in the model size, which indicates that a nearly symmetric intrinsic permeability matrix would present.Comparing with the ellipse at the side length of 1 m, an ideal ellipse can be reached at the side length of 8 m (as shown in Figure 27, the dashed line represents average 1/ √ , and the solid line represents the calculated 1/ √ ).

Results of the Non-Darcy
Flow. Figure 23 shows the curves of calculated results, Min, Max, and mean of intrinsic permeability components,   ,   ,   , and   , which are scattered when the side length is small.With the increase in the model side length, the intrinsic permeability components tend to be steady, indicating the existence of the REV in the non-Darcy flow.The non-Darcy coefficient component (  ,   ,   , and   ) variation curves are shown in Figure 24.
Like the intrinsic permeability component curves, the non-Darcy coefficient component curves tend to be steady with the increase in the side length.Figures 25 and 26 show the values of "CV" of   ,   ,   , and   with the increase in the model side length.As for   and   , at the side length of 1 m, the values of "CV" are beyond 30%.The values of "CV" tend to decrease with the increase of the model side length and are below 20% at the side length of 6 m.At the side length of 8 m, the values of "CV" are below 10%.As for   and   , the values of "CV" are both beyond 40% at the side length of 8 m.With the increase of the model side length, the values of "CV" decrease and are below 20% at the side length of 7 m.The values of "CV" are below 10% at the side length of 9 m.
The intrinsic permeability tensor ellipses of the non-Darcy flow are drawn in Figure 28, from which it can be found that the tendency of approaching an ellipse by the simulation data becomes stronger as the side length of the model increases.As for the non-Darcy coefficient tensor ellipse, an ideal ellipse can be achieved at the side length of 7 m as shown in Figure 29, which means that the non-Darcy coefficient has the quality of the tensor.

Scale Effect and the REV Size.
According to the calculated results above, the following conclusions are mainly drawn.(I) As for the Darcy flow, if an acceptable value of "CV" of 20% is set, then the side length of 6 m can be approximated for the REV size.If the value of "CV" of 10% is accepted, then the side length of 8 m can be approximated by the REV size.(II) As for the non-Darcy flow, two parameters (the intrinsic permeability  and the non-Darcy coefficient ) are necessary to determine the REV size.For the intrinsic permeability , if the value of "CV" of 20% is accepted, then the side length of 6 m can be approximated by the REV size.Besides, if the value of "CV" of 10% is accepted, then the side length of 8 m can be approximated by the REV size.For the non-Darcy coefficient , if the value of "CV" of 20% is accepted, then the side length of 7 m can be approximated by the REV size, and if the value of "CV" of 10% is accepted, then the side length of 9 m can be approximated by the REV size.
Besides, it is worth noting that, for the flow of Darcy and non-Darcy, the variation of the CV value of the equivalent permeability coefficient  is consistent, or, in other words, both flow regimes share the same REV.When the shapes of the intrinsic permeability ellipse are compared, they are almost identical.Thus, it is proved that the intrinsic permeability  is a physical quantity that reflects the inherent properties of fractured rock and has nothing to do with flow regimes.
Above all, what we can conclude is that if the value of "CV" of 20% is accepted, REV Darcy = 6 m and REV non-Darcy = 7 m.If the value of "CV" of 10% is accepted, REV Darcy = 8 m and REV non-Darcy = 9 m.Therefore, the REV size of the non-Darcy flow is a little larger than the REV size of the Darcy flow.

Discussions and Conclusions
Above all, the pore pressure cohesive element (PPCE) with a subroutine is developed, and the simulation of the non-Darcy flow in the fracture is realized.Then the simulation method of the Darcy and non-Darcy flow in the dual-porosity model, like fractured rock, is built.Main conclusions are as follows.
( simulation to validate its correctness and effectiveness.A new and more convenient method is provided to simulate the flow in the dual-porosity model. (II) The Monte Carlo Simulation technique is adopted for the generation of the fractured network.Based on the dual-porosity model with the PPCE, generated fracture networks with different side lengths and rotation angles are simulated in both Darcy and non-Darcy flow.On the one hand, the intrinsic permeability  calculated by Darcy's law and Forchheimer's law has identical convergence, and the determined REVs with the same criterion are identical.On the other hand, the value of "CV" of non-Darcy coefficient  has poor convergence, indicating that the REVs determined in the Darcy flow and the non-Darcy flow are inconsistent.Thus, the REV determined in the Darcy flow, the traditional way, fails to suit the case of the non-Darcy flow.
(III) According to the study on the determination of the REV in the Darcy flow and the non-Darcy flow, the REV determined in the non-Darcy flow is a little larger than that in the Darcy flow.Thus, it is appropriate to propose that the REV should be determined by the non-Darcy coefficient , which can satisfy not only the case of the Darcy flow, but also the case of the non-Darcy flow.
In this paper, based on the dual-medium model, the Monte Carlo fractured rock mass is studied tentatively using the non-Darcy coefficient to determine the REV value under the non-Darcy flow.The REV of the Darcy and non-Darcy flow is compared.Some new understandings have been    The inconsistency between the REVs mainly results from the REV of the non-Darcy flow determined by non-Darcy coefficient , and the REV of the Darcy flow is determined by inherent permeability , when the non-Darcy coefficient of porous media is considered to be only related to intrinsic permeability  [49][50][51].Then the REV determined by non-Darcy coefficient  and the REV determined by inherent permeability  should be consistent, but more studies have shown that the non-Darcy coefficient is not only related to the intrinsic permeability  of media.Evans and Civan [52], Geertsma [53], and Macdonald et al. [54] found that non-Darcy coefficient  is also related to the effective porosity of the pore medium; Li and Engler [55] and Veyskarami et al. [56] proposed that, for isotropic media, several correlations relate  to petrophysical parameters (in particular , k, and the tortuosity ); Knupp and Lage [57] also pointed out that  is related to the inherent permeability  and inertia coefficient   , where   represents the microform resistance exerted by the solid porous substrate; that is to say, it depends on the geometry of the solids in each basic representative volume of the porous medium.It is noteworthy that   and  are independent parameters.Therefore, it is not only related to the intrinsic permeability  but also related to other physical properties of the medium, which in turn leads to the inconsistency between the REV determined by  and the REV defined by , which may also be the reason for the disparity of the REV between the Darcy flow and non-Darcy flow in the fractured rock mass.
The conclusion that the REV of the non-Darcy flow is larger than that of the Darcy flow is only based on the results shown in the model calculation in this paper.Whether this conclusion is universal and whether this difference can be quantitatively described and which physical parameters of the fractured rock mass influence this outcome need to be further explored and studied in the future work.
Besides, it should also be mentioned that the PPCE has a wide range of applications in flow in fractured rock and soil, and the method proposed in this paper can also be applied in the simulation of on-site engineering.What is more, when it comes to the intersecting fractures in the non-Darcy flow, the "preference flow" is ignored in this study.It is suggested that similar studies involve this factor.For more work on the REV and its intrinsic permeability tensor of the fractured rock mass in the non-Darcy flow, factors like roughness, stress, and aperture also should be taken into consideration.

Figure 4 :
Figure 4: Relationship between the flow rate and pore pressure.

Figure 5 :Figure 6 :
Figure 5: Relationship between the flow rate and pore pressure in the incremental method.

Figure 9 :
Figure 9: The pressure profile for a vertical section taken through the center of the fractured media.

Figure 10 :
Figure 10: Dimensions and boundaries of the dam model.

Figure 12 :
Figure 12: Comparison of the results obtained by different methods.

Figure 13 :Figure 14 :Figure 15 :
Figure 13: Dimensions and boundaries of the model with one fracture.

Figure 16 :
Figure 16: Generation of the fracture network and the rotated models.

Figure 17 :
Figure 17: A certain fracture network model inside the length of 2 m and its meshes.
The left side of the model is set as the pressure inlet boundary, while the right side of the model is set as the pressure outlet boundary, and the exit boundary is set to zero pressure.The upper and lower boundaries of the fracture model are set as boundary conditions with no fluid flow and no slip.Water density and dynamic viscosity take the corresponding value of 20 ∘ C. Results are shown in Figure19.Fit the numerical data with (3), and  = 1.186913 kg/m −8 can be obtained with  2 = 0.999.Thus, (5) can be used to calculate , and  = 85.75 m −1 and all the non-Darcy flow parameters are determined.

Figure 20 :Figure 21 :
Figure 20: The pore pressure field at the side length of 9 m and the rotation angle of 0 ∘ .

Figure 22 :
Figure 22: The values of coefficients of variation of   and   with increasing side lengths in the Darcy flow.

Figure 23 :
Figure 23: The results, Max, Min, and mean of calculated intrinsic permeability components in the non-Darcy flow.(a)   , (b)   , (c)   , and   .

Figure 28 :
Figure 28: The equivalent intrinsic permeability ellipse with the increasing model size in the non-Darcy flow.

Figure 29 :
Figure 29: The equivalent non-Darcy coefficient ellipse with the increasing model size in the non-Darcy flow.

Table 1 :
Fracture parameters used for the generation of the fracture network.