Biu Victor T^{*} and Zheng ShiYi  
London South Bank University, United Kingdom  
Corresponding Author :  Biu Victor T London South Bank University, United Kingdom Tel: 009647505641986 Email: [email protected] 
Received: July 03, 2015; Accepted: July 16, 2015; Published: July 26, 2015  
Citation: Victor BT, ShiYi Z (2015) Distinctive Flow Regions in Crossform Fracture Model in Shale Gas Reservoir Using Numerical Density Derivative Part 3. J Pet Environ Biotechnol 6:245. doi:10.4172/21577463.1000245  
Copyright: © 2015 Victor BT, et al. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.  
Related article at Pubmed, Scholar Google 
Visit for more related articles at Journal of Petroleum & Environmental Biotechnology
The numerical density derivative approach is used to measure fluid densities around the wellbore and to generate pressure equivalent for each phase using simplified pressuredensity 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. This study introduces a new mathematical model for interpreting pressures behavior of a vertical well with cross form fracture in shale gas reservoir using numerical density approach. In this case, the imposed fractures can be longitudinal and transverse but symmetrical to a reference point (the wellbore). The major advantage is that it simplified the complex fracturematrix flow equation by applying ordinary laplace transform model OLTM to formulate linear, bilinear and trilinear flow model. The model is tested for constant pressure and constant rate conditions with the generated average fluid phase pressuredensities equivalent displaying the distinctive fractures flow fingerprint. It also indicates that the dimensionless rate or pressure derivative response and distinctive flow regions are influenced by mostly fracture’s conductivities, dimensions and reservoir’s boundaries. A new flow region have been added with the first as 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 is the 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. In this paper, it has been demonstrated that for constant rate solution, the smaller the fracture aperture, the reduction in the number of flow regions to be seen.
Keywords 
Pressure transient analysis; Drawdown test flow regions; Linear; Bilinear; Tri linear derivatives 
Nomenclature 
L_{fD}=Dimensionless fracture length 
Cs=Wellbore storage constant 
wf=Fracture width ft 
C_{A}=Area compressibility 1/psi 
Ct=Total compressibility 1/psi 
Δp=Change in pressure psia 
A=Drainage Area acres 
A_{mf}=Fracture crosssectional area to flow ft2 
Bgi=Formation volume factor at initial reservoir pressure, rcf/scf 
ct=Liquid total compressibility,1/psi 
dz=Well position in reservoir, dimensionless 
D=Diameter, fracture spacing, ft 
h=Reservoir thickness, ft 
k=Homogeneous reservoir permeability, md 
kf=Fracture permeability of dual porosity models, md 
km=Matrix permeability, md 
kx=Permeability in the Xdirection, md 
ky =Permeability in the Ydirection, md 
kz =Permeability in the Zdirection, md 
l=Half of fracture spacing, ft 
l−1=Inverse Laplace space operator 
L=General fracture spacing, ft 
m(p)=Pseudopressure (gas), psi2/cp 
Pi=Initial reservoir pressure, psia 
Pwf=Wellbore flowing pressure, psia 
P_{f}=Fracture pressure psia 
P_{flD}=Dimensionless pressure in the fracture 
Pm=Matrix pressure psia 
P_{mlD}=Dimensionless pressure in the matrix 
q_{lwD}=Dimensionless well rate based on matrixfracture 
qg=Gas rate, Mscf/day 
Q=Cumulative production, STB 
r=Radial geometry coordinate 
rw=Wellbore radius , ft 
s=Laplace space variable 
S=Skin dimensionless 
t=Time, hrs 
t_{DXf}=Dimensionless time coordinate 
t_{DA}=Dimensionless time based on fracture matrix geometry 
T=Temperature, 
y_{D}=Dimensionless reservoir length (rectangular geometry) 
x_{D}=Dimensionless reservoir length (rectangular geometry) 
z=Coordinate, zdirection (matrix) 
zD=Dimensionless coordinate, zdirection 
xw =XCartesian coordinates of the production point 
yw =YCartesian coordinates of the production point 
zw =ZCartesian coordinates of the production point 
f(s)= Relation used in Laplace space to distinguish matrix geometry types 
Abbreviation 
PEDNA: Pressure Equivalent of Density Weighted Average 
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 unsteadystate 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. 
CincoLey et al. [3] developed the concept of finite flow capacity and applied semi analytical approach to illustrate the importance of finite fracture when the FCD<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 constantrate flow test or a pressurebuildup test, depicting vertical hydraulic fracture model in an infiniteacting reservoir. In their study, they introduced a relationship between dimensionless time and pressure behavior which depends on time, and dimensionless fracture conductivity, FCD: 
Bennett et al. [4] established finite conductivity typecurves 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, pseudoradial, second linear and pseudosteady 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 loglog pressure derivative diagnostic plot and is used to determine fracture halflength, 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 finiteconductivity 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 loglog pressure derivative diagnostic plot and is used to determine the fracture conductivity [7]. 
Thirdly, pseudoradial 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 nonconventional 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 
• 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=L_{f1} and secondary fracture=L_{f2}. 
• 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/Trilinear 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: 
Dimensionless pressure: 
(2) 
Dimension time: 
(3) 
Dimensionless length coordinate (L), assuming isotropic properties 
(4) 
(5) 
(6) 
Where and are lengths of primary and secondary fractures, =distance along the fracture path. 
For isotropic system: 
(7) 
(8) 
The dimensionless variable rescaling the anisotropy system to an equivalent isotropic system is given as; 
(9) 
h=reservoir thickness, k_{x} and k_{y} are permeabilities along x and y axis If the fractures are of same length, then 
=equivalent fracture length 
The dimensionless fracture conductivity is defined as: 
(10) 
Interporosity flow parameter: 
(11) 
as defined by Warren and Roof, 1963 
Dimensionless storativity: (11a) 
For fracture one (Primary fracture), the diffusivity equation is given as: 
(12) 
And for fracture two (Secondary fracture) 
(13) 
If the fractures are of same length, then 
First the diffusivity equation for the matrix is given as: 
(14) 
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: 
(15) 
Initial boundary: 
(16) 
Outer boundary 
(17) 
Resolving the matrix diffusivity, equation (14) into dimensionless form using equation (2, (3) and (6), we have: 
(18) 
(19) 
Therefore 
(20) 
Where interporosity flow parameter [8] 
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: 
(21) 
Resolving the equation in dimensionless form using equation (2), (3) and (6) 
(22) 
Substitute foralong x axis, we have: 
(23) 
For secondary fracture resolve in y axis from equation (13) 
(24) 
Converting the equation dimensionless form using equation (2), (3) and (6) 
(25) 
Substituting for along y axis 
(26) 
Assuming the flux is uniform along the fracture and the pressure at the wellbore is a summation of pressure along the fractures segment, hencex axisy axis 
Combining equation (23) and (26) 
(27) 
Resolving in x axis 
(28) 
See full detail at the Appendix. 
The general solution combining the matrix and fractures differential equations with different BCs is given as: 
(29) 
Resolving the above equation by differentiating and applying BC, the final solution for crossform fractures model is given as: 
(30) 
(31) 
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. 
At x > 4.5 and coth(x) = 1 
Hence 
if 
Therefore 
(32) 
Where 
Case (I) bilinear flow: 
If ω=0 (33) 
(34) 
If 
Then 
Therefore 
If 
The general solution is given as; 
(35) 
This equation is due to bilinear flow period. 
However considering the assumptions, the flow regime is limited by; 
(35a) 
Therefore equation (35) is limited to: 
Case (II) linear flow: 
ω=1 (36) 
Therefore 
(37) 
This is the general solution for; 
This equation is due to the linear flow period with assumption limiting to t. 
Equation (38) is limited to the region: (38a) 
Case (III)radial homogeneous flow: For homogenous case, 
Recall that for matrix slab 
(39) 
Therefore 
(40) 
Substitute into equation 
(41) 
If 
(42) 
The general solution is given as; 
(43) 
This equation represents the homogenous phase. 
Also this equation is limited by 
(43a) 
Case (IV)trilinear flow: 
If and 
Asymptons 
Recall that for matrix slab 
If 
Therefore 
Recall 
Substitute equation 44 and the asympton 
into above equation 
(45) 
Substitute into equation 30 
(46) 
If 
(47) 
Converting to time dependent function using Laplace inverse 
(48) 
Where = Gamma Function 
Γ0.875=1.456 
Γ 0.75=1.225 
Therefore 
(49) 
This equation is due to Trilinear flow period. Also this equation is limited by 
(49a) 
Summary of the matrix and fractures diffusivity PDEs and the generated equations for each flow regimes is shown in Figures 6 and 7. 
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 buildup 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. 
(50) 
Svjetlana [7] Where: 
w – Fracture width 
we – Equivalent fracture width 
– Fracture porosity, fraction 
– Equivalent fracture porosity, fraction 
Fracture permeability is the function of the dimensionless fracture conductivity. 
Svjetlana [7] Where: 
F_{CD} – Dimensionless fracture conductivity 
k – Formation permeability 
x_{f} – Fracture halflength 
w – Fracture width 
Equivalent fracture permeability 
[7] Where 
w_{fe} – equivalent fracture width 
Summary of all reservoir, fracture and fluid properties used for the modeling are listed in 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 (Table 4); 
• Constant rate solution 
• Constant Pressure solution 
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. 
Five flow regions: Linear  Pseudoradial  Transition  Bilinear Trilinear were identified with this model. Region 1 which is linear due to transient flow only in the fractures. Region 2 is the response for a homogeneous reservoir which is dominated by transient matrix drainage and is the transient flow regime of interest. At this point, a small transition region dominated by a mix of linear and bilinear flow effect (region 3) [10,11]. 
Region 4 is bilinear flow and occurs when the matrix drainage begins simultaneously with the transient flow in the fractures. 
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 9 shows the regions with the dimensionless derivative rate generated from the developed fracture flow equations (35), (38) and (43) with limitation using equation (35a), (38a) and (43a). To visualize the effect of the fractures conductivities on the number of flow regions on the crossform fracture, F_{CD} ranges from 1.0 to 1000 mD was simulated. 
Result from Figure 10 support the five flow regions identified with the new developed fracture flow mathematical equations and flow model. Without the limiting equation of (35a), (38a) and (43a), only three flowing regions: LinearTransitionLinear can be identified as shown in Figure 11. Likewise, the result from sensitivities on gas production rate as shown in Figure 12 which indicates the loglog plot of gas rate versus time depicts three flow regions such as Linear BilinearLinear for F_{CD}<5 and two flow regions LinearLinear 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: LinearBilinear 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 we 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 LinearBilinearLinear 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 fracturematrix flow model applying ordinary laplace transform model OLTM to formulate linear, bilinear and trilinear 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: 
(i) 
Dimension Time: 
(ii) 
Dimensionless Length Coordinate (L) 
Assuming Isotropic Properties 
(iii) 
(iv) 
(v) 
For isotropic system: 
(vi) 
(vii) 
(viii) 
The dimensionless fracture conductivity is defined as: 
(ix) 
Interporosity flow parameter: 
(x) 
as defined by Warren and Root, 1963 Dimensionless storativity: 
(xa) 
For fracture one (Primary fracture) the diffusivity equation is given as: 
(xi) 
For fracture two {Secondary fracture} 
(xii) 
If the fractures are of same length, then 
First the diffusivity equation for the matrix is given as: 
(xiii) 
For the matrix and fracture interflowing period, the diffusivity equation is same for both fractures, therefore the following BC are applicable 
Initial condition: 
(xiii(a) 
Initial boundary: 
xiii(b) 
Outer boundary 
xiii(c) 
Resolving equation (xiii) into dimensionless form using equation (i), (ii) and (v) 
(xiv) 
Therefore 
(xv) 
Where [Warren and Root Interporosity flow parameter] 
Introducing the boundary conditions xiii(a), xiii(b), xiii(c) into equation xv 
(xvi) 
Since dimensionally 
Therefore 
(xvii) 
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 
Taking into consideration the inner boundary condition 
Substitute into equation xviii to obtain B 
Therefore B = 0 
Substitute B into equation xvii to obtain A with outer boundary condition. 
Outer boundary condition 
Therefore 
(xi) 
Therefore the general solution for the matrix flow solution is given as: 
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: 
(xxi) 
Resolving the equation in dimensionless form using equation (i), (ii) and (v) 
Substitute for l_{f1} = 4R^{2}cos^{2}θ along x axis, we have: 
(xxii) 
For secondary fracture resolve in y axis from equation xxii 
Converting the equation dimensionless form using equation (i), (ii) and (v) 
Substituting for l_{f2} = 4R^{2}sin^{2}θ along y axis 
(xxiii) 
Assuming the flux is uniform along the fracture and the pressure at the wellbore is a summation of pressure along the fractures segment, hence x axis = 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; 
Where 
Substitute into the above equation 
From Warren and Root interporosity flow parameter 
and 
Therefore 
Since the flow region modelling is within the fracture tips, therefore from equation vi: 
From dimensionless length coordinate (vi) and (vii) 
(xxiv) 
Mathematically 
Fracture (i) (Primary) 
 Fracture (ii) (Secondary) 
Invariably 
Substitute into equation (xxiv) 
Case a 
If the crosssectional area of the wall face is far away the fracture face, 
(xxv) 
Case b 
Assuming the crosssectional area of the fracture wall face 
Substitute into equation. Also apply equation (vi) and (vii) 
(xxvi) 
For case a 
Resolving this equation in Laplace form, 
The crosssectional area of the wall face 
Taking into account the boundary conditions 
Initial BC 
Inner BC 
(Craig 2006) formulate inner boundary condition describing transient flow in a finite conductivity fracture oriented along x axis. 
The dimensionless Laplace domain is given as r. 
(xxviii) 
For 
It is also written as 
Where the dimensionless variables are defined as 
(x,5) = Laplace domain flow rate per unit length into fracture 
q_{w} = Total well flow rate. 
K_{f?} = fracture permeability 
Therefore for constant rate, the inner BC for fracture face is given as: 
(xxix) 
and the outer BC for no flow through the fracture hp in Laplace form is given as: 
(xxx) 
Substitute the initial BC into equation (xxv) 
(xxxi) 
Recall equation (xx) 
Differentiate w.r.t Z_{D} 
Substitute into equation (xxxi) 
Where 
(xxxii) 
The characteristic equation for this differentiate equation and its roots are 
The General solution for this differential equation is given as: 
(xxxiii) 
Differentiate w.s.t. X_{D} 
(xxxiv) 
Applying the inner BC, 
Substitute into equation above; 
Apply outer BC; 
Introduce into equation xxxiv 
Substitute A and B into equation xxxiii 
All the Wellbore condition XD = 0 
Also in Laplace domain, the constant pressure case at the Wellbore can be obtained from the solution of the constant rate using the equation. 
Therefore 
Where 
(xxxv) 
This is general solution for 2w fracture connecting at the Wellbore. This equation can be inverted to obtain time dependant solution using Laplace inversion such as Stehfest’s inversion algorithm. 
At x > 4.5 coth(x) =1 
Hence 
Therefore 
Where 
case i 
If ω = 0 
Therefore 
And 
Hence 
Therefore 
Converting to time dependent function using Laplace inverse 
Where Γ = Gamma Function 
Γ 0.75 = 1.225 
(xxxvii) 
This equation is due to bilinear flow period. 
However considering the assumptions, the flow regime is limited by; 
Converting by Laplace inverse function 
Therefore 
 Condition (1) 
Also 
Therefore equation (xxxvi) is applicable if 
Case (ii) 
ω =1 
m⇒1+ 0 
m⇒1 
Therefore 
Where 
(xxxviii) 
This equation is due to the linear flow period with assumption limiting to t. 
Equation (xxxviii) is limited to the region 
Case (iii) 
f (s)⇒1.0 , This is for homogenous case. 
Recall that for matrix slab 
Therefore 
Substitute into equation (xxxv) 
If 
Converting to time dependent function using Laplace inverse function. 
This equation represents the homogenous phase. 
Also this equation is limited by 
Case (iv) 
If f (s)⇒0 and ω is a function of the Laplace parameter as defined below: 
Asymptons 
Recall that for matrix slab 
If 
Therefore 
(xxxx) 
Recall 
Substitute equation 35 and the asympton into above equation 
(xxxxi) 
Substitute into equation 
(xxxxii) 
If 
Converting to time dependent function using Laplace inverse 
Where Γ = Gamma Function 
Γ 0.875 = 1.456 
Γ 0.75 = 1.225 
Therefore 
(xxxxiii) 
This equation is due to trilinear flow period. Also this equation is limited by 
For case b 
Assuming the crosssectional area of the wall face 
And from equation xxvi 
Taking into account the boundary conditions as in case aand applying the same steps from equation xxvi to xxxvii, the resolve equation is given as : 
Where 
The general solution for this differential equation is given as: 
Differentiate w.s.t. XD 
Applying the inner and outer BC, the Laplace solution is given as: 
This is general solution for 2w fracture connecting at the Wellbore. This equation can be inverted to obtain time dependant solution using Laplace inversion such as Stehfest’s inversion algorithm. 
At x > 4.5 coth(x) = 1 
Hence 
Hence 
Where 
Case i 
If ω = 0 
Therefore 
Therefore 
Converting to time dependent function using Laplace inverse 
(xxxxiv) 
This equation is due to bilinear flow period. The flow regime is limited by; 
Case (ii) 
ω = 1 
m⇒1 
Therefore 
Where 
(xxxxv) 
Equation (xxxxv) is limited to the region 
Case (iii) 
f (s)⇒1.0 , This is for homogenous case. 
Recall that for matrix slab 
Therefore 
Substitute into equation (xxxv) 
Converting to time dependent function using Laplace inverse function. 
((xxxxvi) 
This equation represents the homogenous phase. 
Also this equation is limited by 
Case (iv) 
If f (s) ⇒ 0 and w is a function of the Laplace parameter as defined below: 
Asymptons 
Recall that for matrix slab 
Recall 
(xxxxvii) 
Converting to time dependent function using Laplace inverse 
This equation is due to trilinear flow period. Also this equation is limited by 
References 

Table 1  Table 2  Table 3  Table 4 
Figure 1  Figure 2  Figure 3  Figure 4 
Figure 5  Figure 6  Figure 7  Figure 8 
Figure 9  Figure 10  Figure 11  Figure 12 
Figure 13  Figure 14  Figure 15  Figure 16 