Regular Article
Fractal analysis of shape factor for matrixfracture transfer function in fractured reservoirs
^{1}
Key Laboratory of Tectonics and Petroleum Resources, Ministry of Education, China University of Geosciences, Wuhan 430074, PR China
^{2}
Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, PR China
^{*} Corresponding authors: cughzhang@163.com, caijc@cug.edu.cn
Received:
17
February
2020
Accepted:
27
May
2020
As the core function of dualporosity model in fluids flow simulation of fractured reservoirs, matrixfracture transfer function is affected by several key parameters, such as shape factor. However, modeling the shape factor based on Euclidean geometry theory is hard to characterize the complexity of pore structures. Microscopic pore structures could be well characterized by fractal geometry theory. In this study, the separation variable method and Bessel function are applied to solve the singlephase fractal pressure diffusion equation, and then the obtained analytical solution is used to deduce onedimensional, twodimensional and threedimensional fractal shape factors. The proposed fractal shape factor can be used to explain the influence of microstructure of matrix on the fluid exchange rate between matrix and fracture, and is verified by numerical simulation. Results of sensitivity analysis indicate that shape factor decreases with tortuosity fractal dimension and characteristic length of matrix, increases with maximum pore diameter of matrix. Furthermore, the proposed fractal shape factor is effective in the condition that tortuosity fractal dimension of matrix is roughly between 1 and 1.25. This study shows that microscopic pore structures have an important effect on fluid transfer between matrix and fracture, which further improves the study on flow characteristics in fractured systems.
© L. Mei et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
q : Fluid transfer rate between matrix and fracture, [M/(TL^{3})]
N : Number of sets of fractures (1, 2 or 3)
L : Characteristic length of matrix, [L]
P : Matrix pressure, [M/(LT^{2})]
P _{0} : Initial pressure of matrix, [M/(LT^{2})]
: Average pressure of matrix, [M/(LT^{2})]
P _{ f } : Fracture pressure, [M/(LT^{2})]
K _{ m } : Matrix permeability, [L^{2}]
K _{ a } : Permeability constant, [L^{2}]
D _{ f } : Fractal dimension of matrix
D _{ T } : Tortuosity fractal dimension of matrix
λ _{max} : Maximum pore diameter of matrix, [L]
c : Total compressibility of reservoir, [(LT^{2})/M]
1 Introduction
Naturally fractured reservoirs account for a large proportion of the world’s resources and play an important role in the world’s energy structure. Unconventional resources have significantly transformed the landscape of oil and gas industry in these years. Hydraulic fracturing is a primary technology for economical and effective exploitation of unconventional reservoirs (Li et al., 2015; Wang et al., 2017; Behnoudfar et al., 2019; Vishkai and Gates, 2019), which aims at creating largescale fracture network to increase production. Therefore, fluids flow mechanism in fractured systems has aroused the strong interest of petroleum engineers, geothermal engineers and geologists, etc.
The dualporosity model is one of common simulation methods to study the flow properties in fractured reservoirs. The concept of dualporosity media in fractured reservoirs was firstly proposed by Barenblatt et al. (1960), and then was applied by Warren and Root (1963) to the field of petroleum engineering. As shown in Figure 1, fractured reservoirs are usually supposed to consist of a continuous fracture region that serves as the primary flow channels (high permeability and low porosity), and a discontinuous matrix region which acts as the primary hydrocarbon storage area (low permeability and high porosity). They hypothesized that singlephase of slightly compressible fluid transfer between matrix and fracture occurred under pseudosteady state conditions. The transfer function was originally formulated as (Warren and Root, 1963):(1)where q is the fluid transfer rate between matrix and fracture, [M/(TL^{3})]; σ is the shape factor, [1/L^{2}]; ρ is the fluid density, [M/L^{3}]; μ is the fluid viscosity, [M/(LT)]; K_{m} is the matrix permeability, [L^{2}]; is the average pressure of matrix, [M/(LT^{2})]; P_{f} is the fracture pressure, [M/(LT^{2})].
Fig. 1 Schematic of idealization of fractured reservoirs (Warren and Root, 1963). 
The shape factor model for cubic matrix can be written as (Warren and Root, 1963):(2)where N is the number of sets of fractures (1, 2 or 3), L is the characteristic length of matrix which is the distance between the fracture surface and the center of matrix, [L]. The shape factors are 12/L^{2}, 32/L^{2} and 60/L^{2} for onedimensional (1D), twodimensional (2D) and threedimensional (3D) spaces, respectively.
Furthermore, many scholars have carried out studies on fluids exchange between matrix and fracture based on Warren and Root (1963)’s model. Kazemi et al. (1976) directly extended singlephase transfer of Warren and Root (1963) to twophase transfer and firstly applied it to the numerical simulation of 3D isotropic dualporosity porous media. When the length of matrix is the same in each direction, the shape factors are 4/L^{2}, 8/L^{2} and 12/L^{2} for 1D, 2D and 3D spaces, respectively. By including the effects of gravity, capillary force and viscous force on transfer function, Thomas et al. (1983) established a 3D threephase model to simulate fluids flow in fractured reservoirs. By comparing singleporosity finegrid with the dualporosity models, they concluded that shape factor was 25/L^{2} when the two models have good agreement. Coats (1989) solved the diffusion equation of the slightly compressible fluids and derived a transfer function. They found that their shape factor was exactly twice the value of Kazemi et al. (1976). By comparing finegrid models with experiments, Ueda et al. (1989) believed that the shape factor was twice and three times that of Kazemi et al. (1976) in 1D and 2D spaces, respectively.
However, the application of these shape factors to explain the transient characteristics of matrix pressure in the initial stage of fluids transfer process usually produces large deviation. Therefore, some scholars have conducted many studies on the description of matrix pressure transient behavior. Zimmerman et al. (1993) developed a singlephase dualporosity model in fractured reservoirs and solved the nonlinear ordinary differential equation to obtain matrixfracture transfer function, which was more accurate than the linear Warren and Root (1963) equation to calculate the flux in the early and late stages. Lim and Aziz (1995) derived shape factor by solving matrix pressure diffusion equation, and used finegrid numerical simulation to validate the presented shape factor. They considered that both the geometric of matrix and pressure gradient of the whole system have a great impact on shape factor. Noetinger and Estebenet (2000) used continuous time random walks methods to simulate matrixfracture exchange term and obtained shape factor involved in dualporosity formulations of fluid flow through fractured reservoirs. Good agreement were found between their results and numerical simulation. Sarma and Aziz (2004) believed that fracture networks were nonorthogonal, and solved the nonorthogonal singlephase pressure diffusion equation to deduce 2D and 3D shape factors. Their results were verified by comparing the discrete fracture model with the dualporosity model. For different geometries and boundary conditions, Hassanzadeh and Pooladidarvish (2006) introduced the shape factor, and indicated that both the geometric of matrix and pressure change of fracture have an important effect on the shape factor. Ranjbar et al. (2011) solved the singlephase nonlinear diffusivity equation of different pressure depletion regimes to derive the shape factor, and investigated the impact of fracture pressure depletion regime on the shape factor. They further applied the finegrid numerical simulation technique to verify their shape factor. By solving the saturation diffusion equation in fractured reservoirs, SaboorianJooybari et al. (2015) developed the shape factor accounting for both capillary and gravity forces, which was verified using finegrid numerical simulation. For low permeability fractured reservoirs, He et al. (2017) proposed a new shape factor by solving singlephase pressure diffusion equation. Furthermore, their results shown that the new shape factor could predict production accurately. For tight fractured reservoirs with the existence of the boundary layer and the heterogeneous pressure distribution of matrix, Cao et al. (2019) presented the singlephase matrixfracture transfer function, which was verified by the numerical simulation and the experimental data.
Fractal geometry theory has been widely used to study the transport properties of porous media. The application of fractal geometry theory in the field of petroleum engineering has received extensive attentions. Both of fracture network (Chang and Yortsos, 1990; Acuna et al., 1995) and matrix (Katz and Thompson, 1985; Krohn, 1988; Adler, 1996) can be represented effectively and conveniently by fractal geometry. Zhou et al. (2000) studied the fractal characteristics of natural fractures, and developed a fractal transfer function model to analyze the effect of pore size distribution of matrix on the fluids transfer between matrix and fracture in fractured reservoirs. Based on the fractal characteristics of fractured reservoirs, FlamencoLopez and CamachoVelazquez (2001) established the fractal dualporosity equation, and obtained the analytical solution of pressure by Laplace transform method. The importance of transient and pseudosteady pressure characteristics of fractured reservoirs with fractal geometry in a singlewell was discussed. Kong et al. (2009) described the pore characteristics of porous and fractured media through fractal theory and presented the new flow velocity, porosity and permeability models in both porous and fractured media. In their study, the pressure diffusion equations of porous and fractured media were established and their transient pressure characteristics were analyzed. Yao et al. (2012) believed that fractures have fractal characteristics and established two dualporosity equations based on circular matrix and cylindrical matrix, respectively. They used Laplace method to get the pressure analytical solution, and analyzed the influences of matrix shape and fractal parameters on the transient pressure characteristics of fractured reservoirs. Fan and Ettehadtavakkol (2017) indicated that fractures have fractal characteristics and established induced fractures network distribution model, which was a function of density distribution and permeability/porosity distribution of induced fractures. Lian et al. (2018) used fractal geometry to characterize the different scale distribution of pores and fractures, and presented singlephase radial fluid flow equation by involving nonequilibrium adsorption to study the flow properties in shale gas reservoirs. Wang et al. (2018) combined the dual porosity model with the trilinear flow model, and used fractal geometry to characterize the heterogeneous distribution of fracture networks and the nonuniform of both matrix/fracture porosity and permeability. By applying Bessel Function and Laplace transform techniques, the flow of fluid through primary hydraulic fractures, stimulatedreservoirvolume and unstimulated reservoir matrices were integrated to derive analytical solutions. The accuracy of the analytical solutions in their study were verified by a synthetic finegrid numerical simulation example.
In the previous literature, the effect of matrix microstructure on the shape factor is not considered during reservoir development. Fractal geometry theory can relatively accurately characterize the effect of matrix microstructure on fluids transfer between matrix and fracture. Fractal geometry theory has been widely applied in flow characteristics of fractured reservoirs, but its applications in the transfer function and shape factor are still rare. In this paper, we consider the effect of matrix microstructure on fluids transfer between matrix and fracture by means of fractal geometry theory, from which a new shape factor is derived. Singleporosity finegrid and dualporosity models are then applied to validate the new shape factor. Finally, the influences of microstructure and characteristic length of matrix on the shape factor are accordingly discussed.
2 Mathematical model
It has been reported that pore surface and pore size distribution of porous media follow fractal features (Katz and Thompson, 1985; Yu and Cheng, 2002). The fractal geometry theory has great advantages in accurate representation of porous media with complex microstructures (Costa, 2006; Erol et al., 2017). Considering the fractal capillary bundle model and HagenPoiseuille equation of curved capillary, the fractal permeability and fractal porosity of matrix can be respectively expressed by (Kong et al., 2007):(3)(4)where ϕ_{m} is the matrix porosity; A_{0} is the unit area, [L^{2}]; x is the distance, [L]; D_{f} is the pore fractal dimension, 0 < D_{f} < 2 (in 2D space) and 0 < D_{f} < 3 (in 3D space); λ_{max} is the maximum pore diameter, [L]; D_{T} is the tortuosity fractal dimension, 1 < D_{T} < 2 (in 2D space) and 1 < D_{T} < 3 (in 3D space), D_{T} = 1 represents a straight capillary, D_{T} = 2 represents a plane is completely filled with tortuous line, and D_{T} = 3 represents a 3D space is completely filled with tortuous line; K_{a} is the permeability constant, [L^{2}]; ϕ_{a} is the porosity constant. It is noted that the symbols D_{f} and D_{T} are respectively used for pores and streamlines in matrix in this work.
Substituting equations (3) and (4) into the continuity equation, the matrix fractal pressure diffusion equation can be obtained by (Kong et al., 2007):(5)where ξ corresponds to 1D flow, 2D flow and 3D flow, which are 1, 2 and 3, respectively. When the capillary is straight, equation (5) can be simplified as the classical pressure diffusion equation.
Before establishing mathematical model, some assumptions need to be considered as follows:

The fractured reservoirs are seen as a dualporosity model (Fig. 2), where fluids in matrix flow through fractures instead of directly to the bottom of well, and matrix is isotropic.
Fig. 2 Schematic of dualporosity media.
The fluid is slightly compressible.
The rock is compressible.
The matrix is represented by the fractal tortuous capillary bundle model.
2.1 Transfer function and shape factor of 1D flow
Supposing matrix has been surrounded by two fractures, which have a spacing L (Fig. 3), fluids flow from matrix to fracture surface, and therefore the fractal pressure diffusion equation in matrix can be given by Kong et al. (2007):(6)where P is the matrix pressure, [M/(LT^{2})]; ϕ_{ao} is the value of ϕ_{a} under the reference pressure P_{0}; c is the total compressibility of reservoir, [(LT^{2})/M]; t is the transfer flow time, [T].
Fig. 3 Schematic of 1D matrixfracture system. 
Initial and boundary conditions are given by:(7)
Using separation variable method and Bessel function and combining equation (6) with equation (7) (see ^{Appendix A} for details):(8)where is a constant.
The partial derivative of equation (8) with respect to time t gives:(9)
Equation (9) can be also written as:(10)
The 1D transfer function can be expressed as:(11)
Substituting equation (10) into equation (11) yields:(12)
Compared to the transfer function (Eq. (1)), shape factor in 1D flow can be expressed as:(13)
From equation (13), the proposed shape factor considers the effect of flow channel bending and micro pore structure of matrix. Specifically, shape factor is clearly related to tortuosity fractal dimension, maximum pore diameter and characteristic length of matrix.
2.2 Transfer function and shape factor of 2D flow
In 2D model, the matrix is surrounded by four fractures, with two groups of parallel fractures. Fracture space is L. For the 2D matrixfracture system, the pressure diffusion behavior from matrix to fracture can be similar to a circular region (Fig. 4). The equivalent radius R of matrix is:(14)
Fig. 4 Schematic of 2D matrixfracture system. 
The 2D matrix fractal pressure diffusion equation can be given by (Kong et al., 2007):(15)where r is the distance, [L].
Initial and boundary conditions are given by:(16)
Using separation variable method, the analytical solution of equation (15) can be obtained in the form of Bessel function (see ^{Appendix B} for details):(17)where , J_{v} (x) is the Bessel function, v is the order of Bessel function, I_{1} is the positive zero point of Bessel function (the solution code of I_{1} is shown in ^{Appendix C}), it depends on the value of v, and .
Assuming that the area of matrix A_{av} is equal to πr^{2}. The average pressure of matrix is as follows:(18)
Substituting equation (17) into equation (18) yields:(19)
Both sides of equation (19) are subtracted by P_{0} and then divided by P_{f} − P_{0}. Thus, equation (19) can be expressed as:(20)
The partial derivative of equation (20) with respect to t can be expressed as:(21)
Equation (21) also can be written as:(22)
The 2D transfer function can be expressed as:(23)
Substituting equation (22) into equation (23) yields:(24)
Consequently, the 2D shape factor can be expressed as:(25)
Substituting equation (14) into equation (25) yields:(26)
The presented 2D shape factor also considers the effect of tortuosity fractal dimension, maximum pore diameter and characteristic length of matrix. Moreover, the value of the positive zero point of Bessel function is related to the value of tortuosity fractal dimension.
2.3 Transfer function and shape factor of 3D flow
In 3D model, the matrix is surrounded by six fractures, with three groups of parallel fractures, which have a spacing L. For the 3D matrixfracture system, the pressure diffusion behavior from matrix to fracture can be approximately regarded as a sphere (Fig. 5). The equivalent radius R of matrix is:(27)
Fig. 5 Schematic of 3D matrixfracture system. 
The 3D matrix fractal pressure diffusion equation can be given by (Kong et al., 2007):(28)
Initial and boundary conditions are given by:(29)
Using separation variable method, the analytical solution of equation (28) can be obtained in the form of Bessel function (see ^{Appendix B} for details):(30)where , I_{2} is the positive zero point of Bessel function (the solution code of I_{2} is shown in ^{Appendix C}), it depends on the value of v, and .
Assuming that the volume of matrix V_{av} is equal to , the average pressure of matrix is as follows:(31)
Substituting equation (30) into equation (31) yields:(32)
Both sides of equation (32) are subtracted by P_{0} and then divided by P_{f} − P_{0}. Thus, equation (32) can be expressed as:(33)
The partial derivative of equation (33) with respect to t can be expressed as:(34)
Equation (34) can be also written as:(35)
The 3D transfer function can be expressed as:(36)
Substituting equation (35) into equation (36) yields:(37)
Consequently, the 3D shape factor can be expressed as:(38)
Substituting equation (27) into equation (38) gives:(39)
Compared with those shape factors proposed by predecessors, the fractal shape factor accounts for the influences of two additional parameters: tortuosity fractal dimension and maximum pore diameter of matrix. The effect of microstructure on shape factor directly affects fluids transfer in matrixfracture system. When the flow channels in matrix are approximately straight, tortuosity fractal dimension of matrix D_{T} = 1, our shape factor model has the similar value with Lim and Aziz (1995) model. When the flow channels in matrix are curved, according to the experimental observation, Yin et al. (2017) presented that tortuosity fractal dimension of sandstone, coal and shale media is generally ranging from 1.1 to 2.3. Nelson (2009) reported that pore diameters are usually greater than 2 μm in conventional reservoirs, ranging from 2 to 0.03 μm and from 0.1 to 0.005 μm in tightgas sandstone and shale, respectively. Within the reasonable value range D_{T} = 1.1 and λ_{max} = 0.15 mm to calculate the positive zero point of Bessel function, we can get I_{1} = 2.474 and I_{2} = 3.013, and shape factors are 11.35/L^{2.1}, 20.38/L^{2.1} and 24.77/L^{2.1} for 1D, 2D and 3D spaces, respectively. The shape factors predicted by Warren and Root (1963), Kazemi et al. (1976), Lim and Aziz (1995) and our model are presented in Table 1.
The shape factors defined by various models.
3 Model verification
Here, the commercial reservoir simulator ECLIPSE is applied instead of the laboratory experiment to verify the proposed fractal shape factor. Three types of singleporosity finegrid model are developed using ECLIPSE, and offer reference curve for the dualporosity models with different shape factors. In the singleporosity finegrid models, matrix is discretized that grids near fractures are finer than others, and fracture spacing is 10 cm (Fig. 6). The dualporosity models only have one cubic grid with a size of 50 cm, which is mainly controlled by shape factor in transfer function, and fracture spacing is 50 cm. The basic parameters used in these models are shown in Table 2.
Fig. 6 Schematic of singleporosity finegrid models for (a) 1D, (b) 2D, and (c) 3D flow. 
Basic parameters used in numerical simulation.
Curves of cumulative production over time under 1D, 2D and 3D flows are displayed in Figures 7a–7c, respectively. The results of singleporosity finegrid models are taken as the reference curves. Results show that the fractal shape factors have relatively accurate matching degree with the reference curves. Clearly, the curvature of the fluid channels can’t be ignored in the real rocks. Both tortuosity fractal dimension and maximum pore diameter have an important role in the fluid flow in porous media (Yun et al., 2009; Ye et al., 2019; Huang et al., 2020). Thus, there will be a great error in the fluid exchange term if neglecting the effects of tortuosity fractal dimension and maximum pore diameter of matrix on the shape factor. The proposed fractal shape factor can be used to explain the impacts of tortuosity fractal dimension and maximum pore diameter on the fluid exchange rate, and give a better prediction of production rates and recoveries.
Fig. 7 Cumulative production curves of dualporosity model with different shape factors and singleporosity finegrid model, (a) 1D flow; (b) 2D flow; (c) 3D flow. 
4 Results and discussion
The fractal shape factor model presented in this paper (Eqs. (13), (26) and (39)) is in term of tortuosity fractal dimension, maximum pore diameter and characteristic length of matrix. The impacts of parameters on shape factor are studied in this part to understand the sensitivity to shape factor. It is generally believed that the permeability of lowpermeable porous media is less than 10 × 10^{−3} μm^{2}. With the decrease of permeability, the characteristics of nonlinear flow become obvious. Since the derivation of the presented fractal shape factor model follows Darcy’s law, so according to the experimental data (Yin et al., 2017), tortuosity fractal dimension of matrix is roughly less than 1.5. Six values of tortuosity fractal dimension of matrix are selected as 1.0, 1.1, 1.2, 1.3, 1.4 and 1.5, respectively, which reflect the different curvature of flow channels in porous media. A range of maximum pore diameters of matrix that varies from 0.1 μm to 2 mm, covering conventional rocks, tight sandstone and shale. The fractal shape factor model contains the positive zero point of Bessel function, which is related to tortuosity fractal dimension. So, we first solve the positive zero point of the Bessel function at different values of tortuosity fractal dimension, as shown in Table 3.
The positive zero point of Bessel function at different values of tortuosity fractal dimension.
The values of the fractal shape factor decreases as tortuosity fractal dimension increases (Fig. 8). Since the larger tortuosity fractal dimension is, the more curved the flow channel is, which will lead to greater flow resistance, thus reducing the fluid exchange mass between matrix and fracture. However, when tortuosity fractal dimension is larger than 1.25, the shape factor of 2D flow is greater than that of 3D flow. Since the number of fractures around matrix increases, a larger flow crosssectional area is formed, resulting in a greater fluid exchange rate between matrix and fracture. The shape factor of 2D flow should be less than that of 3D flow. Figure 8 indicates the presented fractal shape factor model is suitable for tortuosity fractal dimension less than 1.25.
Fig. 8 The impact of tortuosity fractal dimension of matrix on the fractal shape factor. 
The fractal shape factor growths with the increaseof maximum pore diameter of matrix (Fig. 9). The presented fractal shape factor model is applicable in all range of maximum pore diameter of matrix, indicating that the presented fractal shape factor model is suitable for dualporosity media at different scales and the presented fractal shape factor varies at different scales.
Fig. 9 The impact of maximum pore diameter of matrix on the fractal shape factor. 
The fractal shape factor decreases with the increase of characteristic length of matrix (Fig. 10). The increase of characteristic length of matrix means that the fluid flow time is longer, reducing the rate of fluids transfer between matrix and fracture.
Fig. 10 The impact of characteristic length of matrix on the fractal shape factor. 
In order to further understand the fluids transfer mechanism between matrix and fracture, different shape factor models have been presented in the literature. Warren and Root (1963) combined an integral material balance to derive shape factors are 12/L^{2}, 32/L^{2} and 60/L^{2} for 1D, 2D and 3D, respectively. Moreover, Lim and Aziz (1995) used approximating analytical solution of matrix pressure diffusion equation to derived shape factors are π^{2}/L^{2}, 18.17/L^{2} (or 2π^{2}/L^{2}) and 25.67/L^{2} (or 3π^{2}/L^{2}) for 1D, 2D and 3D, respectively. He et al. (2017) considered the influence of tortuosity (τ) and introduced the modified shape factors are π^{2}/(τL^{2}), 18/(τL^{2}) and 26/(τL^{2}) for 1D, 2D and 3D, respectively. When these shape factors are applied in the transfer function to investigate the fluids transfer between matrix and fracture, the impact of matrix microscopic pore structure is usually neglected, which is deviated from the actual situation. The presented fractal shape factor considers the impacts of tortuosity fractal dimension of matrix, maximum pore diameter and characteristic length. Compared with the previous models, the proposed fractal shape factor model can more accurately characterize the impact of microscopic pore structure on the rate of fluids transfer between matrix and fracture in fractured reservoirs.
5 Conclusion
We established new transfer function and shape factor to investigate the fluids transfer between matrix and fracture in fractured reservoirs. The analytical expression for fractal shape factors in 1D, 2D and 3D spaces are derived via the fractal pressure diffusion equation. The presented fractal shape factor model is authenticated by comparing dualporosity model with singleporosity finegrid model. Furtherly, sensitivity analysis of model parameters is carried out to understand their effects on fractal shape factor. Some conclusions can be obtained:
For different porous media, the fractal shape factor can provide corresponding values. Applying the fractal shape factor to transfer function will help us better understand the impact of microscopic pore structure on fluid exchange law between matrix and fracture.
There is a positive correlation between the fractal shape factor and maximum pore diameter of matrix. The fractal shape factor decreases with the growth of tortuosity fractal dimension and characteristic length .
When tortuosity fractal dimension is roughly less than 1.25, the fractal shape factor has practical application value in different scales.
Acknowledgments
This work is supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDA14010302), the National Natural Science Foundation of China (No. 51904279), and the Fundamental Research Funds for the Central Universities (China University of Geosciences, Wuhan) (No. CUGGC04).
Appendix A
Derivation of 1D fractal pressure diffusion equation
Supposing matrix has been surrounded by two fractures (Fig. 3), equation (A.1) for 1D matrix fractal pressure diffusion equation (Kong et al., 2007):(A.1)
Let W = P − P_{f}, and initial and boundary conditions can be rewritten as:(A.2)
Substituting equation (A.2) into equation (A.1), it yields:(A.3)
Using separation variable method, it must have:(A.4)where .
Substituting equation (A.4) into equation (A.3), we can obtain:(A.5)
According to the solution method of Bessel function, the solution of equation (A.5) can be obtained:(A.6)
Substituting equation (A.6) into equation (A.4), it yields:(A.7)
Combining equations (A.2) with (A.7), we can obtain:(A.8a)(A.8b)(A.8c)
Consequently, the solutions of equations (A.8a) and (A.8b) can be expressed as:(A.9)
Substituting equation (A.9) into equation (A.8c), equation (A.8c) can be rewritten as:(A.10)
Both sides of equation (A.10) are multiplied by and integrated with respect to x. Then, the value of B_{m} can be obtained as:(A.11)
Substituting equations (A.9) and (A.11) into equation (A.7), thus, equation (A.7) becomes:(A.12)
Equation (A.12) can be converted to:(A.13)
Assuming that the area of matrix A_{av} is equal to hx, h is the height which is constant, and x is distance. The average pressure of the matrix is as follows:(A.14)
Substituting equation (A.13) into equation (A.14), it yields:(A.15)where .
Both sides of equation (A.15) are subtracted by P_{0} and then divided by P_{f} − P_{0}. Thus, equation (A.15) can be expressed as:(A.16)
Appendix B
Derivation of 2D and 3D fractal pressure diffusion equation
The matrix fractal pressure diffusion equation is given by (Kong et al., 2007):(B.1)where ξ corresponds to 2D flow and 3D flow, which are 2 and 3, respectively.
Let G = P − P_{f}, and initial and boundary conditions can be rewritten as:(B.2)
Substituting equation (B.2) into equation (B.1), it yields:(B.3)
Using separation variable method, we can see that:(B.4)
Substituting equation (B.4) into equation (B.3), we can obtain:(B.5)
According to the solution method of Bessel function, the solution of equation (B.5) can be obtained:(B.6)where , v corresponds to 2D flow and 3D flow, which are and , respectively.
Substituting equation (B.6) into equation (B.4), it yields:(B.7)
Combining equations (B.2) and (B.7), we can obtain:(B.8a)(B.8b)
Equation (B.8b) can be converted to:(B.9)
In equation (B.9), when r → 0, we have , but P_{0} − P_{f} << ∞, so we can obtain:(B.10)
Then, equation (B.9) can be reduced to:(B.11)
Substituting equation (B.10) into equation (B.8a) yields:(B.12)
In equation (B.12), due to and , so we can obtain:(B.13)where J_{v} (x) is the Bessel function, v is the order of Bessel function, I_{m} is the positive zero point of Bessel function (the solution code of I_{m} is shown in ^{Appendix C}), it depends on the value of v, I_{m} corresponds to 2D flow and 3D flow, which are I_{1} and I_{2}, respectively.
Both sides of equation (B.12) are multiplied by and integrated with respect to x. Then, the value of c_{m} can be obtained as:(B.14)
Substituting equations (B.10), (B.13) and (B.14) into equation (B.7), thus, equation (B.7) becomes:(B.15)
Equation (B.15) can be converted to:(B.16)
Appendix C
Solution code of the positive zero of Bessel function
In this paper, we need to get the positive zero of Bessel function which is solved by MATLAB software. The specific code is as follows:
clear, clc;
format long
x=(0:0.2:1)′;
y_0=fzero(@(x)besselj(0.045,x),3);
References
 Acuna J.A., Ershaghi I., Yortsos Y.C. (1995) Practical application of fractal pressure transient analysis of naturally fractured reservoirs, SPE Form. Eval. 10, 3, 173–179. [CrossRef] [Google Scholar]
 Adler P.M. (1996) Transports in fractal porous media, J. Hydrol. 187, 195–213. [CrossRef] [Google Scholar]
 Barenblatt G.I., Zheltov I.P., Kochina I.N. (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, J. Appl. Math. Mech. 24, 5, 1286–1303. [CrossRef] [Google Scholar]
 Behnoudfar P., Asadi M.B., Gholilou A., Zendehboudi S. (2019) A new model to conduct hydraulic fracture design in coalbed methane reservoirs by incorporating stress variations, J. Pet. Sci. Eng. 174, 1208–1222. [Google Scholar]
 Cao R., Xu Z., Cheng L., Peng Y., Wang Y., Guo Z. (2019) Study of single phase mass transfer between matrix and fracture in tight oil reservoirs, Geofluids 2019, 1–11. [Google Scholar]
 Chang J., Yortsos Y.C. (1990) Pressuretransient analysis of fractal reservoirs, SPE Form. Eval. 5, 1, 31–38. [CrossRef] [Google Scholar]
 Coats K.H. (1989) Implicit compositional simulation of singleporosity and dualporosity reservoirs, in: SPE Symposium on Reservoir Simulation, Houston, Texas, 6–8 February, Society of Petroleum Engineers. [Google Scholar]
 Costa A. (2006) Permeabilityporosity relationship: A reexamination of the KozenyCarman equation based on a fractal porespace geometry assumption, Geophys. Res. Lett. 33, 2, L02318. [Google Scholar]
 Erol S., Fowler S.J., HarcouëtMenou V., Laenen B. (2017) An analytical model of porosity–permeability for porous and fractured media, Transp. Porous Media 120, 2, 327–358. [Google Scholar]
 Fan D., Ettehadtavakkol A. (2017) Semianalytical modeling of shale gas flow through fractal induced fracture networks with microseismic data, Fuel 193, 444–459. [CrossRef] [Google Scholar]
 FlamencoLopez F., CamachoVelazquez R. (2001) Fractal transient pressure behavior of naturally fractured reservoirs, in: SPE Annual Technical Conference and Exhibition, 30 September–3 October, New Orleans, Louisiana, Society of Petroleum Engineers. [Google Scholar]
 Hassanzadeh H., Pooladidarvish M. (2006) Effects of fracture boundary conditions on matrixfracture transfer shape factor, Transp. Porous Media 64, 1, 51–71. [Google Scholar]
 He Y., Chen X., Zhang Y., Yu W. (2017) Modeling interporosity flow functions and shape factors in lowpermeability naturally fractured reservoir, J. Pet. Sci. Eng. 156, 110–117. [Google Scholar]
 Huang T., Du P., Peng X., Wang P., Zou G. (2020) Pressure drop and fractal nonDarcy coefficient model for fluid flow through porous media, J. Pet. Sci. Eng. 184, 106579. [Google Scholar]
 Katz A.J., Thompson A.H. (1985) Fractal sandstone pores: Implications for conductivity and pore formation, Phys. Rev. Lett. 5412, 1325–1328. [CrossRef] [PubMed] [Google Scholar]
 Kazemi H., Merrill L.S., Porterfield K.L., Zeman P.R. (1976) Numerical simulation of wateroil flow in naturally fractured reservoirs, SPE J. 16, 6, 317–326. [Google Scholar]
 Kong X., Li D., Lu D. (2007) Basic formulas of fractal seepage and typecurves of fractal reservoirs, Journal of Xi’an Shiyou University (Natural Science Edition) 22, 2, 1–10. [Google Scholar]
 Kong X., Li D., Lu D. (2009) Transient pressure analysis in porous and fractured fractal reservoirs, Sci. China, Ser. E: Technol. Sci. 52, 9, 2700–2708. [CrossRef] [Google Scholar]
 Krohn C.E. (1988) Fractal measurements of sandstones, shales, and carbonates, J. Geophys. Res. 93, 3297–3305. [Google Scholar]
 Li Q., Xing H., Liu J., Liu X. (2015) A review on hydraulic fracturing of unconventional reservoir, Petroleum 1, 1, 8–15. [CrossRef] [Google Scholar]
 Lian P.Q., Duan T.Z., Xu R., Li L.L., Li M. (2018) Pressure behavior of shalegas flow in dual porous medium based on fractal theory, Interpretation 6, 4, SN1–SN10. [CrossRef] [Google Scholar]
 Lim K.T., Aziz K. (1995) Matrixfracture transfer shape factors for dualporosity simulators, J. Pet. Sci. Eng. 13, 3, 169–178. [Google Scholar]
 Nelson P.H. (2009) Porethroat sizes in sandstones, tight sandstones, and shales, AAPG Bull. 93, 3, 329–340. [CrossRef] [Google Scholar]
 Noetinger B., Estebenet T. (2000) Upscaling of double porosity fractured media using continuoustime random walks methods, Transp. Porous Media 39, 3, 315–337. [Google Scholar]
 Ranjbar E., Hassanzadeh H., Chen Z. (2011) Effect of fracture pressure depletion regimes on the dualporosity shape factor for flow of compressible fluids in fractured porous media, Adv. Water Resour. 34, 12, 1681–1693. [Google Scholar]
 SaboorianJooybari H., Ashoori S., Mowazi G. (2015) A new transient matrix/fracture shape factor for capillary and gravity imbibition in fractured reservoirs, Energy Sources Part A 37, 23, 2497–2506. [CrossRef] [Google Scholar]
 Sarma P., Aziz K. (2004) New transfer functions for simulation of naturally fractured reservoirs with dual porosity models, in: SPE Annual Technical Conference and Exhibition, Houston, Texas, 26–29 September, Society of Petroleum Engineers. [Google Scholar]
 Thomas L.K., Dixon T.N., Pierson R.G. (1983) Fractured reservoir simulation, SPE J. 23, 1, 42–54. [Google Scholar]
 Ueda Y., Murata S., Watanabe Y., Funatsu K. (1989) Investigation of the shape factor used in the dualporosity reservoir simulator, in: SPE AsiaPacific Conference, 13–15 September, Sydney, Society of Petroleum Engineers. [Google Scholar]
 Vishkai M., Gates I. (2019) On multistage hydraulic fracturing in tight gas reservoirs: Montney Formation, Alberta, Canada, J. Pet. Sci. Eng. 174, 1127–1141. [Google Scholar]
 Wang W., Yuan B., Su Y., Sheng G., Yao W., Gao H., Wang K. (2018) A composite dualporosity fractal model for channelfractured horizontal wells, Eng. Appl. Comput. Fluid Mech. 12, 1, 104–116. [Google Scholar]
 Wang W., Zheng D., Sheng G., Zhang Q., Su Y. (2017) A review of stimulated reservoir volume characterization for multiple fractured horizontal well in unconventional reservoirs, Adv. GeoEnergy Res. 1, 1, 54–63. [CrossRef] [Google Scholar]
 Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, SPE J. 3, 3, 245–255. [Google Scholar]
 Yao Y., Wu Y., Zhang R. (2012) The transient flow analysis of fluid in a fractal, doubleporosity reservoir, Transp. Porous Media 94, 1, 175–187. [Google Scholar]
 Ye W., Wang X., Cao C., Yu W. (2019) A fractal model for threshold pressure gradient of tight oil reservoirs, J. Pet. Sci. Eng. 179, 427–431. [Google Scholar]
 Yin S., Xie R., Ding W., Shan Y., Zhou W. (2017) Influences of fractal characteristics of reservoir rocks on permeability, Lithologic Reserv. 29, 4, 81–90. [Google Scholar]
 Yu B., Cheng P. (2002) A fractal permeability model for bidispersed porous media, Int. J. Heat Mass Transf. 45, 14, 2983–2993. [Google Scholar]
 Yun M., Yu B., Cai J. (2009) Analysis of seepage characters in fractal porous media, Int. J. Heat Mass Transf. 52, 13, 3272–3278. [Google Scholar]
 Zhou D., Ge J., Li Y., Cai Y. (2000) Establishment of interporosity flow function of complex fractured reservoirs, OGRT 7, 2, 30–32. [Google Scholar]
 Zimmerman R.W., Chen G., Hadgu T., Bodvarsson G.S. (1993) A numerical dualporosity model with semianalytical treatment of fracture/matrix flow, Water Resour. Res. 29, 7, 2127–2137. [Google Scholar]
All Tables
The positive zero point of Bessel function at different values of tortuosity fractal dimension.
All Figures
Fig. 1 Schematic of idealization of fractured reservoirs (Warren and Root, 1963). 

In the text 
Fig. 2 Schematic of dualporosity media. 

In the text 
Fig. 3 Schematic of 1D matrixfracture system. 

In the text 
Fig. 4 Schematic of 2D matrixfracture system. 

In the text 
Fig. 5 Schematic of 3D matrixfracture system. 

In the text 
Fig. 6 Schematic of singleporosity finegrid models for (a) 1D, (b) 2D, and (c) 3D flow. 

In the text 
Fig. 7 Cumulative production curves of dualporosity model with different shape factors and singleporosity finegrid model, (a) 1D flow; (b) 2D flow; (c) 3D flow. 

In the text 
Fig. 8 The impact of tortuosity fractal dimension of matrix on the fractal shape factor. 

In the text 
Fig. 9 The impact of maximum pore diameter of matrix on the fractal shape factor. 

In the text 
Fig. 10 The impact of characteristic length of matrix on the fractal shape factor. 

In the text 