A Two-Dimensional Partitioning of Fracture–Matrix Flow in Fractured Reservoir Rock Using a Dual-Porosity Percolation Model

: Fractures and micropores have varying contributions to the gas permeability of fractured reservoirs. The quantiﬁcation of the contribution of fractures and micropores that form a dual-porosity system for gas permeability is critical when attempting to accurately evaluate gas production. However, due to insufﬁcient knowledge of fracture–matrix ﬂow partitioning in such dual-porosity systems, it is challenging for previous models to quantitatively characterize the fracture heterogeneity and accurately evaluate the gas ﬂow and permeability in fractured rocks. In this study, we propose a dual-porosity percolation model to quantitatively investigate the contributions of fractures and matrix micropores towards the gas permeability of fractured rocks. Using percolation theory, we establish fracture networks with complex heterogeneity, which are characterized by various fracture densities and percolation probabilities within a porous matrix with various fracture/matrix permeability ratios. The compressible Navier–Stokes and Brinkman equations were adopted to describe the gas ﬂow in the fractures and porous matrix, respectively. The simulation results indicate that the gas permeability of the dual-porosity system has an exponential relationship with the fracture density and matrix permeability. The contribution of fractures and matrix micropores toward gas permeability can be classiﬁed by establishing a two-dimensional partitioning of the fracture–matrix ﬂow related to the fracture heterogeneity and fracture/matrix permeability ratio. The contribution of matrix micropores cannot be neglected if the fracture density is lower than a critical value.


Introduction
Unconventional hydrocarbon reservoirs are characterized by low porosity and permeability [1,2]. Fracturing stimulations create a complex fracture network that connects numerous micropores in rock matrices, which dramatically increases the gas permeability of unconventional tight hydrocarbon reservoirs [3,4]. In general, hydraulic fractures and micropores in rock matrices contribute differently to the gas permeability in fractured reservoir rocks. Previous studies have shown that gas flows in fractures and porous matrices behave differently and have effects on each other [5][6][7]. The quantification of fracture and micropore contributions in the formation of dual-porosity systems in terms of gas permeation capabilities is critical to accurately evaluate unconventional gas production [6]. However, due to insufficient knowledge of fracture-matrix flow partitioning in such dualporosity systems, previous studies have been unable to accurately quantify the gas flow and permeability of fractured rocks [5,8,9].
Hydraulic fractures connect the natural microfractures, micropores, and nanopores within a rock matrix in fractured reservoir rocks. They form the flow channels for fractured reservoir rocks after methane desorption from the organic nanopores and micropores [10][11][12]. Figure 1 shows a sketch map of the fractures and micropores within fractured reservoir rocks. Since hydraulic fractures are much larger than both micropores and nanopores in tight reservoir rocks, most of the previous studies have ignored the contributions of micropores to the gas permeability of reservoir rocks [13][14][15] and have mainly considered the fracture fluid flow when investigating the permeability of fractured rocks [16,17]. These studies assume that the reservoir medium is impermeable before fracture reconstruction, whereas fracture networks mainly dominate the permeability after fracturing. For nanopores and micropores, most previous studies have focused on the processes of gas adsorption, desorption, and diffusion in such structures [18][19][20][21]. However, gas flow in fractures and gas diffusion in matrices are different and have effects on each other in real fractured reservoirs. Rasmussen et al. [22] made several improvements by considering fluid flow in the rock matrix and discussed the effects of vertical fracture/matrix permeability ratios on fluid flow using the boundary element method (BEM). They considered fractures as a separate system from the matrices without taking into account the fluid flow from fractures to the matrix. Lough et al. [23] and Lee et al. [24] further improved the method suggested by Rasmussen et al. [22] and extended it to medium-sized fractures or larger fractures using the BEM. The influence of short fractures on matrix permeability was considered, while the effect of long fractures was not taken into account. Teimoori et al. [5] calculated the effective permeability in fractured reservoirs and concluded that the matrix permeability, fracture density, size, aperture, and interconnectivity all contribute to the effective permeability of an individual grid block. Recent studies [6,20,25] also suggested that these factors, especially the hydraulic fracture interconnections with micropores in the matrices, will directly affect the quantitative calculation of gas permeability of fractured rocks. Bai et al. [26] studied the dual-porosity behavior of naturally fractured reservoirs, taking into account the transient flow in the matrix blocks. Cai et al. [27] introduced a new modeling concept to numerically investigate the fracture-matrix interactions in shale gas reservoirs. The applications of new models are restricted to a few shape blocks and the universality of such models requires more tests and verifications. Wang et al. [28] formulated an analytical model to simulate real gas transport in nanopores and complex fractures in shale gas reservoirs and discussed the effect of multiple physics on the shale gas production, however the impact of fracture heterogeneity was not considered in their model. Abbasi et al. [29] proposed an analytical solution for fluid flow in a transient dual-porosity model and focused on investigating the influence of rock matrix block size in fractured formations. These studies attempted to simulate fracture-matrix interactions in tight fractured reservoirs from a variety of perspectives, but did not fully recognize or quantify fracture-matrix flow partitioning and its dependence on hydraulic fractures [3,30]. Nevertheless, Matthäi et al. [31] investigated numerical fluid flow partitioning between fractures and matrices based on field data and revealed the critical fracture aperture values that mark the transition from the matrix to fracture-dominated flow. Actually, the gas permeability in fractured rocks is not only related to the fracture aperture, but it intrinsically depends on the relative value of the fracture permeability to the matrix permeability. The three characteristic flow regimes were divided by their respective fracture-matrix permeability ratios (k f /k m ), named fracture-dominated flow, fracture-perturbed flow, and fracture-negligible flow. Sanaee et al. [32] also numerically simulated fracture core flooding test data to investigate the flow partitioning between fracture and matrix systems, which is affected by the in situ stress regime in the reservoir. The results also revealed a flow transition of higher fluid permeability under the lower stress magnitudes and a negligible impact of fracture distribution when the fractures are compressed under higher loading pressure. The respective fracture-matrix permeability ratio changes with the fracture closure due to the increase of loading pressure. Although the flow transition was observed in the abovementioned dual-porosity model with fractures and micropores, the partitioning of fracture-matrix flow is ambiguous and there is a lack of quantitative parameters to accurately distinguish the contributions of fractures and matrix flow. It has been demonstrated that reservoir rock has a complex heterogeneous and discontinuous structure, which plays a vital role in fluid flow through fractured geological media [33,34]. Therefore, it is inaccurate to describe and classify the flow transition in the fractured rocks using only one-dimensional partitioning of the fracture/matrix permeability ratio k f /k m . To intuitively characterize the impacts of discontinuity heterogeneity and reveal the controlling mechanism of fluid flow in porous fracture-matrix media, it is crucial to establish two-dimensional partitioning of the fracture-matrix flow regarding the contributions of the gas flow in the matrix. Additionally, previous studies have often ignored Klinkenberg effects and gas compressibility when calculating the gas permeability of fractured reservoirs [5]. abovementioned dual-porosity model with fractures and micropores, the partitioning of fracture-matrix flow is ambiguous and there is a lack of quantitative parameters to accurately distinguish the contributions of fractures and matrix flow. It has been demonstrated that reservoir rock has a complex heterogeneous and discontinuous structure, which plays a vital role in fluid flow through fractured geological media [33,34]. Therefore, it is inaccurate to describe and classify the flow transition in the fractured rocks using only one-dimensional partitioning of the fracture/matrix permeability ratio kf/km. To intuitively characterize the impacts of discontinuity heterogeneity and reveal the controlling mechanism of fluid flow in porous fracture-matrix media, it is crucial to establish two-dimensional partitioning of the fracture-matrix flow regarding the contributions of the gas flow in the matrix. Additionally, previous studies have often ignored Klinkenberg effects and gas compressibility when calculating the gas permeability of fractured reservoirs [5]. Since unconventional oil and gas reservoirs are often buried at depths of more than 3000 m [35], it is extremely difficult to detect and analyze gas flow in fractures and micropores in the matrix utilizing field tests. Additionally, it is also quite difficult to obtain rock core samples and measure gas flow during the fracturing process with exposure to real geological environments in the laboratory. Although the development of numerous laboratory methods has led to attempts to obtain the real dual-porosity structure in fractured reservoir rocks and discuss its influence on rock permeability, such as X-ray computed tomography (CT) [36][37][38] and scanning electron microscopy (SEM) techniques [39,40], it is nearly impossible to obtain a certain regular fracture system via experiments and quantitatively determine the statistical relationships between the pore space, fracture distribution, and equivalent permeability of fractured rocks using these methods. Another challenge to investigations of fractured reservoir permeability is that the fracture apertures and pore sizes are at various scales, which results in difficulties when attempting to quantitatively describe both the fluid flow in micropores and hydraulic fractures in fractured rock samples. Quantitative evaluation studies still face numerous obstacles, such as understanding the contributions of fractures and pores in matrices to gas flow at the core scale in the laboratory. Therefore, previous studies have adopted probabilistic methods and numerical simulation as alternatives, which are time-saving and flexible techniques that allow the construction of a pore-fracture system and an analysis of both the fluid flow in the pores and fractures of the fractured rock mass. Among them, the Monte Carlo method [41] is a widely used technique to simulate the fracture networks of rocks through a series of physical or geometrical statistical parameters. Combined with percolation theory, it is often used to investigate fluid flow through disordered porous media and tight reservoirs [42][43][44][45][46][47]. These studies mostly used simple line segments to obtain a graph of the fracture networks. The generation of most fracture networks requires large amounts of time, and few models accurately reconstruct the fracture while the heterogeneity, fracture length, fracture aperture, and random Since unconventional oil and gas reservoirs are often buried at depths of more than 3000 m [35], it is extremely difficult to detect and analyze gas flow in fractures and micropores in the matrix utilizing field tests. Additionally, it is also quite difficult to obtain rock core samples and measure gas flow during the fracturing process with exposure to real geological environments in the laboratory. Although the development of numerous laboratory methods has led to attempts to obtain the real dual-porosity structure in fractured reservoir rocks and discuss its influence on rock permeability, such as X-ray computed tomography (CT) [36][37][38] and scanning electron microscopy (SEM) techniques [39,40], it is nearly impossible to obtain a certain regular fracture system via experiments and quantitatively determine the statistical relationships between the pore space, fracture distribution, and equivalent permeability of fractured rocks using these methods. Another challenge to investigations of fractured reservoir permeability is that the fracture apertures and pore sizes are at various scales, which results in difficulties when attempting to quantitatively describe both the fluid flow in micropores and hydraulic fractures in fractured rock samples. Quantitative evaluation studies still face numerous obstacles, such as understanding the contributions of fractures and pores in matrices to gas flow at the core scale in the laboratory. Therefore, previous studies have adopted probabilistic methods and numerical simulation as alternatives, which are time-saving and flexible techniques that allow the construction of a pore-fracture system and an analysis of both the fluid flow in the pores and fractures of the fractured rock mass. Among them, the Monte Carlo method [41] is a widely used technique to simulate the fracture networks of rocks through a series of physical or geometrical statistical parameters. Combined with percolation theory, it is often used to investigate fluid flow through disordered porous media and tight reservoirs [42][43][44][45][46][47]. These studies mostly used simple line segments to obtain a graph of the fracture networks. The generation of most fracture networks requires large amounts of time, and few models accurately reconstruct the fracture while the heterogeneity, fracture length, fracture aperture, and random azimuth are simutaneously taken into account [25,45,48]. The complicated pattern of fracture networks is mainly determined by the number of fractures in a rock core, namely the fracture density. Zhang et al. [49] developed and extended the percolation theory to generate heterogeneous fracture networks with various fracture densities, fracture apertures, and random azimuths, which is applied in this study.
Intending to quantitatively investigate fracture-matrix flow partitioning in fractured rocks, in this study we developed a dual-porosity model to evaluate the gas permeability of tight rocks. The percolation theory was used to establish fracture networks with complex heterogeneity, which were characterized by various fracture densities and percolation probabilities. The facture networks were coupled with a porous matrix featuring various permeabilities, thereby leading to reconstructed structures with various fracture/matrix permeability ratios. We considered that gas flow not only occurs through connected fractures but also through the matrix micropores. Fractures and matrices were divided into two solution domains to simulate gas flow in fractured rocks using the compressible Navierstokes (N-S) and Brinkman equations, respectively. The coupling flow in the matrix and fractures was simulated using the variable exchanges between the velocity and pressure fields. We adapted the methods reported by Zhang et al. [49] to build five sets of random fracture models with varying fracture densities using Monte Carlo simulation technology. The fracture connectivity was quantitatively studied using the percolation theory. A singlephase flow simulation was performed to mimic methane flow and calculate the equivalent gas permeability values of fractured rock models, considering gas compressibility and Klinkenberg effects. The same fracture networks (but with three matrix properties) were studied to quantitatively analyze the impacts of the fracture density and fracture/matrix permeability ratio on the equivalent permeability of fractured rocks. The contributions of fractures and matrix micropores to the gas permeability were classified through twodimensional partitioning of fracture-matrix flow in porous media, considering the fracture heterogeneity and fracture/matrix permeability ratio. This study provides a modeling basis for evaluation of the permeability of unconventional oil and gas reservoirs.

Materials and Methods
This study used MATLAB to generate 2-D-connected fracture networks. The procedures included three steps: (1) five groups of fracture models with identical fracture densities but different fracture distributions were generated using the Monte Carlo method in MATLAB; (2) the fracture connectivity was analyzed based on percolation theory, whereby for fracture densities lower than the critical density, the permeable probability was zero; (3) for networks with a fracture density near or above the critical density, the system was considered to be permeable and we performed a compressible flow simulation to compute the methane flow in pores and fractures and determine the equivalent permeability according to Darcy's law.

Fracture Network Generation
Fluid flow in fractured rocks refers to the fluid flow in both porous (matrix and micropores) and fractured media (fractures and faults). Microporous structures have low hydraulic conductivity and their permeability is directly related to porosity. The fracture permeability not only depends on the fracture density but is also determined by the fracture distribution characteristics (i.e., fracture length, fracture apertures, and azimuths). In this study, we built five series of fracture models with varying fracture densities and distributions. While this approach can be expanded for applications to real geological formations, we made the following assumptions to construct the fracture systems: (1) Each system is a two-dimensional (2-D) square domain with sides of length L. Fluid flows into the system from the top of the system. (2) Fractures are randomly positioned. (3) Fractures follow a power law distribution for length, a normal distribution for aperture, and a random distribution for orientation. (4) Only connected fractures are considered. The details of the approach used to generate the fracture models are discussed in the following section.

Fracture Length Distribution
The geometric parameters for real reservoir fractures generally have statistical rules and obey one or several types of probability distribution. Referring to the work by Zhang et al. [49], the power law distribution was used to characterize the fracture length in this study, which is the most widely used model to describe the fault length distribution of reservoir rocks [50][51][52][53][54]. The power law distribution of a fracture system is expressed as: where n(l)dl refers to the fracture number that has a length interval [l, l+dl], α is a proportional coefficient that reflects the fracture density and depends on the fracture system size, and a is an exponent varying from one to three. The fracture density is defined as the cumulative fracture length per unit area in a fracture system, which can be can be calculated as [55]: where ρ means the fracture density (mm −1 ), n is the total number of fractures with a length of l in a fracture system, and L is the system size (mm).

Determining Hydraulic Connections for a Given Fracture Network
The connectivity of individual fractures determines numerous fundamental properties for a fracture network to determine system connectivity. Compared with previous studies, we only needed to perform the network generation process once and the fracture density of the connected networks could be accurately controlled to generate fracture networks with varying connectivity. The connected networks not only satisfied geometric continuity but also conformed to the average degree of connectivity. Here, a percolation parameter, p, which is related to the fracture density, is used as the critical value below which the fracture system is not connected, whereas the fracture system is, on average, connected, i.e., we also consider that the matrix is permeable. For a 2D square fracture system (size L) with N fractures of constant length l, p is calculated as: To generate a fracture network, we must determine three parameters to characterize the individual fracture geometry, i.e., the location, orientation, and length. The fracture center is randomly located in the system. The fracture orientation is also random, i.e., uniformly sampled from all directions. The fracture length (l) is sampled from its power law distribution. We use four geometry parameters to generate the fracture networks, i.e., the normalized system size, L s = L/l min ; the normalized maximum length of the fracture, L maxs = l max /l min (l min and l max refer to the minimum and maximum fracture lengths in a fracture system, respectively); the exponent a; and the ratio of r = a s /a sc . We refer to the fracture generation details reported by Zhang et al. [49]. For a given fracture system, the critical parameter, a sc , is calculated and compared with the actual parameter, a s . When a s is much less than a sc , i.e., r < 1, the system is viewed as non-connected. For systems with a s near or above a sc , i.e., r ≥ 1, the system is viewed as connected (the matrix is also considered to be permeable) and the steps described in the following section must be performed.
The fracture apertures in coal rock follow a normal distribution based on statistical data [9]. Based on their experiments, Wang et al. calculated the fracture apertures in coal rock, which ranged from 0.028 to 0.315 mm, with an average fractal aperture of 0.142 mm in a circular sectional area of 25 mm 2 . Therefore, in this study, we used an average fractal aperture, W av , of 0.284 mm for a square system with a size of L = 50 mm and a standard deviation of σ = 0.1 to generate the fracture networks. We controlled the fracture apertures between 0.03 and 0.63 mm. The fracture azimuths followed a random distribution, such that they were applicable to real geologic formations. The three fracture networks in each group had identical fracture densities. Table 1 lists the initial geometric parameters for the fracture networks. Figure 2 shows the fracture networks, while Figures 3 and 4 show the fracture aperture and model length distributions with various fracture densities, respectively. These geometry parameters for the fracture networks correspond to real geological cases and lay a good foundation for experimental studies in the next step to improve the predictive capability of the numerical methods developed in this study.     . Fracture length frequency distribution of models with various fracture densities. The horizontal axis represents the fracture aperture, while the vertical axis represents the fracture frequency at a certain length for models with various fracture densities. The data are the average results measured from three models at each fracture density.

Dual-Porosity Model for Gas Flow Simulation and Permeability ( ̅ ) Calculations
Here, we developed a 2-D dual-porosity model to simulate the gas flow in fractured rocks. The model treats matrix and pore structures as one component and couples the flow in matrix pores and fractures to compute the equivalent permeability. To study the influences of pore-fracture structures on the permeability properties of fractured rocks, we used an identical fracture network with three rock rank matrix properties from mining areas to simulate gas flow and compute permeability. The three ranks of rocks from mining areas, i.e., Irock, II rock, and IIIrock, were used to set the matrix properties. The properties of Irock, II rock, and IIIrock for the matrix permeability, , and porosity, , were taken from the sandstone used in our experiments, coal [56], and oil reservoirs [57], respectively.
Since fractures have much larger permeability than the matrices in tight reservoirs, solving such a nonlinear problem using traditional numerical methods poses difficulties. In this study, we divided the fractures and matrix into two domains to calculate the gas flow in a rock system. Variable exchanges between the velocity and pressure fields were used to simulate the coupling flows in the matrix and fractures. The compressible N-S equation was used to simulate methane flow in fractures [58] as follows: For the matrix, we used the compressible Brinkman equation, as follows: where p is the pressure, u is the flow velocity, μ is the dynamic viscosity, is the porosity, and is the matrix permeability. The gas flows into the domain from the top . Fracture length frequency distribution of models with various fracture densities. The horizontal axis represents the fracture aperture, while the vertical axis represents the fracture frequency at a certain length for models with various fracture densities. The data are the average results measured from three models at each fracture density.

Dual-Porosity Model for Gas Flow Simulation and Permeability (κ) Calculations
Here, we developed a 2-D dual-porosity model to simulate the gas flow in fractured rocks. The model treats matrix and pore structures as one component and couples the flow in matrix pores and fractures to compute the equivalent permeability. To study the influences of pore-fracture structures on the permeability properties of fractured rocks, we used an identical fracture network with three rock rank matrix properties from mining areas to simulate gas flow and compute permeability. The three ranks of rocks from mining areas, i.e., I rock , II rock, and III rock , were used to set the matrix properties. The properties of I rock , II rock , and III rock for the matrix permeability, k m , and porosity, ε m , were taken from the sandstone used in our experiments, coal [56], and oil reservoirs [57], respectively.
Since fractures have much larger permeability than the matrices in tight reservoirs, solving such a nonlinear problem using traditional numerical methods poses difficulties. In this study, we divided the fractures and matrix into two domains to calculate the gas flow in a rock system. Variable exchanges between the velocity and pressure fields were used to simulate the coupling flows in the matrix and fractures. The compressible N-S equation was used to simulate methane flow in fractures [58] as follows: For the matrix, we used the compressible Brinkman equation, as follows: where p is the pressure, u is the flow velocity, µ is the dynamic viscosity, ε p is the porosity, and k m is the matrix permeability. The gas flows into the domain from the top of the system. The fluid velocity was initially zero and flow was driven by a hydraulic head between the bottom (y = 0) and top (y = L) of the system. The velocities at the grain boundaries equaled zero, i.e., "no slip" conditions. We ensured that the flow was in a laminar regime by checking if the permeability remained constant. According to Darcy's law, the absolute permeability of the fractured rocks was computed using the following equation: where k is the absolute permeability, Q is the total fluid flux, L is the system length, A is the cross-sectional area of flow (i.e., here we use the system size L), ∆p is the pressure difference at the inlets and outlets, respectively; v i and s i are the section velocity and section area at the outlets, respectively. Due to Klinkenberg's effects, gas permeability is usually greater than liquid permeability in tight reservoirs [59]. This is especially the case when the pressure difference is low, and therefore the rock permeability measured by gas is much larger than that measured by water [60]. In this study, the following corrected formula was used to calculate the equivalent gas permeability κ g , which considers Klinkenberg effects based on previous studies [60][61][62]: where p av represents the average gas pressure and b refers to the coefficient of Klinkenberg effects, which is modified using the following equation: Figure 5 illustrates the domain settings and simulation boundary and Table 2 lists the model data used in this study. The fluid properties at normal temperature and pressure are used in the simulations.
Energies 2021, 14, x FOR PEER REVIEW 10 of 23 of the system. The fluid velocity was initially zero and flow was driven by a hydraulic head between the bottom (y = 0) and top (y = L) of the system. The velocities at the grain boundaries equaled zero, i.e., "no slip" conditions. We ensured that the flow was in a laminar regime by checking if the permeability remained constant. According to Darcy's law, the absolute permeability of the fractured rocks was computed using the following equation: where is the absolute permeability, Q is the total fluid flux, L is the system length, A is the cross-sectional area of flow (i.e., here we use the system size L), ∆ is the pressure difference at the inlets and outlets, respectively; and are the section velocity and section area at the outlets, respectively.
Due to Klinkenberg's effects, gas permeability is usually greater than liquid permeability in tight reservoirs [59]. This is especially the case when the pressure difference is low, and therefore the rock permeability measured by gas is much larger than that measured by water [60]. In this study, the following corrected formula was used to calculate the equivalent gas permeability ̅ , which considers Klinkenberg effects based on previous studies [60][61][62]: where represents the average gas pressure and b refers to the coefficient of Klinkenberg effects, which is modified using the following equation: Figure 5 illustrates the domain settings and simulation boundary and Table 2 lists the model data used in this study. The fluid properties at normal temperature and pressure are used in the simulations.

Dual-Porosity Model Verification
We performed flow simulations of a series of simple fracture networks to verify the accuracy of the dual-porosity model. We compared our numerical results with the theoretical solutions and numerical solutions from previous studies [63]. The previous studies assumed that the single fracture flow follows the cubic law and used the flow-equivalent principle to calculate the equivalent fracture permeability [64], which is expressed with the following equation: where l e is the element length of a square, w is the fracture aperture, and γ is the coefficient related to the roughness of the fracture surface, i.e., when the fracture is straight and smooth γ = 1/12, otherwise γ < 1/12. This formula reflects the effects that the fracture aperture, fracture surface roughness, and numerical method (element size) have on the equivalent permeability of the fractures. The verification fracture system uses an identical fracture structure to that reported in [63] to compare the numerical results with their numerical solutions. The domain is a square (20 m × 20 m), as shown in Figure 6a. For a single fracture network, we set nine sets of single fractures with different slopes ranging from 50 • to 90 • , with an interval of 5 • . A single fracture passes through the coordinate (0, 0) and the fracture aperture is equal to 5 mm.

Dual-Porosity Model Verification
We performed flow simulations of a series of simple fracture networks to verify the accuracy of the dual-porosity model. We compared our numerical results with the theoretical solutions and numerical solutions from previous studies [63]. The previous studies assumed that the single fracture flow follows the cubic law and used the flowequivalent principle to calculate the equivalent fracture permeability [64], which is expressed with the following equation: where is the element length of a square, w is the fracture aperture, and is the coefficient related to the roughness of the fracture surface, i.e., when the fracture is straight and smooth = 1/12, otherwise < 1/12. This formula reflects the effects that the fracture aperture, fracture surface roughness, and numerical method (element size) have on the equivalent permeability of the fractures. The verification fracture system uses an identical fracture structure to that reported in [63] to compare the numerical results with their numerical solutions. The domain is a square (20 m × 20 m), as shown in Figure 6a. For a single fracture network, we set nine sets of single fractures with different slopes ranging from 50° to 90°, with an interval of 5°. A single fracture passes through the coordinate (0, 0) and the fracture aperture is equal to 5 mm.
where is the water density (i.e., = 998.4 kg/m 3 ), is the gravitational acceleration (i.e., = 10 m/s 2 , μ is equal to 1.005 × 10 −3 Pa·s. We also assume that the fracture surface is smooth, with a value for of 1/12. Table 3 lists the simulation results. where ρ w is the water density (i.e., ρ w = 998.4 kg/m 3 ), g is the gravitational acceleration (i.e., g = 10 m/s 2 , µ is equal to 1.005 × 10 −3 Pa·s. We also assume that the fracture surface is smooth, with a value for γ of 1/12. Table 3 lists the simulation results.  Table 4 lists a comparison of our numerical results for nine single fractures with different dip angles with the analytical solutions. Our model's numerical solutions are slightly smaller than the theoretical solutions, whereby the deviation fluctuates within 5%. Table 5 provides a comparison of our numerical results with the corresponding results from the study by Wang et al. [63] and indicates that our calculation results are slightly smaller than the corresponding results reported in the latter, whereby deviations fluctuate within 6.5%. This is because different from the quadrilateral mesh used in the modified equivalent permeability model (MEPM) of Wang et al. [63], our model uses a free triangle mesh. The fracture and matrix were divided into two domains for meshing and the fracture area was meshed with a size ratio of fracture apertures to mesh elements that was larger than 0.5 [65] (see Figure 7). This type of meshing can ensure that each fracture is accurately characterized by the fractured mesh, such that the actual flow paths are very close to the fracture length. Although the meshing accuracy is directly related to the actual flow paths, results with different mesh accuracies all fluctuate around the theoretical solution, which indicates the accuracy of the proposed dual-porosity model.      Table 4 shows the intersecting fracture parameters and downstream outlet flow flux. The results demonstrate that the downstream outlet flow (0.1417 m 2 /s) is quite similar to the MPEM numerical solution (0.1424 m 2 /s) from the study by Wang et al. [63], as well as the theoretical solution (0.1426 m 2 /s). Our model can also achieve good simulation results for intersecting fractures and non-through fractures, verifying the validity of the dual-porosity model. Figure 8 demonstrates the COMSOL © Multiphysics results predicted using the methane velocity and pressure distribution values of the three rocks with fracture densities of 0.30/mm and 1.20/mm. We observed that as the fracture density increased, the fracture connectivity improved, the main flow paths increased, and the flow velocity gradually increased. The maximum flow velocity in the A1 model was 3.89 mm/s, while the flow velocity in the E model increased to 6.24 mm/s, i.e., 0.65-fold. The fracture number increased 3-fold and the main flow paths also significantly increased. The flow had clear hydraulic paths, which consisted of several connected fracture clusters. Less or even flow did not occur in the majority of connected fractures.   Table 4 shows the intersecting fracture parameters and downstream outlet flow flux. The results demonstrate that the downstream outlet flow (0.1417 m 2 /s) is quite similar to the MPEM numerical solution (0.1424 m 2 /s) from the study by Wang et al. [63], as well as the theoretical solution (0.1426 m 2 /s). Our model can also achieve good simulation results for intersecting fractures and non-through fractures, verifying the validity of the dual-porosity model. Figure 8 demonstrates the COMSOL © Multiphysics results predicted using the methane velocity and pressure distribution values of the three rocks with fracture densities of 0.30/mm and 1.20/mm. We observed that as the fracture density increased, the fracture connectivity improved, the main flow paths increased, and the flow velocity gradually increased. The maximum flow velocity in the A1 model was 3.89 mm/s, while the flow velocity in the E model increased to 6.24 mm/s, i.e., 0.65-fold. The fracture number increased 3-fold and the main flow paths also significantly increased. The flow had clear hydraulic paths, which consisted of several connected fracture clusters. Less or even flow did not occur in the majority of connected fractures. The fracture permeability values were much larger than the matrix perm values, i.e., 10 4 -10 9 -fold larger (see the velocity at the bottom of models A1 an Figure 9). For the low-density fracture model A1, the gas flow velocity and p gradient changed slightly with increasing matrix permeability. For the highfracture model E1, there were no significant changes in the gas flow velocity or p gradient. The matrix flow velocity of IIIrock increased, i.e., approximately 10 3larger than that of Irock and IIrock. The low fracture density model A1 had a small n of fracture outlets at the bottom and the matrix contribution to the rock permeabi relatively large. In contrast, the high fracture density model E1 had a large num fracture outlets at the bottom and the matrix made nearly no contribution permeability relative to fracture flow. We suggest that for rocks with low densities, the matrix will significantly contribute to rock permeability. These indicate that the proposed dual-porosity model is capable of quantitatively and v predicting the flow velocity paths and the pressure distribution of tight fracture r The fracture permeability values were much larger than the matrix permeability values, i.e., 10 4 -10 9 -fold larger (see the velocity at the bottom of models A1 and E1 in Figure 9). For the low-density fracture model A1, the gas flow velocity and pressure gradient changed slightly with increasing matrix permeability. For the high-density fracture model E1, there were no significant changes in the gas flow velocity or pressure gradient. The matrix flow velocity of III rock increased, i.e., approximately 10 3 -10 6 -fold larger than that of I rock and II rock . The low fracture density model A1 had a small number of fracture outlets at the bottom and the matrix contribution to the rock permeability was relatively large. In contrast, the high fracture density model E1 had a large number of fracture outlets at the bottom and the matrix made nearly no contribution to rock permeability relative to fracture flow. We suggest that for rocks with low fracture densities, the matrix will significantly contribute to rock permeability. These results indicate that the proposed dual-porosity model is capable of quantitatively and visually predicting the flow velocity paths and the pressure distribution of tight fracture rocks. Energies 2021, 14, x FOR PEER REVIEW

Estimation of Fractured Rock Permeability
Based on a qualitative analysis, flow paths increase with the fracture numb rock permeability increases as the fracture density becomes larger. Table 6 l simulation results for the equivalent permeability and total porosity.  Figure 10a,b plot the variations in equivalent permeability and total porosity w fracture density and different rock rank types, respectively. Here, the total porosit to the total volume fraction of pores and fractures in the fracture system. The perm and porosity of fractured rocks increased with increases in the fracture density, was consistent with field examples [66]. Based on Figure 10c, one can observe t permeability of the three rock ranks has the same variation trend as with an incr fracture density. Elevated fracture densities result in larger increases in the equ permeability with the growth of the fracture density. With a fracture densit 0.30/mm, the fracture density increases by 16.7% and the equivalent perm increases by approximately 24.3%. The fracture density increases by 3-fold at a (mm −1 ) and the equivalent permeability increases by approximately 8.73-fold.

Estimation of Fractured Rock Permeability
Based on a qualitative analysis, flow paths increase with the fracture number, and rock permeability increases as the fracture density becomes larger. Table 6 lists the simulation results for the equivalent permeability and total porosity.  Figure 10a,b plot the variations in equivalent permeability and total porosity with the fracture density and different rock rank types, respectively. Here, the total porosity refers to the total volume fraction of pores and fractures in the fracture system. The permeability and porosity of fractured rocks increased with increases in the fracture density, ρ, which was consistent with field examples [66]. Based on Figure 10c, one can observe that the permeability of the three rock ranks has the same variation trend as with an increase in fracture density. Elevated fracture densities result in larger increases in the equivalent permeability with the growth of the fracture density. With a fracture density ρ of 0.30/mm, the fracture density increases by 16.7% and the equivalent permeability increases by approximately 24.3%. The fracture density increases by 3-fold at a ρ of 1.2 (mm −1 ) and the equivalent permeability increases by approximately 8.73-fold. At the same time, we can observe that the equivalent permeability of fractured rock varies little with increases in the porosity and matrix permeability relative to th influences of the fracture density, which is consistent with field examples [66]. Howeve for low-density fracture networks, our results indicate that there is a considerable chang in the equivalent permeability due to matrix properties. Figure 10d compares the matr contribution ratios of flow with various fracture densities. For Irock and IIrock, there is a ver small matrix contribution, which is near zero, but the matrix contribution from IIIrock ca reach nearly 6.5% of the total flow. As the fracture density increases, the permeabilit which first increases, then gradually decreases with the matrix permeability. When th fracture density was 1.2/mm, the matrix contribution to rock permeability was nearly zer Our results also indicate that the matrix contribution in model A3 can reach 11.98% (se Figure 11). This is mainly due to the fact that despite the good fracture connectivity, ther are both fewer main flow paths and fractures at the bottom boundary (see Figure 12 which ultimately led to a weakened fracture flow and a relative increase in the matr contribution to the equivalent permeability. The matrix velocity (10 −3 m/s) of IIIrock was 3 5 orders of magnitude larger than that of Irock (10 −12 m/s) and IIrock (10 −9 m/s). Additionall the matrix contribution to the total flow, despite being approximately 4 orders o magnitude less than the fracture flow velocity (10 −3 m/s), was relatively significan Therefore, we can conclude that when the fracture/matrix permeability ratio, i.e., kf/km, less than or equal to 10 4 for low fracture density rocks (i.e., ≤ 0.35 (mm −1 ), wav ≤ 0.28 mm), the matrix can contribute significantly to the total flow. Here, fracture permeabilit Figure 10. Variations of κ g (a), ε p (b), ∆κ g, ρ /κ g,ρ=0.30 (c), and κ gm /κ g (d) with fracture density ρ for different rocks. The horizontal axis in each subfigure denotes the fracture density ρ, while the vertical axis represents the ratio of the gas permeability variation ∆κ g, ρ = κ g,ρ − κ g,ρ=0.30 in (c), which is the change of gas permeability with increasing fracture density to the value of rocks with ρ = 0.30, κ g,ρ=0.30 , i.e., ∆κ g, ρ /κ g,ρ=0.30 = (κ g,ρ − κ g,ρ=0.30 )/κ g,ρ=0.30 . The vertical axis represents the ratio of the matrix contribution κ gm /κ g in (d), which is the gas permeability of matrix domains κ gm to the total gas permeability of rocks κ g . Note: The numerical results at each point represent the average values calculated from three models with identical fracture density.
At the same time, we can observe that the equivalent permeability of fractured rocks varies little with increases in the porosity and matrix permeability relative to the influences of the fracture density, which is consistent with field examples [66]. However, for lowdensity fracture networks, our results indicate that there is a considerable change in the equivalent permeability due to matrix properties. Figure 10d compares the matrix contribution ratios of flow with various fracture densities. For I rock and II rock , there is a very small matrix contribution, which is near zero, but the matrix contribution from III rock can reach nearly 6.5% of the total flow. As the fracture density increases, the permeability, which first increases, then gradually decreases with the matrix permeability. When the fracture density was 1.2/mm, the matrix contribution to rock permeability was nearly zero. Our results also indicate that the matrix contribution in model A3 can reach 11.98% (see Figure 11). This is mainly due to the fact that despite the good fracture connectivity, there are both fewer main flow paths and fractures at the bottom boundary (see Figure 12), which ultimately led to a weakened fracture flow and a relative increase in the matrix contribution to the equivalent permeability. The matrix velocity (10 −3 m/s) of III rock was 3-5 orders of magnitude larger than that of I rock (10 −12 m/s) and II rock (10 −9 m/s). Additionally, the matrix contribution to the total flow, despite being approximately 4 orders of magnitude less than the fracture flow velocity (10 −3 m/s), was relatively significant. Therefore, we can conclude that when the fracture/matrix permeability ratio, i.e., k f /k m , is less than or equal to 10 4 for low fracture density rocks (i.e., ρ ≤ 0.35 (mm −1 ), w av ≤ 0.284 mm), the matrix can contribute significantly to the total flow. Here, fracture permeability refers to the permeability of a single fracture and is calculated as k f = w 2 /12, where w is the single fracture aperture.

The Quantitative Relationship between Pore-Fracture Distribution and Permeability
To better understand the impacts that coupled flow in a matrix and fractures has on rock permeability, we calculated the variations in permeability regarding the fracture density and matrix permeability to quantitatively analyze the permeability of fractured rocks. We found that the equivalent permeability ̅ of fractured rocks obeys the linear exponential function of the fracture density, ρ, and matrix permeability, , which is shown in the following equation: ̅ = −109.83 + 106.69 + 3.71 ; = 0.97 (16) Figure 13 illustrates the 3-D relationship between permeability, ̅ (μm 2 ), fracture density, (mm −1 ), and matrix permeability, km (μm 2 ). According to this relationship, for a given fractured rock structure, the permeability can be estimated by the statistical Figure 11. Variation of κ gm /κ g for rocks with low fracture densities: (a) ρ = 0.30 mm −1 and (b) ρ = 0.35 mm −1 . The horizontal axis indicates the matrix permeability k m , while the vertical axis represents the ratio of the matrix contribution κ gm /κ g , which is the gas permeability of matrix domains κ gm to the total gas permeability of rocks κ g .
refers to the permeability of a single fracture and is calculated as = /12, where the single fracture aperture.

The Quantitative Relationship between Pore-Fracture Distribution and Permeability
To better understand the impacts that coupled flow in a matrix and fractures ha rock permeability, we calculated the variations in permeability regarding the fra density and matrix permeability to quantitatively analyze the permeability of fract rocks. We found that the equivalent permeability ̅ of fractured rocks obeys the li exponential function of the fracture density, ρ, and matrix permeability, , whi shown in the following equation: ̅ = −109.83 + 106.69 + 3.71 ; = 0.97 Figure 13 illustrates the 3-D relationship between permeability, ̅ (μm 2 ), fra density, (mm −1 ), and matrix permeability, km (μm 2 ). According to this relationship a given fractured rock structure, the permeability can be estimated by the statis

The Quantitative Relationship between Pore-Fracture Distribution and Permeability
To better understand the impacts that coupled flow in a matrix and fractures has on rock permeability, we calculated the variations in permeability regarding the fracture density and matrix permeability to quantitatively analyze the permeability of fractured rocks. We found that the equivalent permeability κ g of fractured rocks obeys the linear exponential function of the fracture density, ρ, and matrix permeability, k m , which is shown in the following equation: κ g = −109.83 + 106.69exp ρ + 3.71exp k m ; R 2 = 0.97 (16) Figure 13 illustrates the 3-D relationship between permeability, κ g (µm 2 ), fracture density, ρ (mm −1 ), and matrix permeability, k m (µm 2 ). According to this relationship, for a given fractured rock structure, the permeability can be estimated by the statistical distribution characteristics of the fractures and matrix. This relationship provides data references when evaluating fracturing effects. distribution characteristics of the fractures and matrix. This relationship provides data references when evaluating fracturing effects. Figure 13. Equivalent permeability ̅ as a linear exponential function of fracture density and matrix permeability km. Note: The numerical results at each point represent the mean values calculated from three models with identical fracture densities.
Since the contribution of the matrix flow to the total flow is weaker than that of the fracture structures, we observed that the permeability increases with matrix permeability. However, using Equation (16), we speculated that the matrix contribution to the total flow increases with a decrease in the fracture density. When the fracture density, , is less than or equal to 0.30 and the ratio of the fracture-matrix permeability, kf/km, is less than or equal to 10 4 for the fracture system developed in this study (wav = 0.284 mm for the random fracture orientation, while the matrix permeability is between 0.01 and 1000 mD), the matrix contribution to the rock permeability can be significant. In contrast, when the fracture density, ρ, is greater than or equal to 1.2, the matrix contribution is nearly negligible.
Fluid flow partitioning between the matrix and fractures in fractured rocks is proposed in Figure 14. Compared with the study by Matthäi et al. [31], a two-dimensional partitioning of fracture-matrix flow concerning the fracture heterogeneity and the fracture-matrix permeability ratios (kf/km) was established for characteristic flow regimes. This indicated that critical behavior exists when the flow is either dominated by the matrix or fractures. Similar to the results of Matthäi et al. [31], a critical behavior was observed around kf/km ≈ 10 4 . For matrix permeability of 0.01 mD ≤ kf ≤ 1000 mD and a ratio of fracture-matrix permeability of kf/km ≤ 10 4 , the critical behavior occurs near = 0.35 mm −1 (wav = 0.284 mm). However, when the fracture density reaches a critical value, the flow is mainly dominated by fracture flow. This consists of three zones: (i) when the fracture density is extremely low (ρ ≤ 0.35 mm −1 , wav ≤ 0.284 mm), the matrix contributes significantly to the total flow; (ii) when the fracture density is low or medium, the matrix slightly affects the total flow; (iii) when the fracture density is sufficiently high ( ≥ 1.2 mm −1 , wav ≤ 0.284 mm), fracture flow completely dominates the total flow. Since the contribution of the matrix flow to the total flow is weaker than that of the fracture structures, we observed that the permeability increases with matrix permeability. However, using Equation (16), we speculated that the matrix contribution to the total flow increases with a decrease in the fracture density. When the fracture density, ρ, is less than or equal to 0.30 and the ratio of the fracture-matrix permeability, k f /k m , is less than or equal to 10 4 for the fracture system developed in this study (w av = 0.284 mm for the random fracture orientation, while the matrix permeability is between 0.01 and 1000 mD), the matrix contribution to the rock permeability can be significant. In contrast, when the fracture density, ρ, is greater than or equal to 1.2, the matrix contribution is nearly negligible.
Fluid flow partitioning between the matrix and fractures in fractured rocks is proposed in Figure 14. Compared with the study by Matthäi et al. [31], a two-dimensional partitioning of fracture-matrix flow concerning the fracture heterogeneity and the fracturematrix permeability ratios (k f /k m ) was established for characteristic flow regimes. This indicated that critical behavior exists when the flow is either dominated by the matrix or fractures. Similar to the results of Matthäi et al. [31], a critical behavior was observed around k f /k m ≈ 10 4 . For matrix permeability of 0.01 mD ≤ k f ≤ 1000 mD and a ratio of fracture-matrix permeability of k f /k m ≤ 10 4 , the critical behavior occurs near ρ = 0.35 mm −1 (w av = 0.284 mm). However, when the fracture density reaches a critical value, the flow is mainly dominated by fracture flow. This consists of three zones: (i) when the fracture density is extremely low (ρ ≤ 0.35 mm −1 , w av ≤ 0.284 mm), the matrix contributes significantly to the total flow; (ii) when the fracture density is low or medium, the matrix slightly affects the total flow; (iii) when the fracture density is sufficiently high (ρ ≥ 1.2 mm −1 , w av ≤ 0.284 mm), fracture flow completely dominates the total flow.

Conclusions
In this study, we developed a dual-porosity model to investigate the permeability of tight fractured rocks, with a special focus on the quantitative analysis of the coupling contributions of fractures and micropores to the equivalent permeability of fractured reservoirs using variable exchanges between the velocity and pressure fields. Percolation theory was used to establish fracture networks with complicated fracture heterogeneity via self-developed codes, characterized by fracture density and percolation possibility. The rock permeability calculations considered gas compressibility and Klinkenberg effects. The two-dimensional partitioning schemes for three flow regimes were proposed to further demonstrate the fracture-matrix flow partitioning in fractured reservoirs with various fracture densities and fracture/matrix permeability ratios. The main results of this study are as follows: (i) As the fracture density increases, the fracture connectivity gradually improves. The gas permeability significantly increases with an increase in the fracture density, which is due to increases in the fluid velocity and flow paths; (ii) The gas flow in fractured rocks has clear hydraulic paths, which consist of several connected fracture clusters. Lower or even flow did not occur in the majority of the connected fractures; (iii) Three flow regions were divided to demonstrate fracture-matrix flow using the two-dimensional partitioning method. The contribution of fractures and matrix micropores toward gas permeability can be determined using the fracture/matrix permeability ratios and the fracture density. The contribution of matrix micropores cannot be neglected if the fracture density is lower than a critical value; (iv) We derived the quantitative relationships between the fracture density, porosity, and equivalent permeability for the three ranks of rocks. Permeability is characterized by exponential growth with increases in fracture density and rock porosity. This study intends to provide a research basis and numerical reference for the quantitative analysis

Conclusions
In this study, we developed a dual-porosity model to investigate the permeability of tight fractured rocks, with a special focus on the quantitative analysis of the coupling contributions of fractures and micropores to the equivalent permeability of fractured reservoirs using variable exchanges between the velocity and pressure fields. Percolation theory was used to establish fracture networks with complicated fracture heterogeneity via self-developed codes, characterized by fracture density and percolation possibility. The rock permeability calculations considered gas compressibility and Klinkenberg effects. The two-dimensional partitioning schemes for three flow regimes were proposed to further demonstrate the fracture-matrix flow partitioning in fractured reservoirs with various fracture densities and fracture/matrix permeability ratios. The main results of this study are as follows: (i) As the fracture density increases, the fracture connectivity gradually improves. The gas permeability significantly increases with an increase in the fracture density, which is due to increases in the fluid velocity and flow paths; (ii) The gas flow in fractured rocks has clear hydraulic paths, which consist of several connected fracture clusters. Lower or even flow did not occur in the majority of the connected fractures; (iii) Three flow regions were divided to demonstrate fracture-matrix flow using the two-dimensional partitioning method. The contribution of fractures and matrix micropores toward gas permeability can be determined using the fracture/matrix permeability ratios and the fracture density. The contribution of matrix micropores cannot be neglected if the fracture density is lower than a critical value; (iv) We derived the quantitative relationships between the fracture density, porosity, and equivalent permeability for the three ranks of rocks. Permeability is characterized by exponential growth with increases in fracture density and rock porosity. This study intends to provide a research basis and numerical reference for the quantitative analysis and visual description of the relationships between rock permeability and pore-fracture structures.
It is noteworthy that although we were able to validate the dual-porosity model based on the tests of a series of single-fracture and cross-fracture models, there is still a lack of direct experimental verification of the proposed dual-porosity model. To overcome this, we intend to used three-dimensional printing (or additive manufacturing) technology [67] to fabricate a transparent dual-porosity model for future experimental investigations and to improve the predictive capability of the current numerical model. Moreover, because of the extremely high heterogeneity and anisotropy of the fracture structures in the real rocks, more factors need to be considered to better characterize the real rock structures and describe the 3D flow behavior in future studies and to explore the controlling mechanism of rock permeability in a dual-porosity model with micropores and fractures. Considering the primary purpose of this study and the limited length of research articles, we will discuss this issue in our follow-up studies.