Fractal analysis of shape factor for matrix-fracture transfer function in fractured reservoirs

As the core function of dual-porosity model in fluids flow simulation of fractured reservoirs, matrix-fracture 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 single-phase fractal pressure diffusion equation, and then the obtained analytical solution is used to deduce one-dimensional, two-dimensional and three-dimensional 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.

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 large-scale 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 dual-porosity model is one of common simulation methods to study the flow properties in fractured reservoirs.
The concept of dual-porosity 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 single-phase of slightly compressible fluid transfer between matrix and fracture occurred under pseudo-steady state conditions. The transfer function was originally formulated as (Warren and Root, 1963): where q is the fluid transfer rate between matrix and fracture, [ The shape factor model for cubic matrix can be written as (Warren and Root, 1963): 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 one-dimensional (1D), two-dimensional (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 single-phase transfer of Warren and Root (1963) to two-phase transfer and firstly applied it to the numerical simulation of 3D isotropic dual-porosity 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 three-phase model to simulate fluids flow in fractured reservoirs. By comparing single-porosity fine-grid with the dual-porosity 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 fine-grid 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 single-phase dual-porosity model in fractured reservoirs and solved the non-linear ordinary differential equation to obtain matrix-fracture 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 fine-grid 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 matrix-fracture exchange (a) Actual reservoir (b) Idealized model reservoir Fig. 1. Schematic of idealization of fractured reservoirs (Warren and Root, 1963). term and obtained shape factor involved in dual-porosity 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 non-orthogonal, and solved the nonorthogonal single-phase 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 single-phase 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, Saboorian-Jooybari et al. (2015) developed the shape factor accounting for both capillary and gravity forces, which was verified using fine-grid numerical simulation. For low permeability fractured reservoirs, He et al. (2017) proposed a new shape factor by solving single-phase 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 single-phase matrix-fracture 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, Flamenco-Lopez and Camacho-Velazquez (2001) established the fractal dual-porosity equation, and obtained the analytical solution of pressure by Laplace transform method. The importance of transient and pseudo-steady pressure characteristics of fractured reservoirs with fractal geometry in a single-well 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 dual-porosity 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 single-phase 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 tri-linear flow model, and used fractal geometry to characterize the heterogeneous distribution of fracture networks and the non-uniform of both matrix/fracture porosity and permeability. By applying Bessel Function and Laplace transform techniques, the flow of fluid through primary hydraulic fractures, stimulated-reservoir-volume and unstimulated reservoir matrices were integrated to derive analytical solutions. The accuracy of the analytical solutions in their study were verified by a synthetic fine-grid 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. Single-porosity fine-grid and dual-porosity 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.

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 Hagen-Poiseuille equation of curved capillary, the fractal permeability and fractal porosity of matrix can be respectively expressed by (Kong et al., 2007): 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); k 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): where n 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: 1. The fractured reservoirs are seen as a dual-porosity model (Fig. 2), where fluids in matrix flow through fractures instead of directly to the bottom of well, and matrix is isotropic. 2. The fluid is slightly compressible. 3. The rock is compressible.
4. The matrix is represented by the fractal tortuous capillary bundle model.

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): 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]. Initial and boundary conditions are given by: Using separation variable method and Bessel function and combining equation (6) with equation (7) (see Appendix A for details): where The partial derivative of equation (8) with respect to time t gives: oP ot  Equation (9) can be also written as: The 1D transfer function can be expressed as: Substituting equation (10) into equation (11) yields: Compared to the transfer function (Eq. (1)), shape factor in 1D flow can be expressed as: 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.

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 matrix-fracture 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: The 2D matrix fractal pressure diffusion equation can be given by (Kong et al., 2007): where r is the distance, [L]. Initial and boundary conditions are given by: Using separation variable method, the analytical solution of equation (15) can be obtained in the form of Bessel function (see Appendix B for details): where Þ 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 v ¼ D T À1 2D T . Assuming that the area of matrix A av is equal to pr 2 . The average pressure of matrix is as follows: Substituting equation (17) into equation (18) yields: 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: The partial derivative of equation (20) with respect to t can be expressed as: oP ot Equation (21) also can be written as: The 2D transfer function can be expressed as: Substituting equation (22) into equation (23) yields: Consequently, the 2D shape factor can be expressed as: Substituting equation (14) into equation (25) yields: 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.

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 matrix-fracture 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: The 3D matrix fractal pressure diffusion equation can be given by (Kong et al., 2007): Initial and boundary conditions are given by: Using separation variable method, the analytical solution of equation (28) can be obtained in the form of Bessel function (see Appendix B for details): where -¼ 1 , 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 v ¼ D T À2 j j 2D T . Assuming that the volume of matrix V av is equal to 4 3 pr 3 , the average pressure of matrix is as follows: Substituting equation (30) into equation (31) yields: 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: The partial derivative of equation (33) with respect to t can be expressed as: Fig. 4. Schematic of 2D matrix-fracture system.
oP ot Equation (34) can be also written as: The 3D transfer function can be expressed as: Substituting equation (35) into equation (36) yields: Consequently, the 3D shape factor can be expressed as: Substituting equation (27) into equation (38) gives: 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 matrix-fracture 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 lm in conventional reservoirs, ranging from 2 to 0.03 lm and from 0.1 to 0.005 lm in tight-gas sandstone and shale, respectively. Within the reasonable value range D T = 1.1 and k 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.

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 single-porosity fine-grid model are developed using ECLIPSE, and offer reference curve for the dual-porosity models with different shape factors. In the single-porosity fine-grid models, matrix is discretized that grids near fractures are finer than others, and fracture spacing is 10 cm (Fig. 6). The dual-porosity models only have one cubic grid with a size of 50 cm, which   Lim and Aziz (1995) p 2 =L 2 18:17=L 2 25:67=L 2 This paper (D T = 1.1, k max = 0.15 mm) 11:35=L 2:1 20:38=L 2:1 24:77=L 2:1 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.
Curves of cumulative production over time under 1D, 2D and 3D flows are displayed in Figures 7a-7c, respectively. The results of single-porosity fine-grid 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  Single-porosity (1D flow) Grid dimensions 12 Â 1 Â 1 Grid spacing x = 10, 1.5, 1.5, 3, 6, 13, 13, 6, 3, 1.5, 1.5, 10 cm y = 50 cm z = 50 cm Single-porosity (2D flow) Grid dimensions 12 Â 12 Â 1 Grid spacing x = 10, 1.5, 1.5, 3, 6, 13, 13, 6, 3, 1.5, 1.5, 10 cm y = same as x z = 50 cm Single-porosity (3D flow) Grid dimensions 12 Â 12 Â 12 Grid spacing x = 10, 1.5, 1.5, 3, 6, 13, 13, 6, 3, 1.5, 1.5, 10 cm y = z = same as x Dual-porosity Grid dimensions 1 Â 1 Â 1 Grid spacing 50 cm Â 50 cm Â 50 cm Common data for all models Matrix porosity 0.6 Matrix permeability 1 mD Fracture porosity 0.07 Fracture permeability 600 mD Matrix initial pressure 30 MPa Fracture initial pressure 20 MPa 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.

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 low-permeable porous media is less than 10 Â 10 À3 lm 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 lm 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 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 cross-sectional 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.
The fractal shape factor growths with the increase of 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.
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. 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 p 2 =L 2 , 18:17=L 2 (or 2p 2 =L 2 ) and 25:67=L 2 (or 3p 2 =L 2 ) for 1D, 2D and 3D, respectively. He et al. (2017) considered the influence of tortuosity (s) and introduced the modified shape factors are p 2 = sL 2 À Á , 18= sL 2 À Á and 26= sL 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, maximum pore diameter and characteristic length of matrix. 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.

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 single-porosity fine-grid model. Furtherly, sensitivity analysis of model parameters is carried out to understand their effects on fractal shape factor. Some conclusions can be obtained: 1. 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. 2. 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 of matrix. 3. When tortuosity fractal dimension is roughly less than 1.25, the fractal shape factor has practical application value in different scales.