Rate and pressure behavior considering the fractal characteristics of structurally disordered fractured reservoirs

The main objective of this paper is to understand the impact of the fractal characteristics of fractured reservoirs on their pressure behavior, flow rate decline, and productivity index. The paper proposes a new methodology for developing several analytical models for describing the wellbore pressure distribution and the flow rate decline trend. The proposed models consider including the fractal characteristics such as the mass fractal dimension, conductivity index of anomalous diffusion flow mechanism, fractal-network parameters, fractional-derivative order, and matrix/fracture-interaction index as well as dual-porosity media characteristics such as the storativity and interporosity flow coefficient in the analytical models of the pressure, rate, and productivity index. The study has found that: (1) Some of the fractal characteristics have a significant impact on reservoir performance, while others may not have a significant impact. (2) Fractal reservoirs exhibit better performance than the standard geometry reservoirs of single and dual-porosity media.


Introduction
Conventional and unconventional fractured reservoirs have been the focus of the petroleum industry and have received a considerable attention over the past few decades as they contribute by more than 60% of the total world oil and gas supply. Naturally fractured reservoirs (NFR) such as most of the carbonate reservoirs represent the conventional type of fractured reservoirs where the two-scale porous media (fracture/matrix) exist, while triple porous media exist in the unconventional fractured reservoirs where hydraulic fractures propagate in the matrix that has already a network of embedded natural fractures. The two types are well known by the complex structures, the heterogeneity, and the irregular flow paths or the nonuniform flow distribution controlled by the distribution of fractures. Moreover, reservoir fluids may undergo different flow patterns such as Darcy, non-Darcy flow, and diffusion flow, especially in the unconventional reservoirs with either classic or anomalous diffusion mechanisms.
There are two approaches for describing the geometry of these reservoirs. The first is the Euclidean geometry wherein the fracture networks are assumed to be connected and equivalent to a homogenous porous media, while the porous media in the second approach is characterized by multiple property scales and the non-Euclidean geometry (Chang and Yortsos 1990). The main concepts in the first approach are the existence of the dual-porosity and permeability media, the matrix and the fracture networks are Euclidean objects within which the fractures are uniformly distributed in the matrix, and the matrix is not interconnected i.e.,the reservoir fluids flow throughout the perfectly connected fracture networks toward the hydraulic fractures or the wellbores. However, a lot of evidence from the outcrops, well logs, and production behavior observations have indicated that some of the reservoirs may not follow the above-mentioned concepts (Camacho-Velazquez et al. 2008). Accordingly, the second approach wherein the non-Euclidean geometry may be more appropriate for describing these reservoirs. In this geometry, the fractal theory that considers a nonuniform distribution of fractures and the presence of fractures at different scales is more applicable. This theory states that fractal characteristics of non-Euclidean geometry such as fracture intensity and conductivity as well as the interaction between the matrix, and these fractures may have a significant impact on reservoir performance. It is also defined, according to Flamenco-Lopez and Comacho-Velazquez (2001), as a method to describe the network of fractures in a rock that connects the heterogeneities in small scales to those in bog scales. While Emanuel et al. (1989) explained that the fractal theory represents the reservoir heterogeneity between wells as a random variation superimposed on the smooth interpolation of well-logging values. The concept of describing reservoir heterogeneity and characterizing the spatial correlation structures of the porosity and permeability using the fractal theory was adopted also by Perez and Chopra (1997) who confirmed that the well logs of vertical and horizontal wells may help in this task.
Not all the oil and gas reservoirs could have the non-Euclidean geometry where the fractal theory could be applicable and fractal characteristics may have a significant impact on reservoir performance. Most of the carbonate reservoirs, typically classified as naturally fractured reservoirs, are most suitably described by the fractal theory than the Warren and Root's dual-porosity model approach. The formation heterogeneity, irregular fluid flow paths, and the presence of fracture clusters in some regions of the porous media of these reservoirs, and the emptiness of other regions from these clusters as well as the spatial changes in the matrix bock sizes and the fractures dimensions are all indicate the spatial variability of the petrophysical properties in the carbonate reservoirs. Because of the complexity of the unconventional reservoir structures where three porous elements are always existed; the matrix, natural fractures, and hydraulic fractures as well as the induced fractures by the fracturing process, they are considered a good candidate for the applications of the fractal theory. The petrophysical properties in these reservoirs undergo a significant spatial change with the distance from the hydraulic fracture face to the point where the no-flow boundaries are present (Al-Rbeawi 2018). Not only the petrophysical properties that could have the spatial changes in the unconventional reservoirs, the size of the matrix block, natural fractures, and induced fractures may change from nanopores scale to micro and macroscale (Ozcan 2014;Ozcan et al. 2014). The matrix block and fracture size as well as the matrix permeability distributions in the unconventional reservoirs may follow the power-law model (Camacho Velazquez et al. 2008) or the exponential distribution (Fuentes-Cruz and Valko 2015), while Al-Rbeawi 2018 stated that the bivariate lognormal distribution could be most fit the distributions of the matrix block size and the matrix permeability in these reservoirs as it is shown in Fig. 1. The power-law behavior of fracture size distribution, matrix permeability, and the fractal characteristics of the porous media in the structurally complex reservoirs have also investigated by Laubach and Gale (2006), and Ortega et al. (2006). Both groups have emphasized that the distribution of attributes such as length, height, and aperture of fractures can be frequently expressed by the power-law model, and hence, the scaling analysis could be used to infer other fracture attributes such as fracture strike, number of fracture sets, and fracture intensity.
The spatial changes in the matrix and fracture sizes and the matrix permeability set an environment for the diffusive flow mechanisms rather than the viscous or slip flow (Javadpour et al. 2007). Wu et al. (2014) stated that the fluid flow in the unconventional reservoirs, characterized by the discontinuity and heterogeneity that could exist in different locations, may undergo different influential processes such as adsorption/desorption process, non-Darcy flow, rock/fluid interaction, rock deformation within nanopores or microfractures, and multiscaled heterogeneity. Accordingly, the complex structure reservoirs might be dominated by the anomalous diffusion flow mechanisms wherein the fractional diffusivity equation with a power of derivative that is not equal to one is used in the flow equations (Raghavan 2011;Chen and Raghavan 2015;Albinali and Ozkan 2016;Raghavan and Chen 2017).
Classic diffusivity equation (the order of fractional derivative equals to 1.0) is used typically to describe pressure distribution in the porous media where classic diffusion exists, while anomalous diffusion (Ozcan 2014;Ozcan et al. 2014) refers to the diffusion flow phenomenon in heterogeneous reservoirs where disordered and complex structure characterizes porous media of these reservoirs.
Due to the difficulties accompanied the process of resembling the fractal theory and the flow equations in the heterogeneous reservoirs, neither the pressure transient nor rate transient models have considered chronically the fractal dimensions. The first serious attempt for this type of resembling was conducted by Chang and Yortsos 1990 and a new model was presented. The proposed model describes the pressure distribution in naturally fractured reservoirs that have different scales with poor fracture connectivity, and disorderly spatial distribution. The general formulation of this model was derived for a single-phase flow in a system consisting of a fractal object (fracture network) embedded in a Euclidean object (matrix). At the same time, Beier (1990) used a set of field data taken from Grayburg/San Andres formations in Southeastern New Mexico to indicate that applying the conventional models of homogenous reservoirs may not match the pressure records, while the models that consider the heterogeneity of the reservoir with the fractal structure give more realistic results that could match the records. Four years before, Hewett (1986) used a Euclidean reservoir with a permeability/porosity spatial distribution that obeys the fractal statistics (fractional Brownian motion). Different models for the pressure distribution in the fractal reservoir were introduced by Beier (1994) for a vertical well with a vertical hydraulic fracture. He considered the permeability distribution in the porous media obeys the bimodal distribution as it has either zero or constant value at any location in the reservoir. The proposed model by Chang and Yortsos was used by Acuna et al. (1995) as an alternative approach for the classical interpretation in well test analysis of naturally fractured reservoirs where the double-porosity model is in doubt. It is also adopted by Abbaszadeh (1995) for Euclidean porous media with embedded fractal objects and permeability/porosity variation that obeys the power-law function of the radial distance from a well. The characterization of this type of reservoir with the fractal geometry, according to Flamenco-Lopez and Camacho-Velazquez (2003), may need the analysis of pressure data representing transient and pseudo-steady state flow, otherwise, the fractal parameters should be determined by the porosity well logs or another type of information resources. Raghavan and Chen (2013a, b), and Al-Rbeawi and Owayed (2019) presented several models for the reservoir performance considering the anomalous diffusion flow in the fractal reservoirs that are depleted by hydraulic fractures. These models include the temporal and spatial anomalous diffusion flow coefficients i.e., the diffusivity equation has two fractional pressure derivatives with an order that is not equal to (1.0), but equals to the above-mentioned anomalous diffusion flow coefficients.
Unlike pressure transient analysis, resembling fractal theory with rate transient models is very rarely conducted and presented in the literature even though the rate transient analysis might be more favorable for characterizing the unconventional reservoirs because of the low and extralow permeability that could cause some uncertainties in analyzing the pressure behavior of these reservoirs. Camacho-Velasquez et al. (2008) presented analytical solutions during transient and boundary-dominated flow regime that show the possibilities to characterize the naturally fractured reservoir with fractal network using the production decline data. These solutions are applicable for the radial flow in the finite and infinite reservoirs with and without matrix participation. Recently, Raghavan and Chen (2017) presented several models for the rate decline considering the powerlaw behavior under the subdiffusion flow mechanism in the fractured reservoirs. This paper proposes several models for the pressure distribution, flow rate decline behavior, and cumulative production trend for a rectangular reservoir depleted by multiple hydraulic fractures. These models consider the fractal characteristics of the fractured reservoirs where fractal object represented by the fractures networks of multiple property scales i.e., hydraulic fractures, induced fractures by hydraulic fracturing stimulation, and naturally induced fractures are embedded in the Euclidean object represented by the matrix. All fractures networks are treated as non-Euclidean objects. The reservoir of interest is assumed consisting of stimulated volume (SRV) represents the drainage area between hydraulic fractures, and unstimulated volume represents the drainage area away from the fracture tips. The novel point introduced in this paper is represented by the fact that each SRV and the hydraulic fracture assigned to it have the same fractal characteristics. Furthermore, the SRV porous media is considered having a square-shaped configuration with a length and width equal to the hydraulic fracture half-length, and therefore, it could be treated as a single wellbore with a radius equals to the hydraulic fracture half-length. Two expected outcomes of the above-mentioned two points are validated by the results of the proposed models in this paper. The first is developing a pseudo-radial flow regime for the reservoirs with big unstimulated reservoir volume, while the second is observing the boundary-dominated flow regimes after a very short production time for the reservoirs with stimulated reservoir volume only.

Mathematical formulation of fluid flow in fractal reservoirs
Considering a rectangular reservoir with a length 2x e and width 2y e is depleted by a number of hydraulic fractures (n) of width 2w f and fracture half-length x f as it is shown in Fig. 2. These fractures are assumed to fully penetrate the formation in the vertical direction h f = h , symmetrically propagate in the porous media, and have the same spacing between them as it is demonstrated in Fig. 3. The reservoir of interest has two drainage areas. The first is the area between hydraulic fractures and stimulated by these fractures. It is called stimulated reservoir volume and characterized by a complex structure consisting of the hydraulic fracture, induced fractures by the stimulation process, and naturally induced fractures. Accordingly, this volume is best described by the fractal theory, while the second is the area away from the hydraulic fracture tips toward the reservoir boundary. It is called unstimulated reservoir volume and consists of the natural fractures embedded in the matrix; therefore, it is considered also fractal porous media as it is characterized by dual-porosity media or naturally fractured porous media.
The following assumptions are utilized for the mathematical formulation proposed in this study: 1. The fractal object is defined by each hydraulic fracture and its SRV. The fractal objects are similar. 2. Single-phase flow either oil or gas. 3. The SRV has a square a square-shaped area i.e., y = x f . 4. The no-flow boundary between hydraulic fractures is located at the mid-point of the spacing between these fractures. 5. Hydraulic fracture width and conductivity do not change along the fracture half-length.
While reservoir inner and outer boundary conditions are: where P Df is the dimensionless pressure drop inside hydraulic fractures. The definitions of the dimensionless parameters used in the mathematical models presented in this paper are given in "Appendix".
In this study, two fractal porous media are assumed in the drainage area. The first is the fractal porous media of the stimulated reservoir volume and determined by the hydraulic fracture length and the no-flow boundary as shown in Fig. 3. For this fractal object, the fractal characteristics are similar for the hydraulic fractures and the fractures embedded in the matrix of the stimulated reservoir volume. While the second is the fractal porous media of the unstimulated reservoir volume and determined by the area between fracture tips and reservoir boundaries where the fractal characterstics are controlled by the naturally induced fractures embedded in the matrix.
Furthermore, the stimulated reservoir volume assigned to each hydraulic is considered as a square-shaped drainage area, and reservoir fluid flows radially from the unstimulated reservoir volume to the stimulated volume as shown in Fig. 4.
The pressure distribution of the fractal system represented in Fig. 4 can be generated using the generalization of the  Warren and Root (1963) model and the concept of the fractal theory as follows (Chang and Yortsos 1990): Since three flow regimes might be developed during the early and intermediate production time before pseudo-steady state flow regime, the Euclidean dimension (D) in this study may have different values according to the flow geometry during these flow regimes. For example, the Euclidean dimension (D = 2.0) is best fit to the cylindrical flow geometry that would be used during pseudo-radial flow regime wherein reservoir fluid flows radially from the unstimulated porous media toward stimulated part in the reservoirs with large drainage areas. While (D = 1.0) is more applicable at early production time when the reservoir fluid flows linearly inside a rectangular cross-sectional area toward the wellbore during hydraulic fracture linear flow regime. Later, at intermediate production time, the bi-linear flow regime becomes the dominant flow wherein the Euclidean dimension is close to (D = 1.5) (Chang and Yortsos 1990).
The mass fractal dimension (d) in Eq. (1) represents the density of the fractures, naturally and hydraulically induced, in the matrix, while the parameter ( ) refers to the conductivity index of these fractures and responsible for the anomalous diffusion flow mechanisms in the porous media. Normal diffusion flow is typically associated with the fluid flow in homogenous reservoirs where the random Brownian motion of diffusive particles is obeying Gaussian probability density in which the variance is proportional to the first power of time. In this type of diffusion flow, the conductivity index is approaching ( = 0.0) , while for anomalous diffusion flow, it could be ( = 1.0) . O'Shaughnessy and Procaccia (1985) stated that the fracture conductivity index is related to the mass fractal dimension and suggested the following relationship: They also stated that the mass fractal dimension is correlated with the Euclidean dimension by the following model: The pressure distribution in the matrix can be written as follows (Chang and Yortsos 1990): In the above-mentioned models: while ( ) in Eq. (7) is the storativity. For the law permeability matrix block wherein the matrix does not participate in the rate response, the storativity could reach ( = 1.0) . The storativity is given by: where r in Eq. (12) is the matrix/fracture storativity ratio and given by: The generalized interporosity flow coefficient in Eq. (10) is developed using the same concept used in Warren and Root (1963) model. It is given by: where while the parameter (m) in Eq. (14) represents an equivalent property to the permeability and given by: The configuration of stimulated and unstimulated reservoir volumes while (A) in Eq. (14) is the drainage area assigned to each hydraulic fracture: The constant ( ) in Eq. (14) represents a parameter for the distance from the hydraulic fracture face to the no-flow boundary. It is an equivalent to the characteristic length used in Warren and Root (1963) model. In fractal reservoirs, this constant is a function of the fractal exponent D ′′ . For the Euclidean matrix where the fractal exponent D �� = 0.0 , the constant ( ) equals the spacing between hydraulic fractures i.e., (y) . It is calculated by: The hydraulic diffusivity coefficient used in Eq. (14) is given by:

Pressure behavior of fractal reservoirs-constant flow rate
The solution of Eq. (7) for the wellbore pressure drop in Laplace domain and dimensionless parameters considering closed rectangular reservoirs where i the matrix blocks participate in the production can be obtained for the early time production as.
where while the solution for late time production is given by: and (K) and (I) in Eqs. (22) and (26) are the modified Bessel functions.
The wellbore pressure drop and pressure derivative behaviors of a fractal reservoir of a square-shaped drainage area x eD = y eD = 2.0 with Euclidean dimension (D = 1.5) , and the conductivity index ( = 0.0) for different mass fractal dimensions (d) are shown in Fig. 5. Two different impacts for the mass fractal dimension can be recognized in Fig. 5. At early production time, the pressure drop increases with the increase of (d) , while it decreases with the increase of (d) at late production time. This would be explained by the fact that at early and intermediate production time, the production mostly comes from the stimulated reservoir volume where the fractal porous media are the common geometry. Therefore, the mass fractal dimension (d) may have more impact on the pressure drop than the Euclidean dimension (D) because of the complex structure in the stimulated reservoir volume characterized by the severe heterogeneity and the high disorder in the property distribution, while the unstimulated reservoir volume contributes most of the production during the late time wherein pseudo-steady state flow regime becomes the dominant flow. In this volume, the impact of the mass fractal dimension is not significant, and the flow is more controlled by the Euclidean dimension.
It is also inferred from Fig. 5 that the increase of the mass fractal dimension could cause a change in the dominant flow regime at early production time. The hydraulic fracture linear flow regime of a slope (1∕2) typically seen at very early production time might be terminated faster and bi-linear flow regime of a slope (1∕4) might be developed instead of it when the mass fractal dimension increases.
The impact of the Euclidean dimension (D) on the wellbore pressure drop and pressure derivative is represented in Fig. 6. It is worth saying that this impact is not significant and different than the behavior of the Euclidean system of dual-porosity reservoirs (D = 2.0) that could be obtained by applying the Warren and Root (1963) model. The higher the Euclidean dimension is; the lower wellbore pressure drop is observed. The impact of the Euclidean dimension is similar for all reservoir configurations and all other conditions rather than those used in preparing Fig. 6 wherein the mass fractal dimension (d = 2.0) and the conductivity index is ( = 0.0).
Neither the conductivity index ( ) nor the fractal-network parameter (m) have a significant impact on the pressure behavior of fractal reservoirs as it is depicted in Figs. 7 and 8. However, the drainage area size and the hydraulic fracture half-length have a significant impact on the wellbore pressure drop and the pressure derivative behavior as it is demonstrated by Figs. 9 and 10, respectively. Increasing the drainage area for the same hydraulic fracture characteristics (length and width) means increasing the spacing between these fractures and increasing the volume of the stimulated and unstimulated parts of the reservoirs. Accordingly, pseudo-radial flow might be developed for the big drainage areas wherein reservoir fluid flows radially from the unstimulated reservoir volume toward the stimulated volume. The pressure and pressure derivative curves in Fig. 9 confirm the fact that reaching pseudo-steady state flow may need a  for developing bi-linear flow regimes and hydraulic fracture linear flow regimes as it is indicated by the red, green, and light blue color curves in Fig. 10. While using a hydraulic fracture half-length equals the distance to the reservoir boundary x f = 1000.0 ft , only stimulated reservoir volume exists in the reservoir, causes developing the bi-linear flow regime for a very short time followed by long-time pseudosteady state flow regimes as it is shown in Fig. 10. Figures 9 and 10 indicate that the fractal and regular reservoirs with single and dual-porosity media may have the same response (in the shape but not in the magnitude) for the wellbore pressure and pressure derivative regardless of the drainage area size and hydraulic fracture  Fig. 10 Wellbore pressure and pressure derivative of different hydraulic fracture half-lengths half-length. This fact means that the fractal characteristics may not cause different impacts for the drainage area size and hydraulic fracture half-length on the wellbore pressure drop and pressure derivative than those known for regular reservoirs. The same conclusion is reached for the impact of the two main parameters in the Warren and Root (1963) model i.e., the storativity ( ) and interporosity flow coefficient ( ) on the pressure and pressure derivative behavior of the fractal and regular reservoirs. Figure 11 shows the impact of the storativity for fractal reservoirs with (D = 1.5) , (d = 2.0) , and ( = 0.5) . The wellbore pressure drop reaches the minimum value when the storativity equals to ( = 1.0) . Physically this would be explained by the fact the main contribution of the production is given by the fractures rather than the matrix. While decreasing the storativity indicates more participation for the matrix in the production, and thereby, more pressure drop is required.

Flow rate and cumulative production behaviors of fractal reservoirs-Constant wellbore pressure
Sandface flow rate behavior assuming constant wellbore pressure can be predicted using Duhamel's formula (van Everdingen 1949) in which: therefore, the solution for the early production time presented by Eq. (22) can be used for estimating the flow rate in dimensionless form and Laplace domain as follows: while the solution for the late production time is obtained from Eq. (26): Similarly, the cumulative production, in dimensionless form and Laplace domain, is calculated by Helmy and Wattenbarger (1998): therefore, for the early production time, it is: and for the late production time, it is: The observable trends for the flow rate and cumulative production, shown Fig. 12, show two different behaviors with time for different mass fractal dimensions. The first exhibits declining behavior in the flow rate and cumulative production during transient state flow, while increasing flow rate and cumulative production are seen during pseudosteady state flow when the reservoirs are characterized by high mass fractal dimension. This would be an indication of the negative impact of the fractal characteristics during transient early production time while these characteristics may have a positive impact at late production time.
Euclidean matrix with single fracture (D = 1.0) gives the minimum flow rate and minimum cumulative production compared to the matrix where naturally and hydraulically induced fractures are embedded and complex structures are formed in the porous media, especially the stimulated reservoir volume between hydraulic fractures (D > 1.0) . Unlike the changes caused by the mass fractal dimension, there is always a very slight increase in the flow rate and cumulative production with the increase of the Euclidean dimension. This increase, shown in Fig. 13, is seen more during the pseudo-steady state flow at late production time than the transit state flow at early and intermediate production time.
Similar to the regular reservoirs with Euclidean porous media, the flow rate and cumulative production increase significantly as the drainage area size increases during late (36) For the high values of the conductivity index to ( ≥ 0.5) , there is no impact for the fractal characteristics on the flow rate and cumulative production during transient state flow and only late production pseudo-steady state flow is affected by the fractal characteristics, as it is demonstrated by Fig. 15, wherein a very slight increase in the flow rate and cumulative production is seen compared to the small values of the conductivity index ( ≤ 0.5).

Productivity index of fractal reservoirs
Productivity index of constant sandface flow rate is calculated by: where P wD is the pressure obtained by converting Eqs. (22) or (26)  while productivity index of constant wellbore pressure is calculated by: where q D is dimensionless flow rate calculated by converting Eqs. (32) or (33) to the real-time flow rate and N PD is the cumulative production calculated from Eqs. (35) or (36). (39)

3
The behaviors of the productivity index for constant sandface flow rate J Dq and constant wellbore pressure J DP of fractal reservoirs with different mass fractal dimensions are shown in Fig. 16. Both indices exhibit a very slight decrease with the increase of the mass fractal dimension during the early transient state flow, while they demonstrate an increasing trend with the mass fractal dimension during the late pseudo-steady state flow. Figure 17 shows the two indices for fractal reservoirs with different Euclidean dimensions. Unlike the impact of the mass fractal dimension, there is only slightly increasing behavior for these two indices when the Euclidean dimension increases in the transient and The impact of the conductivity index ( ) on the two productivity indices of fractal reservoirs is presented in Fig. 18. The comparison of the results shown in Figs. 17 and 18 indicates that the productivity indices for the conductivity index

Flow regimes of fractal reservoirs
Four flow regimes might be observed during the entire production life of the fractal reservoirs. In a time-sequence manner, these are: hydraulic fracture, bi-linear, pseudo-radial, and pseudo-steady state flow regime. Pseudo-radial flow regime is seen only in the reservoirs characterized by very big drainage areas, while pseudo-steady state flow regime might not be reached in unconventional reservoirs or it could need for a very long production time because of the ultralow permeability. The hydraulic fracture linear flow regime represents reservoir fluid flow linearly inside hydraulic fractures toward the wellbore. This flow develops at very early production time and controlled by hydraulic fracture characteristics. There is no significant impact for the matrix characteristics on this flow regime. It is characterized by the slope of (1∕2) on the pressure derivative curves. Analytically, the mathematical models for this flow regime can be obtained by considering that the reservoir matrix does not contribute to the flow during early production time i.e., the production comes from the fluid inside hydraulic fractures only, therefore ( = 1.0) . It is also reasonable considering that the conductivity index of the fractures that is responsible for anomalous diffusion flow during the early production time is ( = 0.0) . Accordingly: At the early production time (t → 0.0) , therefore, (s → ∞) . Accordingly, the wellbore pressure drop assuming constant sandface flow rate can be obtained in dimensionless form in Laplace domain from Eq. (22): and in real-time: while the dimensionless pressure derivative is given by: The two mathematical models of wellbore pressure drop and derivative given by Eqs. (44) and (45) indicate that the hydraulic fracture linear flow regime is not affected by the characteristic of the fractal reservoirs. This fact is confirmed by the plot of pressure and pressure derivative shown in Fig. 19 for different fractal dimensions (d) . However, the impact of the fractal dimensions at early production time causes changing the slope of this flow regime from (1∕2) to (1∕4) i.e., hydraulic fracture linear flow regime could be changed to a bi-linear flow regime when the mass fractal dimension increases. Similarly, the flow rate and cumulative production during hydraulic fracture linear flow regime assuming constant wellbore pressure can be obtained from Eqs. (32) and (35), respectively in dimensionless form and in Laplace domain as follows: and in real-time: The results of the flow rate and cumulative production shown in Fig. 20 demonstrate the same indication obtained by the pressure behavior depicted in Fig. 19 wherein no impact for the fractal characteristics is seen on reservoir performance. Analytically, as it is represented in the two models of Eqs. (48) and (49) and graphically as it is shown in Fig. 19, it can be easily inferred that the slope of the flow rate decline curve is (−1∕2) and the slope of the cumulative production curve is (1∕2) during hydraulic fracture linear flow regime. Mathematically, the cumulative production and wellbore pressure drop have the same model in dimensionless form during hydraulic fracture linear flow regime as it is indicated from Eqs. (44) and (49); however, there is no physical reason for this similarity.
The second flow regime is the bi-linear flow regime that represents the flow linearly from the fractures embedded in the matrix toward the hydraulic fractures and the flow linearly from these hydraulic fractures toward the wellbore. This flow regime develops during the intermediate production time and is controlled by the characteristic of the stimulated reservoir volume (SRV) as well as the characteristics of hydraulic fractures. The fractal characteristics of the reservoirs may have a significant impact on this flow regime. It is characterized by a slope of (1∕4) on the pressure derivative curves as it is shown in Fig. 21. The bi-linear flow regime is very well developed for the reservoirs with high fractal Unlike the hydraulic fracture linear flow regime, the matrix may contribute to the production during the bi-linear flow regime, therefore ( ≠ 1.0) . Accordingly, the analytical model for the dimensionless wellbore pressure drop in realtime, assuming constant pressure, during this flow regime is obtained from Eq. (26) as follows: while the dimensionless flow rate and cumulative production, assuming constant wellbore pressure, is given, respectively, by: and where (50)  Figure 22 shows the flow rate and cumulative production behavior of fractal reservoirs during bi-linear flow regime. There is a significant impact for the fractal dimension on the flow rate and cumulative production where both of them (1 − ) increase with the increase of the fractal dimension (d) . Physically, increasing the fractal dimension means increasing the density of fractures in the porous media, and thereby, most production comes from the reservoir fluid in these fractures. The pseudo-radial flow regime might be observed after bi-linear flow regime and before pseudo-steady state flow. This flow regime is very well established in the reservoirs with big drainage area or when hydraulic fractures may not have long hydraulic fractures as it is shown in Fig. 23. In the two cases, there is enough space for reservoir fluids to flow radially in the horizontal plane toward the hydraulic fractures. The pseudo-radial flow regime represents the flow from unstimulated reservoir volume toward the volume. For Euclidian reservoirs or when the fractal dimension is close  Fig. 25 Flow rate and cumulative production during the pseudo-radial flow regime for different fractal dimensions respectively. Figure 24 depicts the pressure and pressure derivative of fractal reservoirs with different fractal dimensions during the pseudo-radial flow regime, while Fig. 25 shows the flow rate and cumulative production of pseudo-radial flow regimes. Unlike hydraulic fracture linear flow and bi-linear flow regime, the maximum flow rate and cumulative production are obtained for the fractal reservoir with fractal dimension (d ≥ 2.0).
The last seen flow regime is the pseudo-steady state flow regime. It develops when the production pulse reaches the boundary of the reservoir at a very late production time. It is characterized by a unit slope line for  Fig. 27 Flow rate and cumulative production behaviors during pseudo-steady state flow regime for different fractal dimensions the pressure and pressure derivative curves and both of them perfectly overleaped as it is shown in Fig.26. The analytical wellbore pressure drop during this flow regime can be calculated using Eq. (50), while the flow rate and cumulative production can be calculated using Eqs. (51) and (52), respectively. Figure 27 demonstrates the flow rate and cumulative production of fractal reservoirs during the pseudo-steady state flow regime considering different fractal dimensions. It is very well indicated that bigger fractal dimension leads to more decrease in wellbore pressure drop for constant flow rate production or maximum flow rate and maximum cumulative production for constant wellbore pressure. The analytical models for pseudo-steady state flow regime in fractal reservoirs, given by Eq. (50), can be approximated to:

Application
The proposed models in this study are examined twice. They are verified by the comparison with the previous works and validated by applying field data. In the verification process, the pressure drops for a rectangular reservoir using the solutions of the short and long production time given by Eqs. (22) and (26) are compared with the pressure drops obtained by Flamenco-Lopez and Camacho-Velazquez 2003. In the verification process, a hypothetical formation is considered with a rectangular drainage area of different configurations. This reservoir is assumed hydraulically depleted by multiple fractures. Several Euclidean dimension (D) , mass fractal dimension (d) , and the conductivity index ( ) . Because the reservoir configurations used by Flamenco-Lopez and Camacho-Velazquez (2003) are not given, type-curve matching process has been applied for several reservoir configurations and fractal parameters. An excellent matching, shown Fig. 28, has been observed for the following reservoir configuration and fractal parameters: The proposed models in this study for the wellbore pressure drop, flow rate, and cumulative production are validated also by a set of field data taken from Fuentes- Cruz and Valko (2015). The set of data belongs to a gas well producing in the Barnett shale, USA wherein the bottom hole flowing pressure and gas production rate are shown in Fig. 29. The formation and fluid data are given in Table 1.
Several attempts are conducted to match the field data shown in Fig. 29. The excellent matching is obtained for a fractal reservoir with fractal dimension (d ≅ 1.5) , while Euclidean dimension is used (D = 2.0) and the conductivity index, reflecting the anomalous diffusion flow, is ( = 0.2) . The reservoir in this case is assumed consisting of stimulated volume only i.e., there is no unstimulated reservoir volume x eD = 1.0, y eD = 1.0 . The wellbore pressure represented by the red dots is calculated assuming constant flow rate of x eD = 2.2 y eD = 1.85 D = 2.0 d = 1.7 = 0.08

Fig. 31
The flow regimes of the field example shown by the wellbore pressure derivative Q sc = 4.0MMScf ∕D , while the flow rate represented by the blue dots is calculated assuming constant bottom hole flowing pressure of P wf = 1000.0psi . Excellent matching is seen between the calculated flow rate by the proposed models in this study and the actual data, while good matching is seen between the calculated wellbore pressure and the actual recorded pressure. Figure 30 shows the comparison between the calculated wellbore pressure and flow rate using the proposed models in this study and the field data. Wellbore pressure derivative is calculated using the results of the calculated wellbore pressure by the proposed models for the case of constant flow rate as it is shown in Fig. 31. Three flow regimes can be easily recognized. The first is the hydraulic fracture linear flow regime. This flow regime develops at the early production time and finishes after 2 days. The second is the bi-linear flow regime that develops for almost 2 years followed by the pseudo-steady state flow regime. The reason for observing the pseudosteady state flow regime after short production time even though the extra-low permeability is the existence of stimulated reservoir volume only and the big flow rate used in the production Q sc = 4.0MMScf ∕D .

Conclusions
1. The fractal dimension (d) may have a significant impact on the wellbore pressure behavior, flow rate, cumulative production, and productivity index of fractal reservoirs. However, the conductivity index ( ) , and the fractalnetwork parameter (m) may not have significant impacts. 2. The wellbore pressure drop at the early production time increases as the fractal dimension increases, and however, at the late production time, the wellbore pressure drop shows increasing trend with the increasing of the fractal dimension. 3. Drainage area size and hydraulic fracture length have the same impact on reservoirs performance regardless of the geometry type of these reservoirs whether it is Euclidean or fractal porous media. 4. Bi-linear flow regime and pseudo-steady state flow regime are most affected by the fractal dimensions compared to the hydraulic fracture linear flow regime. 5. Fractal reservoirs exhibit a better stabilized pseudosteady state productivity index than the Euclidean reservoirs.
Funding There is no funding party for the research of interest.
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/.