Calculation method of equivalent permeability of dual-porosity media considering fractal characteristics and fracture stress sensitivity

In the petroleum industry, the accurate calculation of equivalent permeability of dual medium is very important for the estimation of reserves and the design of oil and gas production. At present, there are many methods to calculate the equivalent permeability of dual medium. These methods have their own advantages and disadvantages, and are suitable for different medium characteristics and physical problems, but they cannot calculate the equivalent permeability of dual medium under reservoir conditions. Thus, the emergence of fractal theory provides a theoretical basis for us to study the law of fluid transport in randomly distributed pores and fractures. Many scholars have deduced the analytical solution of the fractal model of double-porosity media, and the fractal dimension of the rock matrix and fracture network has been widely used in the permeability model of dual-porosity media. It has been proven that fracture permeability and pore matrix permeability considering fractal properties can be realistically applied. In this work, the fractal dimension of the rock matrix and fracture network aperture surface tortuosity are introduced to establish a dual-permeability calculation model considering the fracture closure effect. Taking the Yanchi Chang 8 area of the Ordos Basin as an example, comparing the results of the numerical simulation model and fractal model, the proposed two-porosity medium model considering the stress closure effect is more accurate. The sensitivity analysis of each parameter in the equivalent permeability model shows that considering the effect of stress closure reduced the fracture and pore aperture. The equivalent permeability is significantly influenced by the fracture inclination, the fractal dimension of the fracture aperture, and the fractal dimension of the matrix pore diameter. The rougher fracture surface and more tortuous capillary path forces particles to move longer distances, thereby reducing the permeability of the fracture network.


D T
Fractal dimension for the tortuosity of the matrix aperture (dimensionless), Dimensionless D x Fractal dimension for the tortuosity of the fracture aperture after σ eff is applied (dimensionless), Dimensionless E Young's modulus, GPa H f Fracture aperture, μm k Proportionality constant, Dimensionless K f Permeability of the fracture network, μm 2 K m Permeability of the rock matrix, μm 2 K t Equivalent permeability of dual-porosity media, μm 2 k x Fluid-solid coupling coefficient, Dimensionless L The length of the characterizing unit, μm L f The length of the fracture unit, μm L fo The initial length of the fracture unit, μm L 0 The initial length of the characterizing unit, μm L tf Tortuous length of a fracture, μm The total number of pores in porous media, Unit Q f The total flow rate in fracture network, m 3 Q m The total flow rate in porosity matrix, m 3 r Porosity matrix radius, μm r f Fracture radius, μm r fmin Minimum fracture radius, μm r fmax Maximum fracture radius, μm r f0min Minimum closure fracture radius, μm r f0max Maximum closure fracture radius, μm r min Minimum porosity matrix radius, μm r max Maximum porosity matrix radius, μm △p Differential pressure, MPa θ Fracture plane dip concerning flow direction in the front plane, °υ Poisson's ratio, Dimensionless

Introduction
The fracture shape in real reservoirs varies; that is, the distribution, length, orientation, and capillary radius of a fracture network in nature are usually disordered, there are few methods to find its analytical solution, and it is difficult to accurately characterize it in a dual-porosity media seepage model (Yu et al. 2017;Liu et al. 2018a, b;Zheng and Yu 2012;Li et al. 2016;Wei and Xia 2017;Berkowitz et al. 2022). Therefore, a large number of researchers believe that fractal geometry is an effective method to describe dual-porosity media (Lei et al. 2019;Xu et al. 2016;Luo et al. 2021;Sun et al. 2015). In addition, whether it is an artificial pressure fracture or natural fracture, under the action of the original stress, the fracture surface will reveal a closed state or tendency (Ghanbarian 2021; Jiang et al. 2013;Zhu and Cheng 2018;Zhang et al. 2021). With the continuous reduction in the pressure of the oil and gas reservoir, the effective closure stress of the fracture surface will also increase. Therefore, it is particularly important to calculate the equivalent permeability of fractures under different stress states (Cai et al. 2017;Xia et al. 2021Xia et al. , 2019Mandelbrot 2006;Sarkheil et al. 2012). Meanwhile, the spontaneous imbibition model in porous media and the model of unidirectional flow and multiphase flow in porous media with capillary forces have been solved by researchers (Sanei et al. 2015;Rigby et al. 2004;Xu et al. 2012). Later, Yu et al. (2003) proposed a new method of fractal-Monte Carlo simulation and then used this method to study the permeability of porous media. However, the spherical seepage in porous media is also significant (Yu et al. 2003;Hassani et al. 2011;Sarkheil et al. 2013a, b); thus, the spherical flow in porous media is implied by employing fractal theory, and the fractal solution of two-phase flow considering capillary force is obtained by Miao et al. (2018).
Most unconventional reservoirs have dual-media seepage problems. Therefore, many scholars have used fractal theory to describe the fracture network and pore structure of the porosity matrix and constructed expressions for the equivalent permeability of dual media (Chen et al. 2012;Miao et al. 2018;Liu et al. 2018a, b;Sarkheil et al. 2013a, b;Sarkheil et al. 2009). Chiles and Marisly (1993) proposed the relationship between the linear density, surface density, and volume density of fractures. Berkowiz and Hada (1997) believe that both natural and synthetic fracture networks have fractal characteristics, in which fracture surface density and fractal dimension are closely related. Meng et al. (2011) found that the surface density of fractures is positively correlated with the fractal dimension of fractures, and the surface density and fractal dimension of fractures jointly determine the development degree of fractures. Jafari and Babadagli (2011) obtained the fractal permeability expression of random fractures by multiple regression analysis based on logging observation data, which contained many empirical constants.
However, these expressions did not include the microstructure parameters of fracture tortuosity. Miao et al. (2015) obtained only the fractal permeability expression, including the tortuosity of the porosity matrix, according to the basic fractal theory of dual-porosity media. Therefore, to describe the pore structure of dual-porosity media more accurately, the fractal dimension of the tortuosity of the capillary wall in the fracture network should be introduced.
Although there are many results on the influence of fracture parameters on permeability, there is a lack of research on the quantitative influence of the fracture tortuosity fractal dimension and fluid-solid coupling coefficient on permeability. In this work, based on fractal geometry theory and models, a fractal model of equivalent permeability considering the tortuosity of fractures and the porosity matrix is established. The equivalent permeability of fractured porous media is calculated through the equation of fluid-solid coupling between fractures and the porosity matrix. The accuracy of the proposed model is verified by comparing the flow rate obtained by the model prediction with the numerical simulation results, and a sensitivity analysis of each parameter in the model is carried out.

Theories and methodologies
In this section, we focus on the derivation of permeability models of the porosity matrix and fracture network, considering the fracture closure effect and fractal dimension of fracture tortuosity. The basis for calculating the equivalent permeability of dual-porosity media is presented in the next section.

Fractal theory for a porosity matrix
Many researchers have reported the fractal theory of porous media, which assumes that the matrix porous medium is composed of a bundle of curved capillaries, and the cumulative number of capillary diameters greater than or equal to r in the porous medium satisfies the following fractal scaling law (Yu et al. 2003): The total number of pores in porous media is: In general, the number of capillary tubes or pores in porous media is very large, so Eq. (1) can be regarded as a continuous and differentiable function. By differentiating both sides of Eq. (1), we can obtain: The probability density function of the porosity size distribution can be obtained by dividing Eq. (3) by Eq. (2): The density function f(r) satisfies the normalization condition: If the above formula is true, the condition r min /r max < < 1 is satisfied, and the fractal dimension and ratio (r min /r max ) of the porosity of fractal porous media meet the following relation: The tortuosity of a curved flow path in porous media is defined as τ = L t /L 0 . The trace of fluid flow in a matrix curved capillary tub is usually curved and has fractal characteristics. The capillary diameter and length satisfy the following fractal scaling law: The total pore area representing the cross section of the unit is: The cross sectional area of the unit is: If the characterizing unit is assumed to be a cube and the length of the characterizing unit is L 0 = √ A m , according to Poiseuille's law, the flow rate in a single curved capillary ) is: Substituting Eq. (7) into Eq. (10) and integrating between the minimum and maximum pore sizes of all capillary tubes on the whole element body, we can obtain: Based on r min < < r max , then The permeability of the porosity matrix can be obtained from Darcy's law:

Fractal theory for a fracture network
Previous studies have shown that the fracture network can be analogized to countless capillary bundles that satisfy fractal scaling. According to basic fractal theory (Liu et al. 2018a, b), its fractal power-law relation is: Due to the large number of fractures in the fracture network, by differentiating Eq. (14), the number of fractures in the fracture network with capillary radii between [r, r + dr] can be obtained as: The fracture probability density formula is: where N t represents the total number of fractures, and the probability density of fractures is normalized as: When r fmin < < r fmax , Eq. (17) can be expressed as: When r fmin /r fmax ≤ 10 -2 , natural fracture networks generally meet the following requirement. The cross sectional area of all the capillaries in the fracture network ) is: Based on r min < < r max , then To accurately describe the change in degree of fracture opening during the development process, it is assumed that the bulge of the fracture surface will not rupture during the effect of closure stress. In this case, the capillary diameter considering closure stress can be written as: The streamline length increases with increasing effective stress, as obtained by Gere and Goodno (2012): The relationship between the velocity and the pressure gradient of a single capillary tube in the fracture can be obtained by the Navier-Stokes equation If the total flow rate of the fluid through a component fracture collection is needed, the flow rate of the fluid through a single fracture should be obtained first, and then the minimum radius to the maximum radius of the capillary in all fractures on the entire characteristic element cross section should be calculated.
In the two-dimensional plane, 0 < D f < 2, (r min /r max ) 3+D x −D f < < 1, and according to Darcy's law, the above equation can be simplified as: Substituting Eq. (25) into Eq. (30), and if the fracture plane dip concerning the flow direction is θ, Eq. (30) can be written as:

Equivalent permeability model of dual-porosity media based on fractal theory
According to the equivalent permeability of the porosity matrix and fracture network presented in the previous section, the flow rate expression of the porosity matrix and fracture network can be obtained as follows: Porosity matrix: Fracture network: The characteristic units of the fracture network are larger than those of the pore matrix and contain multiple matrix units. Assuming A f = n A m , Eq. (9) and (26) can be obtained as follows: Ignoring the interaction between the fluid in the pores of the porosity matrix and the fluid in the fracture network, the inflow and outflow of the fracture are equal, so the total flow rate of the medium with dual-porosity can be expressed as: The Newtonian fluid satisfies Darcy's law; then, the permeability expression of Newtonian fluid with dual porous media is: However, when the fluid-solid coupling between the fracture network and porosity matrix is considered, the large error of the calculated permeability is compared with the actual permeability, resulting in an inaccurate reservoir description. Therefore, the permeability coupling factor is proposed to characterize the fluid-solid coupling between the fracture network and porosity matrix, as shown in Eq. (38): Then, the total permeability of the dual medium system is: According to the lattice Boltzmann method, the equation of the permeability coupling factor is obtained by Yan (2019) Ultimately, the total permeability of the dual medium system is:

Results and discussion
Taking the Yanchi Chang 8 area of the Ordos Basin as an example, the equivalent permeability of double porous media is calculated. In this section, we use the Darcy seepage equation to calculate the flow rate. To verify the accuracy of the model, the model in this paper and the numerical simulation model were calculated. The parameters of the fractal model are shown in Table 1. Figure 1 shows the predictions by the present model compared with those of the numerical simulation model and the fractal model. From Fig. 1, the fractal model predictions do not fit the numerical simulation results. In contrast, the present model predictions agree well with the numerical simulation data. Among them, the numerical simulation model is simulated by tNavigator, a commercial software developed by the Russian RFD company. Because the pore structure of rock satisfies the fractal geometry theory, the fractal model is Yu (Yu et al. 2003) numerical model proposed by A, which adopts the fractal distribution of pores. In this model, the porous medium is regarded as a space composed of fractal geometry.
By deriving the permeability equation of fractured reservoirs considering fracture closure, a series of flow rate charts were calculated and plotted to conduct sensitivity analysis. Below, the influence of fracture network parameters on flow rate is presented first, and then the effects of some key parameters on flow rate, including the fractal dimension of fracture tortuosity, maximum radius of the fracture network, and the physical parameters of the fracture network, are further discussed. For this case, consider fracture closure by Young's modulus (E), Poisson's ratio (ν), and effective stress (σ eff ) on a fracture network. Second, influence of the porosity matrix parameters on the flow rate is discussed, including the fractal dimension of the porosity matrix, fractal dimension of the porosity matrix tortuosity, maximum capillary radius, and matrix porosity.  of fracture tortuosity has a larger impact on flow rate than the fractal dimension of fracture, and then the decreased range of flow rate decreases with the increase in the fractal dimension of fracture tortuosity and fracture. When the fractal dimension of the fracture tortuosity increased from 1.0 to 1.5, the flow rate decreased by 68.7%; when the fractal dimension of the fracture increased from 1.0 to 1.5, the flow rate decreased by 43.7%. For a fixed fractal dimension of fracture, Fig. 3 shows the effect of the maximum fracture radius on the flow rate. It can be seen in Fig. 3 that the fracture radius has a significant effect on the flow rate. The fracture radius is positively correlated with the flow rate, and an increase in flow rate can be seen as the fracture radius increases.

Influence of fracture parameters on the flow rate of the permeability model
Subsurface porous media rocks are composed of rock particles and pore fluids. During the development stage of the reservoir, compression of the rock will generate compound stress, which will lead to the deformation and displacement of the rock skeleton (Fig. 4). When the pore fluid migrates with the rock skeleton, relative seepage will also occur (Cao 2016;Wang et al. 2019).
As shown in Fig. 4, the reservoir rock can be equivalently regarded as a single hexahedral unit (Dai 2006;Feng et al. 2018). The opposite sides of the regular hexahedral elements are parallel to each other. Along the x, y, and z axes, the lengths of each side of the hexahedral unit correspond to dx, dy, and dz, respectively. The nine components on the reservoir rock are shown in Eq. (42). Moreover, the size of the nine components is related not only to the direction of the coordinate axis but also to the force of the point. The stress state at a certain point inside the rock can be represented in the form of a matrix: Chen and Ewing (Chen and Ewing 1999) used the elastic physical equation of porous media to describe the relationship between the effective stress and strain of rock, and its general expression is:  (Dai., 2006). Notes: σ x , σ y , σ z are normal stresses in the x, y, and z directions; τ xy and τ yx are shear stresses on the x and y planes; τ xz and τ zx are shear stresses on the x and z planes; τ yz and τ zy are shear stresses on the y and z planes where {σ eff } is the stress vector of the rock skeleton, m; and {ε} is the strain vector of the rock skeleton, m. Meanwhile, as displayed in Fig. 5, Young's modulus has a significant influence on the flow rate. The flow rate increases significantly under the same Poisson's ratio, and when Young's modulus is between 1 and 10 GPa, the range of the flow rate change is relatively obvious. As shown in Fig. 5, the effect of Poisson's ratio on the flow rate is negligible.
Similarly, Fig. 6 shows that the flow rate gradually increases with increasing fracture inclination angle and slowly decreases in the initial stage, but with increasing fracture inclination angle in the later stage, the flow rate reduction rate is accelerated. According to the result analysis presented in Fig. 6, in this model, when the fractal dimension of the fracture network is 1, there is no correlation between the fracture inclination angle and flow rate.
Influence of porosity matrix parameters on the flow rate of the permeability model Figure 7 displays the effect of the fractal dimension of the porosity matrix and porosity matrix tortuosity on the flow rate. The fractal dimension of the porosity matrix and porosity matrix tortuosity are negatively correlated with the flow rate. For a fixed fractal dimension of porosity matrix tortuosity, a decrease in flow rate can be observed as the fractal dimension of the porosity matrix increases. Meanwhile, the decreased range of the flow rate decreases with the increase in the fractal dimension of the porosity matrix tortuosity and porosity matrix. When the fractal dimension of the porosity matrix tortuosity increased from 1.0 to 1.5, the flow rate decreased by 0.9%; when the fractal dimension of the porosity matrix increased from 1.0 to 1.5, the flow rate decreased by 0.8%.
As shown in Figs. 8 and 9, the maximum porosity matrix radius and matrix porosity are positively correlated with the flow rate, but these two parameters have little impact on the flow rate.
Through comparison, we found that the fracture parameters have a significant influence on the porosity matrix parameters and the flow rate. In particular, the maximum   Fig. 7 Effect of the fractal dimension of the porosity matrix and porosity matrix tortuosity on the flow rate fracture radius, the fractal dimension of the fracture network, and fracture tortuosity have a significant influence on the flow rate, while the effect of the fractal dimension of the porosity matrix and porosity matrix tortuosity on the flow rate is weak. The difference between the fracture network and porosity matrix may demonstrate that the permeability of the fracture network is dominant in the equivalent permeability of dual-porosity media.

Conclusions
(1) The equivalent permeability of dual-porosity media can be described more accurately by considering the effect of stress closure and the fluid-solid coupling coefficient.
(2) The fractal dimension of the fracture diameter and tortuosity has a great effect on the equivalent permeability, and the fractal dimension of the fracture diameter and tortuosity is negatively correlated with the flow rate and equivalent permeability. (3) The fracture closure effect has a weak influence on the flow rate, while the maximum fracture radius has a significant effect on the flow rate. With the increase in the fractal dimension of the fracture surface tortuosity, the rising trend of the flow rate slows with the increase in the maximum fracture radius. (4) When the fracture inclination is less than 45°, it has a weak impact on the dual-porosity medium equivalent permeability; when the fracture inclination is greater than 45°, it has a significant impact on the dual-porosity medium equivalent permeability. (5) With the increase in the fracture inclination angle, there is a significant decrease in the flow rate under the same fracture dimension of fracture surface tortuosity. When the fractal dimension of the fracture network is 1, there is no correlation between the fracture inclination angle and flow rate. (6) The maximum matrix porosity radius and the matrix porosity are both positively correlated with the flow rate and equivalent permeability. The fractal dimensions of the matrix surface tortuosity and porosity matrix are both negatively correlated with the flow rate and equivalent permeability.

Conflict of interest
We declare that we have no financial or personal relationships with other people or organizations that can inappropriately influence our work. There is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in the manuscript. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.