Distinctive Flow Regions in Crossform Fracture Model in Shale Gas Reservoir Using Numerical Density Derivative Part 3

The numerical density derivative approach is used to measure fluid densities around the wellbore and to generate pressure equivalent for each phase using simplified pressure-density correlation. While statistical derivative method determines fluid phase permeabilities and also average effective permeability for a given reservoir system with new empirical model. Both methods were only tested in conventional oil and gas reservoir system.


Introduction
Pressure transient analysis (PTA) is the industry's most recognized and acceptable method for assessing well deliverability, skin, nearwellbore permeability and characterize reservoir heterogeneities for hydraulically fractured wells. For these wells, several flowing regions may occur in or around due to the 3D nature of formations flow geometry for which the radial flow symmetry do not often exists. These flow regions are difficult to define by basis of pressure transient data because of near wellbore and formation factors, such as penetration ratio (the ratio of the fractures height to the formation height), inclination angle from the vertical direction, the spacing between fractures, heterogeneities such as vertical and horizontal permeability's and anisotropy [1]. These parameters influence the well sand face pressure and derivative response.
Since early seventies, PTA industry's experts and researchers have developed several models considering different well, reservoir and boundaries conditions to describe the pressure transient behavior with or without hydraulic fractures in vertical or horizontal wells. These models were developed based on the source solution and Green's function to solve unsteady-state flow problem in the reservoir which was presented by Gringarten and Ramey [2]. Also, the Newman product method and source function have been used for solving transient flow problem interpreting pressure behaviors.
Cinco-Ley et al. [3] developed the concept of finite flow capacity and applied semi analytical approach to illustrate the importance of finite fracture when the F CD <300 which is similar to long fractures and low capacity fractures. Their idea facilitated the evaluation of massive hydraulic fracturing programs, although with limitation applicable to systems with small, constant compressibility. Also, their type curve is presently the reference for data analysis from a constant-rate flow test or a pressure-build-up test, depicting vertical hydraulic fracture model in an infinite-acting reservoir. In their study, they introduced a relationship between dimensionless time and pressure behavior which depends on time, and dimensionless fracture conductivity, F CD : Bennett et al. [4] established finite conductivity type-curves to distinguish the linear and bilinear flow region with a straight line for multi layered reservoirs using analytical solution for cases of constant pressure and rate. They concluded that this approach is applicable only if the productive interval is within the fracture and that the fracture conductivity is dependent on depth.
Zerzar et al. [5] integrated the boundary element method and Laplace transformation to publish a comprehensive solution for multiple vertical fractures in horizontal wells. In this study, seven flow regions were identified which include bilinear, first linear, elliptical, radial, pseudo-radial, second linear and pseudo-steady state.
Several studies on modeling the fractured flow patterns in hydraulically fractured wells have been done by researchers over the last four decades with well documented results in various engineering and mathematical research journals. In all of these researches, four flow regions have been reliably in dentified to occur in the reservoir with hydraulically fractured well. These flow regions are highlighted below: Firstly, linear flow: which is due to flow from fluid expansion along the fracture parallel to the wellbore. Occasionally, the wellbore storage effect could mask its response. Its occurrence depends on the length of the test and the fracture conductivities. This flow regime is recognized as a 1/2 slope in the log-log pressure derivative diagnostic plot and is used to determine fracture half-length, channel or reservoir width if vertical permeability is known [6]. Secondly, bilinear flow: a combination of a combination of two simultaneous linear flows in perpendicular directions. This only occurs for finite-conductivity fracture where linear flow exists both in the fracture and to the fracture plane. This flow regime is recognized as a 1/4 slope in the log-log pressure derivative diagnostic plot and is used to determine the fracture conductivity [7].
Thirdly, pseudo-radial flow with fractures of all conductivities and in most cases as late time features. It does occur after sufficient long flowing period.
Lastly, trilinear flow model has been developed over the last decade to account for flow from dual fracture features. Notably, research on this topic has centered on modeling trilinear flow in finite conductivity fractures in tight gas formation. Further researches are been carried out to ascertain the flow source and sink. This paper focuses on developing mathematical model for different (linear, bilinear and trilinear) flow regions in a non-conventional reservoir completed with vertical well and within a cross form fractures.

Cross form fracture flowing region model derivation
For cross form fractures symmetrical to a reference point (the wellbore) at an arbitrary angle to the horizontal axis X as shown in Figure 1, the following assumptions are made to formulate the general solution: • The reservoir system is dual porosity in nature • The system is naturally fracture reservoir consisting of natural fractures.
• Flow are via the two fractures into the wellbore • Pressure at the wellbore is the sum of combining pressure units along the fractures Dimensionless pressure: Dimensionless length coordinate (L), assuming isotropic properties • The matrix part act as a uniformly source of flow distribution into the fracture • Viscosity is constant and slightly compressible fluid • Reservoir is on a rectangular shape with producing well located at the centre.
• Transient interporosity flow model is adopted for the matrix and fracture transient flow.
• The fractures are modelled as homogeneous slab porous median with primary fracture=Lf 1 and secondary fracture=Lf 2 .
• The fractures have width W f and fully penetrate the entire pay zone.
• Flows into the fractures are along the fractures and no flow via the fractures tips.
• The fractures centreline is resolved along X axis and Y axis by virtue of the angle of inclination.
• Figure 2 represent a 2D pictorial view of the fracture pathway along the wellbore, then the expected flow behavior around and away from the wellbore. Also a 3D view of a centered secondary fracture intersecting two primary fractures is shown in Figure 3.
In terms of flowing pressure behavior; the following which is not in sequential order are possible: • Single/Dual linear flow into the well.
• Possible Pseudo radial flow close to the well but could be masked by wellbore storage • Bi/Tri-linear flow away from the well • Pseudo radial flow at the end of the fractures The benefit of the cross form fracture model is that it creates a high conductive path close and some distance away from the well bore which allows large surface area exposure of low permeability formation resulting to more flow into the wellbore. In this case, large volume of fluid per unit of time is produced into the wellbore resulting in an increased production rate without drilling another well.
The following dimensionless parameters are defined for the formulation:   The dimensionless variable rescaling the anisotropy system to an equivalent isotropic system is given as; The dimensionless fracture conductivity is defined as: Interporosity flow parameter: as defined by Warren and Roof, 1963 Dimensionless storativity: For fracture one (Primary fracture), the diffusivity equation is given as: And for fracture two (Secondary fracture) If the fractures are of same length, then First the diffusivity equation for the matrix is given as: For the matrix and fracture interflowing period, the diffusivity equation is similar for both fractures, therefore the following boundary condition BC are applicable Initial condition: Initial boundary: Resolving the matrix diffusivity, equation (14) into dimensionless form using equation (2, (3) and (6), we have: Where interporosity flow parameter 2 12 m The fracture solution for the crossform fracture model is formulated as follows: The primary fracture The diffusivity equation as stated in equation (12) is given as: Resolving the equation in dimensionless form using equation (2), (3) and (6) cos R l f = along x axis, we have: For secondary fracture resolve in y axis from equation (13) Converting the equation dimensionless form using equation (2), (3) and (6) Substituting for Assuming the flux is uniform along the fracture and the pressure at the wellbore is a summation of pressure along the fractures segment, hence fD x axis fD P = y axis Combining equation (23) and (26) Resolving in x axis See full detail at the Appendix.
The general solution combining the matrix and fractures differential equations with different BCs is given as: This is general solution for the crossform fractures model connecting at the wellbore. This equation can be inverted to obtain time dependant solution using Laplace inversion such as Stehfest's inversion algorithm. Figure 4 shows the flow path and expected flow regions for the crossform fracture in a vertical well. The flow geometry phase system for crossform is summarized in Figure 5 below.

Case (I) -----bilinear flow:
If ω=0 (33) The general solution is given as; This equation is due to bilinear flow period.
However considering the assumptions, the flow regime is limited by; This is the general solution for; 1 2 5.57 This equation is due to the linear flow period with assumption limiting to t.

Asymptons
Recall that for matrix slab ( ) Substitute equation 44 and the asympton Converting to time dependent function using Laplace inverse

Example
Numerical simulation is performed with a synthetic rock and fluid data. First, the reservoir is discretized into blocks using the Bennett et al. [4] empirical guideline for designing the fractures pathway. A single layer reservoir is discretized into 10,000 blocks with distribution as x:y:z=100:100:1 with the crossform fracture and well modeled in such a way that there are no boundary effects (Table 1). Bennett's [4] recommendation for designing x and y grids in a fracture model as shown in Table 2 is used to determine the block dimensions and fractures aperture. To adhere to the physics of fluid flow, the following assumption are considered: • Isothermal condition with no diffusion and dispersion process.
• No chemical reactions (thermodynamically equilibrium) • Single phase Tables 1 and 3 present a summary of the well and reservoir synthetic data used for the build-up and drawdown simulated scenarios with additional information given below. It is require generating pressuredensity equivalent and derivative for each fluid phase, comparing their diagnostic signatures and also identifying the cross form fracture flowing regions.

Assumption
• Shale Gas reservoir, completed with one well.
• Model with Bennett's [4] recommendation on grid sizes close to the well to account for density changes.
• Only flowing condition is simulated.
Gas, Oil, Water densities and WBHP around the local grid refinement (wellbore) are output using the simulator keywords. The following scenarios were evaluated: As stated earlier, the fracture blocks are modeled based on Bennett's [4] recommendation in which the fracture's block dimension in x and y axis are increasing to maximum value for each well grid block with different block dimensions used to model the x direction. In this case, the fracture half length in x direction is the distance between the tip of the fracture (minimum x dimension) and the well. Also the adjacent grids blocks dimensions are increased until the maximum value, then all the next grid blocks are assigned the same dimension. However, the well is completed in the minimum dimension of the grid block in y direction [7].
As shown in Figure 8, there is an uneven distribution of grid block in both x and y directions. The grid block are modeled in such a way that the dimensions of the adjacent grid blocks increases to the maximum and then have constant dimension. This indicates that there is uneven distribution of grid blocks in the reservoir. However, since the reservoir is symmetrically related to the fracture position and the well, a quarter of the reservoir has been observed. For the crossform fracture model, the fractures are reoriented in the x and y direction.
The vertical well is located in the center of the square reservoir and the grid block has both minimum x and y dimensions. Since the fracture is modeled along x and y direction, the finite conductivity fracture is parallel to the x and y axis and totally intersects the well symmetrically.
The best value for modeling the fracture width is 2ft, which is the dimension of the smallest grid block with the well. The equivalent fracture porosity is calculated from equation 50 below since the fracture porosity of 35% corresponds to the fracture width of 0.5ft.

Parameters
Design Value       Tables 3 and 4. Calculating fracture permeabilities, the following input data were used are listed below: The two approaches for analyzing the transient flow behavior for this study include ( For the constant rate production, the gas rate is fixed at 500Mscf/ day while for the constant pressure production, the BHFP is assumed to be 3500psi. The bottomhole flowing rate and fluid densities are determined from the result of the numerical simulation for constant pressure and constant rate production respectively. For constant bottomhole pressure solution, the multiphase fluid distribution is triggered at the wellbore in order to capture the density changes for each phase and calculate fluid pressures densities equivalent at bottomhole flowing conditions. Then the derivative dimensionless rate is calculated from the Pressure Equivalent of Density Weighted Average (PDENDWA→PDENA) by Biu and Zheng [9]. The pressure distribution at each simulated time step along the crossform fracture path is shown in Figure 5.       Region 5 is trilinear flow, accounting for flow from dual fracture features. Its response is similar to the bilinear flow response but a slight deviation from the bilinear flow curve depicts its presence. This Trilinear flow regime is believed to be caused by transient drainage of low permeability matrix blocks into adjoining fractures and parallel flow into the fractures depending on the length of the fractures and permeability distribution. Finally is the flow boundary dominated transient response [12].  Figure 11. Likewise, the result from sensitivities on gas production rate as shown in Figure 12 which indicates the log-log plot of gas rate versus time depicts three flow regions such as Linear-Bilinear-Linear for F CD <5 and two flow regions Linear-Linear for F CD >25 [13].
For constant flowing rate condition, the fracture conductivity, F CD for each dimensionless pressure in Bennett type curves is dependent on the fracture porosity. The fracture porosity is calculated from equation (50), then the real fracture porosity and its equivalent is recalculated due to the change in fracture width.
First, the density change for each phase at the wellbore is obtained and the fluid pressures densities equivalent is calculated at bottom hole flowing conditions. Then, the inverse derivative dimensionless pressure is calculated from the Pressure Equivalent of Density Weighted Average (PDENDWA→PDENA) by Biu and Zheng [9].
In this case, two fracture flow regions: Linear-Bilinear are depicted in this model as shown in Figure 13. Region 1 is the linear flow which is due to transient flow only in the fractures and Region 4 is bilinear flow and occurs when the matrix drainage begins simultaneously with the transient flow in the fractures The fractures apertures were increased from the equivalent fracture width w e of 2ft to 16ft (incremental of 2ft). In this case, it was discovered that for constant rate solution, the smaller the fracture aperture, the lower the number of fracture regions to be seen. At w e >2ft, three flow regions Linear-Bilinear-Linear were depicted as seen in gas production rate sensitivity with constant bottomhole pressure ( Figure 14) [14,15].

Conclusion
The following inferences were drawn from constant rate and pressure solutions reviewed; • A new mathematical model was developed for interpreting pressures behaviour of a vertical well with crossform fracture in shale gas reservoir using numerical density approach.     • A major advantage is that the method simplified the complex fracture-matrix flow model applying ordinary laplace transform model OLTM to formulate linear, bilinear and tri-linear flow mathematical equations.
• With the limit of the fracture model developed, five flowing region is identifiable • It also indicates that pressure responses and distinctive flow regions are influenced by mostly fracture's dimensions, conductivities and reservoir's boundaries.
• It has been demonstrated that for constant rate solution, the smaller the fracture aperture, the lower the number of fracture regions to be seen.
• The first flowing region which is the linear flow region which is the flow along the vertical plane parallel into the wellbore and the second as the Bilinear or Trilinear flow region which accounts for flow along the vertical plane parallel to the wellbore, then into the fracture after the pressure pulse reaches the upper and lower impermeable boundaries depending on the ratio of primary and secondary cross form fracture lengths and conductivities (Figures 15 and 16).

Detail Mathematical Model Derivation
The following dimensionless parameters are defined for the formulation (Figure 16).
Dimensionless Pressure: The dimensionless fracture conductivity is defined as: Interporosity flow parameter: as defined by Warren and Root, 1963 Dimensionless storativity: For fracture one (Primary fracture) the diffusivity equation is given as: For fracture two {Secondary fracture} If the fractures are of same length, then First the diffusivity equation for the matrix is given as: For the matrix and fracture interflowing period, the diffusivity equation is same for both fractures, therefore the following BC are applicable Initial condition: The characteristic equation for this differential equation is and its roots are: The general solution to this differential equation is given as: The constant A and B are obtained by the derivatives of equation xvii is given by The fracture solution for the crossform fracture model is formulated as follows: The primary fracture The diffusivity equation as stated in equation (xi) is given as: Substitute for l f1 = 4R 2 cos 2 θ along x axis, we have: For secondary fracture resolve in y axis from equation xxii Substituting for l f2 = 4R 2 sin 2 θ along y axis Assuming the flux is uniform along the fracture and the pressure at the wellbore is a summation of pressure along the fractures segment, hence fD P x axis fD P = y axis Combining equation (xxii) and (xxiii) Resolving in x axis (Spivey and Lee, 1999) provide a solution for multiple arbitrarily oriented Infinite fracture system with K anisotropy in an infinite slab reservoir.
The dimensionless variable rescaling the anisotropy system to an equivalent isotropic system is given as; 2 2 cos sin Since the flow region modelling is within the fracture tips, therefore from equation vi: 2 cos sin If the cross-sectional area of the wall face is far away the fracture face, The dimensionless Laplace domain is given as r.  K f = fracture permeability Therefore for constant rate, the inner BC for fracture face is given as: and the outer BC for no flow through the fracture hp in Laplace form is given as: Substitute the initial BC into equation (xxv) 2 2 ,0 3