Numerical Simulation of Fluid Flow through Fractal-Based Discrete Fractured Network

In recent years, multi-stage hydraulic fracturing technologies have greatly facilitated the development of unconventional oil and gas resources. However, a quantitative description of the “complexity” of the fracture network created by the hydraulic fracturing is confronted with many unsolved challenges. Given the multiple scales and heterogeneity of the fracture system, this study proposes a “bifurcated fractal” model to quantitatively describe the distribution of induced hydraulic fracture networks. The construction theory is employed to generate hierarchical fracture patterns as a scaled numerical model. With the implementation of discrete fractal-fracture network modeling (DFFN), fluid flow characteristics in bifurcated fractal fracture networks are characterized. The effects of bifurcated fracture length, bifurcated tendency, and number of bifurcation stages are examined. A field example of the fractured horizontal well is introduced to calibrate the accuracy of the flow model. The proposed model can provide a more realistic representation of complex fracture networks around a fractured horizontal well, and offer the way to quantify the “complexity” of the fracture network in shale reservoirs. The simulation results indicate that the geometry of the bifurcated fractal fracture network model has a significant impact on production performance in the tight reservoir, and enhancing connectivity of each bifurcate fracture is the key to improve the stimulation performance. In practice, this work provides a novel and efficient workflow for complex fracture characterization and production prediction in naturally-fractured reservoirs of multi-stage fractured horizontal wells.


Introduction
It has been commonly recognized that significant spatial heterogeneity, characterized by the multi-scale nature, extensively exists in shale reservoirs.The wetting properties and flow behavior of oil and gas change dramatically in different types of porous media.Cross-scale nano-micro porosity plays a dominant role in fluid storage, and micro-scale pore-throat morphology and connectivity strongly affected fluid transport phenomena.In addition, the fracture growth and the occurrence of natural fractures are extremely complex.Figure 1 shows multi-scale natural fractures, namely micrometer-scale (small scale), centimeter-scale (medium scale), and meter-scale (large-scale) fractures, co-exist in the shale reservoirs, of which the distribution pattern and occurrence both exhibit fractal features.Macro-scale favored flow areas are, accordingly, derived from the multi-scale connected/isolated natural fractures, localized accumulation of organic matter, and heterogeneous porous structure of shale reservoirs, which further affect the dynamic changes of pressures and fluid compositions in shale reservoirs.In addition, the complexity associated with reservoir, fracture, and fluid properties make various engineering techniques inaccurate to evaluate reservoir and fracturing properties of shales, such as sampling, well logging, and pressure-transient and rate-transient analysis [1,2].However, compared with individual fractures, characterizing the complexity of multiscale fracture networks still have many challenges.To accurately capture of the fluid flow through the complex fracture network created by stimulation in unconventional reservoirs, previous studies focus on directly improving methodologies characterizing individual fractures for fracture networks.To the best of our knowledge, the existing models describing such complex systems are divided into three categories, including the dual-porosity model, the wire-mesh model, and the unconventional fracture model (UFM).The dual porosity model was first introduced by Warren and Root [3] to characterize the behavior of naturally-fractured reservoirs and is now widely applied in the shale gas fracture modeling [4][5][6][7].Chaudhary [8] and Agboada et al. [9] simulate the fracture network with the help of logarithmic grid refinement, coupled with the dual porosity model, which is more flexible and able to offer a better description of pressure and saturation changes near the fracture.Hence, more accurate characterization of the flow pattern in the complex fracture network in the shale oil reservoir is obtained.By using the commercial simulator CMG, Odunuga [10] studied the self-similar reservoir fracture network fractal pattern and calculated the production capacity and drainage area of the fractured shale gas well with representative fractal patterns.The dual-porosity model is also used to shed light on the flow characteristics in shale gas reservoirs [11][12][13].In the continuum media dual medium-based model, the fracture property is averaged and then assigned to each grid node, and it is generally reflected in the concept of the matrix and fracture, in which it is difficult to accurately characterize every detail of the complex fracture system.Carlson [14] studied the superior advantages of the discrete fracture model in dealing with the explicit fracture description of the non-continuum medium, which can provide more information concerning the inter-porosity flow, given that the fracture spacing is relatively large.A simplified discrete fracture model suitable for general-purpose reservoir simulators is presented to handle both 2D and 3D matrix-fracture systems [15].Both the wire-mesh model and the unconventional fracture model originate from the discrete fracture model, which can explicitly describe the fracture network geometry based on fracture propagation.Xu et al. [16] established the wire-mesh model for hydraulic fracturing simulation, where there are two groups of planar fractures perpendicular to each other in the ellipse of the stimulated area.A new hydraulicfracture model is developed to simulate the complex fracture-network propagation in a formation with pre-existing natural fractures proposed by Weng et al. [17].The interaction, coupling, and deformation of hydraulic fractures and pre-existing natural fractures is taken into consideration in their model.
However, with the multi-scale features in particular shale reservoir, the guiding concept has changed from "bi-wing planar fractures" to "complex fracture networks".To account for the clastic mass after fracturing, Guo et al. [18] conducted several fracture propagation experiments to evaluate the crushability, Figure 2. Al-Obaidy et al. [19] proposed "branched fractal" models that describe the pressure depletion in gas condensate shale reservoirs, and their results show that the fractal models allow for using more realistic shale permeability and improve the match for the observed trend.To capture part of the complexity of the fracture network, Mhiri et al. [20] introduced a novel approach to model the hydraulic fractures using "random-walk" stochastic method.The geometry of a fracture network is related to the fractal controlling parameters, and the multi-level bifurcated geometry of a fractal fracture better represents fractures of different orders and branches.Fuyong et al. [21,22] investigate fluid flow in fractal porous media and propose simulation method for the permeability of porous media based on the multiple fractal models.④ and ⑤ grow only single fractures perpendicular to the wellbore, but ①, ②, ③, and ⑥ split into two branches after propagating a few centimeters, modified from Guo et al. [18] The objective of this work is to apply fractal theory to model a multi-scale fracture system and provide a quantitative characterization of the induced hydraulic fracture networks.The fractalfracture network model will be built to reflect the distribution pattern of the fractal-branching network based on the discrete fracture method.This approach is expected to give a more realistic representation of fracture geometry and multi-scale features observed in the unconventional reservoir.The findings of this research provide theoretical foundations for the characterization of the complex fracture network created by the multi-stage fracturing, while offering guidance for the production prediction.

Physical Model
The fracturing treatment aims at creating the stimulated reservoir volume (SRV) to form the hydraulic fracture-natural fracture system [23][24][25][26].However, without effective fracturing-forced propagation and openness, the natural fractures are very difficult to contribute to the production.Figure 3 shows a schematic plot of the fracture representative element volume (FREV), in which connected natural fractures and induced fractures are the main flow paths in the SRV, where permeability is much higher than the matrix.Given the process and relevant characteristics of the multi-stage hydraulic fracturing in unconventional reservoirs, the fractal fracture pattern is generated by construction theory [27] in Figure 4. Assumptions are made as follows: (1) Only the hydraulic fractures and natural fractures (including each level of induced fractures) connected with each other and isolated natural fractures are not considered.(2) Fractal fractures are generated hierarchically, and they are all perpendicular to the wellbore and grow by means of "bifurcation".(3) Fracture properties stipulate in each level of the fractal fractures by scaling factors.According to the fractal self-similar construction theory [27], the fractal dimension and the scaling factor of fractures can be determined.In terms of the symmetrically-distributed fracture network on each side of the main fracture, the branching fracture in the same level has the same length, while the branching fractures grow further following certain ratio factors.As hydraulic fractures propagate further, the fracture branch continuously increases, and the extended branching fractures become smaller than the previous cracks.Meanwhile, fracture properties, such as fracture width and fracture conductivity, with respect to each level, constantly decline.According to the deterministic self-similar fractal network, the branching fracture network is assumed to consist of N levels of branching fractures from the primary fracture to the secondary fracture, with a maximum fracture number of m.It is defined that the SRV of the secondary branching fractures and fracture network gradually decrease with the fracture propagation and, therefore, we introduced the concept of the stimulated fracture volume (SFV).The SFV of the fractal branching network is defined as being equal to the volume occupied by the fractures, which can be calculated by summarizing the effective stimulated volume of each level's fractures.
where,   =     ℎ represents the volume of the bifurcated fracture.

Mathematical Model of Discrete Fractal-Fracture Network
Following the mathematical model, we assumed a compressible single-phase fluid (i.e., oil) flow in the matrix and fracture systems, following Darcy's law.The fractured horizontal well is located in the center area of the rectangular reservoir; multiple fractures with finite conductivity vertically cross the horizontal wellbore; the stimulated reservoir volume (SRV), distributed on the two sides of the fracture, is composed of the matrix/fracture medium.Dimensionless variables are defined, as follows: =   /,   =   / (3) where  is the horizontal length,   and   represent the distance to the boundary in x and y directions, respectively,   ,   ,   , and   are the dimensionless distances,   is the matrix permeability,   represents the fracture network permeability inside SRV,   represents the hydraulic fracture permeability, and   ,   , and   represent the dimensionless permeability of the matrix, fracture network, and hydraulic fracture, respectively.
=         +     (9) where   is the dimensionless time,   is the dimensionless pressure,  is the flow capacity coefficient, and  represents the storability ratio.
The motion equation for fluid flow in the matrix is: where,   is the flow velocity tensor in the matrix, 10 −3 m/s;   is the matrix pressure, MPa; and  is the fluid viscosity, mPa•s.To account for the compressibility of rock and fluid, it is defined that  =  0  −  (  − ) ,  =  0  −  (  − ) .

Flow in the Matrix
where   is the volume flow rate in unit volume, and t represents the time, s.
To account for the formation matrix properties in both SRV and unstimulated reservoir volume, the surface source in three dimensions turns into a superposition of the line source in twodimensional space.Then, we can obtain the dimensionless mathematical flow equation: The boundary condition and initial condition can be written as follows: (  ,   ,   ;   = 0) where   is the natural fracture compressibility, MPa −1 ;   represents the sink/source term in the complex fracture system;   is the fracture network system pressure;   is the porosity of natural fractures; ( −  ′ ) represents the function of Delta, when  =  ′ , the function equals 1 and, otherwise, equals zero.Initial condition: (, , ;  = 0) =   (, , ;  = 0) =   (16) Boundary condition: (, , ; ) =   (, , ; ) represents the hydraulic fracture pressure, MPa.

Flow in the Hydraulic Fracture
The hydraulic fracture is regarded as a single-continuum medium and the diffusivity equation can be obtained as follows: where   is the hydraulic fracture compressibility, MPa −1 ;   is the hydraulic fracture permeability, and qF represents the sink/source flow term in the hydraulic fracture.Initial condition: Boundary condition: The dimensionless governing equation inside the hydraulic fractures is expressed as: Based on triangular grids and the finite element method, a planar 2-D unsteady flow mathematical model is established in the article, and we solve the model by using the finite element method.Please refer to Appendix A for numerical solution details.

Results and Discussion
Based on the discrete fractal fracture model (DFFM), the effects of the fracture pattern and conductivity on the well performance are analyzed.The numerical model built for the representative element volume of the fracture is applied to estimate the performance following every iteration stage, which assumes that the homogeneous reservoir and each staged fractal fracture are open, and the existing natural fractures are not modeled.Since conductivity declines further from fracture networks [29,30], the multi-scale conductivity is also considered in the simulation.

Fractal Fracture Network Pattern
With the fractal branching fracture propagation, hydraulic fractures encounter natural fractures and divert to form next-level branching fractures.The fracture bifurcation pattern is created via numerical programming, and the fracture geometry is illustrated in Figure 5.The connectivity for each branching fracture varies, with the initial hydraulic fracture conductivity of 200 D•cm.After five iterations, the value decreases to 2.4 D•cm, which only accounts for 12% of the main fracture.Figure 6 shows the fracture width and reservoir contact area, as well as the induced fracture numbers, which indicates that with increasing induced fracture numbers (the fracture network becomes more complex further away), even the fractal bifurcation width decreases evidently, but it leads to exponential growth of the reservoir contact area.The actual production data from the Changqing oil field are used for history matching so as to validate the proposed model, and the result is shown in Table 1.After generating the fractal fracture cases, all five models are found to give similar daily production in the previous two years.Induced fractal fracture geometry varies significantly, and the well performance remains very close (Figure 7a).However, the induced fracture model predicates much higher cumulative volumes in the next 15 years, and the production volume with eight induced fractures estimated from this model is 20% larger than that from a single fracture (Figure 7b), which indicates the great contribution of the induced fracture network to the well performance in the late production stage.Figure 8 shows the stimulated fracture volume (SFV) and cumulative oil production change with respect to the induced fracture numbers.The complexity of the fracture network rises with the stimulated fracture volume of branching fractures.Higher induced fracture numbers can generate more complicated networks, but they also bring strong interference between the connected induced fractures, which can lead to a production decrease in the end.

Multi-Scale Fractal Fracture Network Conductivity
Due to the features of brittle, naturally-fractured shale formations, massive shear failure occurs during hydraulic fracturing.When the rough rock surface encounters shear slippage, it cannot restore to original status and shear-induced permeability is maintained, which is known as the "selfpropped" fracture network.In order to investigate the effect of multi-scale conductivity, fracture conductivity is considered to decrease with the bifurcation, which means the natural decrease of fracture conductivity from the primary fracture to the secondary and tertiary branches (Figure 9).The fully-propped fracture network refers to the case where the proppants can be transported to the secondary fracture network far from the wellbore.The partially-propped case means that after the fracturing fluid is pumped into the formation, the proppants can only migrate near the primary fracture and massive "self-propped" fracture networks are generated at distant fracture systems.Meanwhile, the conductivity of the fracture network gradually declines as proppants break, migrate, or block the formation.In terms of the un-propped fracture network, the main fracture is effectively propped, while the secondary fracture and tertiary branches are only self-propped, which means no more proppant migration in the fracture networks.Based on the fracture geometry in Figure 5, we propose a model to investigate the effects of fracture propped patterns (Table 2) on cumulative oil production.The result shows that cumulative oil production of the fully-propped fracture network greatly exceeds that of the other two propped patterns.This indicates that the models with arbitrarily propped fracture geometry are superior to the single fracture.The partially-propped fractures' cumulative oil volume decreases significantly due to the decline of the fracture permeability and fracture width.Figure 10 also indicates that with the rapid decrease of complex fracture conductivity, the cumulative production increases slightly with growing complexity of the fracture network.Therefore, improving the effective conductivity of the fracture network is the precondition for hydraulic fracture stimulation.Regarding the un-propped fracture model, the cumulative production fluctuation according to the geometry, topological structure, and tortuosity of the fracture network can cause extra fluid flow resistance while the fracture network conductivity is relatively low.Without the propped network in the stimulated reservoir volume, the primary fracture is the main contributor to production.The decrease of the main fracture length is bound to have a great impact on fracture production, and also connected, induced, complex fractures interfering with each other can lead to lower performance.
We also investigate he cumulative production with different induced fracture deviation angles (Figure 11).It is illustrated that the reservoir contact area of the fracture network increases with the bifurcation angle (deviation angle).When fractures are featured by high conductivity, it is inappropriate to create a denser network in the stimulated reservoir volume, and only a few interweaving fractures can result in very good performance.The cumulative oil production of the partially-propped fracture case is fairly close to the others, and the effect of fracture morphology on yield is not great.When it comes to lower conductivity fractures, the effect of the main fracture substantially overweighs the secondary and tertiary fractures.

Conclusions
This paper investigates the fractal characteristics of the complex fracture network created by multi-stage fracturing in shale reservoirs.It has been found that discrete fractal-fracture network (DFFN) models allow for using more realistic fracture values to evaluate the well performance.(1) We establish the numerical simulation model for fractal fracture network characterization.The proposed model meshes with the unstructured grid, which is solved using the finite element method.
(2) We proposed a new method for fractal fracture network characterization according to the construction theory.(3) The correlation between multi-scale fractal fracture patterns, conductivity, and corresponding production has been quantitatively analyzed with the help of the mathematical formula.The proposed research may provide valuable insight into optimal hydraulic fracturing design and unconventional resource recovery maximization.For future extensions, a stochasticbased fractal fracture network combined with micro-seismic events may be coupled to quantify the complex network fractures; this will improve hydraulic fracture design and production performance.

Nomenclature
The approximate solution of differential fluid flow equations is derived using the Galerkin method of weighted residuals.
The tetrahedron characteristic matrix of the reservoir matrix is expressed as: .
In a similar way, the 2D hydraulic fracture characteristic matrix is: ).
We combine the above characteristic matrices of different regions to obtain the matrix of the whole reservoir, and simplify the equation sets according to the obtained matrix elements and the corresponding serial numbers.The equation sets are discrete by backward differentiation on the timescale.In this paper, the triangle mesh grid and adaptive mesh subdivision technology are adopted in the model and we demonstrate that the results converge to a solution and are independent of mesh size.Refinement of the mesh has no effect on the solution.

Figure 1 .
Figure 1.Multi-scale heterogeneity of reservoir system in shale reservoirs.

Figure 2 .
Figure 2. Core fracture network experiment result.A total of six cracks are generated after fracturing.④ and ⑤ grow only single fractures perpendicular to the wellbore, but ①, ②, ③, and ⑥ split into two branches after propagating a few centimeters, modified from Guo et al.[18]

Figure 3 .
Figure 3. Multi-stage fractured horizontal well with representative element volume (dark area).

Figure 4 .
Figure 4. Determined self-similar fractal fracture model by construction theory.(a) represent branching fracture propagation geometry, modified from [23]; (b) represent the fracture representative element volume in Figure 2; (c) is determined Y-shape self-similar fractal fractures [28].

Figure 6 .
Figure 6.Fractal fracture width and reservoir contact area change with respect to induced fracture numbers.

Figure 7 .
Figure 7. Fractal fracture network history matching combined with Changqing oil field data.(a) represent production rate with respect to different induced fracture scenario; (b) represent the cumulative production for all cases.

Figure 8 .
Figure 8. Stimulated fracture volume and cumulative oil production change with respect to induced fracture numbers.

Figure 9 .
Figure 9. Schematic plot for a propped discrete fracture network.

*
DR is declined ratio of fracture permeability, whose value equals the ratio of the secondary fracture permeability over the previous fracture permeability.

Figure 10 .
Figure 10.Cumulative production comparison among fractal bifurcate fracture with respect to different fracture network conductivity.

Figure 11 .
Figure 11.Cumulative production comparison among fractal bifurcate fracture geometries with respect to different fracture network conductivity.(a) reflect the cumulative production changed with induced fracture deviation angle; (b) represent the pressure distribution at the end of simulation.

Table 2 .
Different propped fracture network scenario and parameters.
, and   is the dimensionless fracture width.By applying interpolation function in the element, and the dimensionless pressure in element "s" can be expressed as: ∇ ,  ∇ , Ω           +       ∬  ,   , Ω     is the width of the primary fracture,  , equals [ 1 ,  2 ,  3 ], ∇ ,  =