## A Fully Three Dimensional Semianalytical Model for Shale Gas Reservoirs with Hydraulic Fractures

Abstract:

Two challenges exist for modeling gas transport in shale. One is the existence of complex gas transport mechanisms, and the other is the impact of hydraulic fracture networks. In this study, a truly three dimensional semianalytical model was developed for shale gas reservoirs with hydraulic fractures of various shapes. Using the instantaneous point source solution, the pressure are solved for a bounded reservoir with fully 3D, partially penetrated hydraulic fractures of different strike angles and dip angles. The fractures could have various shapes such as rectangles, disks and ellipses.

Authors:

Yuwei Li^{1,*,†} , Lihua Zuo^{2,*,†}, Wei Yu 2 , * and Youguang Chen^{3}

^{1}Department of Petroleum Engineering, Northeast Petroleum University, Daqing 163318, China. ^{2}Harold Vance Department of Petroleum Engineering, Texas A&M University, 3116 TAMU, College Station, TX 77843, USA. ^{3}Department of Petroleum and Geosystems Engineering, University of Texas at Austin, Austin, TX 78705, USA; [email protected]

^{†}These authors contribute equally to this manuscript.

Received: 26 January 2018; Accepted: 11 February 2018; Published: 15 February 2018

The shale gas diffusion equations considers complex transport mechanism such as gas slippage and gas diffusion. This semianalytical model is verified with a commercial software and an analytical method for single fully penetrated rectangle fracture, and the production results of shale gas are consistent. The impacts of fracture height and strike angles are investigated by five systematically constructed models. The comparison shows that the production increases proportionally with the fracture height, and decreases with the increase of strike angles. The method proposed in this study could also be applied in well testing to analyze the reservoir properties and used to forecast the production for tight oil and conventional resources.

## Introduction

Shale formations are characterized by low permeability, nanopores, high total organic carbon, and natural fractures. In shale, the gas content is mainly composed of free gas and adsorbed gas. Because of its low permeability, shale gas did not gain much success until the application of a horizontal well drilling and hydraulic fracturing technique.

Currently, shale gas accounts for about half of the natural gas produced in the United States of America (U.S.A.) and that ratio will increase to 70% in 2040 according to a prediction by the U.S. Energy Information Agency [1,2] To simulate the gas production process in shales, many complex gas transport mechanisms such as gas slippage and desorption and the presence of natural fractures need to be considered, which makes the simulation of shale gas production a challenging work.

There are mainly three gas coexisting transport mechanisms during the gas production process, including gas diffusion, gas slippage and gas desorption [3]. The complexity of these mechanisms makes the classical free gas diffusivity equations inadequate to model the gas transport in shales. Despite the availability of some existing models for simulation of gas transport in shales [4,5], most of them do not fully capture the complex mechanisms, and a comprehensive gas transport equation is still needed.

Another challenge presented when dealing with shale gas transport simulation is the modeling of complex fracture patterns. The dual-porosity or dual permeability model [6,7] is widely used in numerical simulations of the fractures, but the model is usually over-simplified compared to the actual field complexities, and the resolution of this approach is not enough to capture the details of fractures. In that respect, discrete-fracture models (DFM), using finite-volume or finite-difference methods, have been developed [8].

To be compatible with the complex geometries of fractures, such as nonplanar fractures and fractures with variable aperture, reservoir simulators using unstructured grids have also been designed [9–12] in order to explicitly model the fractures. However, the use of unstructured grids leads to high computational cost so that it is only used in limited applications in real field studies. As a solution, a technique that directly incorporates complex fractures in a conventional structured grid was developed [13,14], and Moinfar et al. [15] named it embedded discrete fracture model (EDFM). The EDFM was further developed for 3D reservoir simulation by Moinfar et al. [16]. Recently, the EDFM has been widely used for naturally and hydraulically fractured reservoirs [17–20]. Despite the accurate capture of the fracture geometries, the EDFM method needs extra computation time to calculate the extra transmissibility caused by the non-neighborhood connections between fractures and matrix cells.

Assuming the rock properties are homogeneous, there are many analytical and semianalytical models derived to efficiently simulate the gas transport in shales. Gringarten et al. [21] used Green’s functions to solve unsteady flow problems for a wide variety of reservoir flow problems. Cinco-Lay and Samaniego [22] presented a new technique for analyzing pressure transient data for wells intercepted by a finite conductivity vertical fracture.

Wan and Aziz [23] described a new semianalytical solution for horizontal wells with multiple fractures with different strike angles and fully/partially penetrating the formation in the vertical direction. With that solution, the authors calculated the well index combined with a numerically computed gridblock pressure.

Lin and Zhu [24] studied the well performance for fractured horizontal wells in an infinite slab reservoir, by using the instantaneous solutions of plane sources. Zhou et al. [25] designed a semianalytical model to simulate the gas production in Barnett shales with fully penetrated planar fractures. Yu et al. [26] extended Zhou’s model to simulate well production for reservoirs with fully penetrated non-planar fractures, which have varying fracture width and fracture permeability. Their model was verified with a numerical reservoir simulator for single fracture case before being applied to several case studies based on Marcellus shale properties.

Yang et al. [27] developed a robust and comprehensive model for real gas transport in shales with complex non-planar fracture network, with gas diffusion and gas slippage effect considered. He et al. [28] proposed a semianalytical method to diagnose the locations of underperforming fully penetrated hydraulic fractures through pressure transient analysis in tight gas reservoirs. Wang et al. [29] constructed a semianalytical model for simulating shale gas reservoirs with complex fully penetrated fractures.

In most of the above studies, due to the complexity of deriving the explicit pressure solutions for partially penetrated fractured in a bounded reservoir, the fractures are assumed to be fully penetrated rectangle shapes in a slab reservoir. Because coupling pressure and flow rate could bring an extra difficulty, most of the previous work concentrate on the pressure profile without being concerned about the production. In real situations, the fractures are most likely partially penetrated and have various shapes, such as disk and ellipse [30–32]. In this study, a novel semianalytical model is designed for bounded reservoir with multiple fully or partially penetrated fractures with various shapes such as rectangles, disks, and ellipses. First, the pressure solutions for fractures with dip angles or strike angles, and various shapes of fractures are derived, including rectangle, ellipse and disk shapes.

Based on the pressure solutions and the mathematical equation developed for shale gas reservoirs with proper assumptions, the semianalytical model is introduced and described. The new model is verified against a commercial reservoir simulator and an analytical model for a bounded reservoir with a horizontal well and three fully penetrated fractures. In order to investigate the impact of fracture heights and strike angles on gas production, this model is then used to model the shale gas production for five bounded reservoirs with different hydraulic fracture patterns, varying from fully penetrated to partially penetrated, from fractures without strike angles to fractures with different strike angles. To the best of the authors’ knowledge, this study is the first to design a semianalytical model for shale gas reservoirs with partially penetrated fractures of various shapes. The methodology proposed in this study could be expanded and implemented in existing reservoir simulation tools to investigate various fracture patterns.

## Materials and Methods

### Model Development

For gas transport from shale matrix into fracture segment, the diffusivity equation is revised for conventional gas reservoirs by considering the important gas transport mechanisms such as gas slippage, gas diffusion, and gas desorption as follows:

where ρ_{g} is gas density, φ is rock porosity, k is reservoir permeability, _{Cg} is the isothermal gas compressibility factor, p is pressure, µ_{g} is gas viscosity, Va is pore volume fraction of adsorbed gas, Kn is Knudsen number, α is a constant and close to 1,8αKn is used to describe the gas slippage effect [33], D_{g} is the Fickian diffusivity of gas component through the pore and used to describe the gas diffusion process, which might be the combination of three distinct mechanisms acting individually or simultaneously: bulk diffusion, Knudsen diffusion and surface diffusion [34]. δ is a dimensionless constrictivity factor (≤1), which accounts for variation of the pore size along its length caused by small pores and can be calculated by the empirical equation developed by Satterfield et al. [35]:

where d_{gas} is gas molecule diameter (methane molecular diameter is 0.38 nm) and d_{pore} is the pore diameter, τ is a dimensionless tortuosity factor (≥1), which can be estimated from porosity and gas saturation by [36]:

where S_{g} is gas saturation, K_{a} is the differential equilibrium partitioning coefficient of gas at a constant temperature, which is used to describe the gas desorption mechanism and defined as [37]:

where ρ_{a} is the adsorbed gas mass per unit shale volume (kilograms of adsorbed gas per cubic meter of solid). A general form of Brunauer–Emmett–Teller (BET) isotherm is used to model gas desorption effect [38], which is given below:

where ν( p) is the gas volume of adsorption at pressure p, p_{o} is the saturation pressure of the gas, ν_{m} is the maximum adsorption gas volume when the entire adsorbent surface is being covered with a complete monomolecular layer, C is a constant related to the net heat of adsorption, n is number of the adsorption layers. Note that when n = 1, Equation (5) will reduce to the standard Langmuir isotherm [39]. For the general form of BET isotherm, the differential equilibrium partitioning coefficient of gas can be expressed as:

For gas flow from fracture to wellbore, the non-Darcy flow effect is considered by using the Forchheimer [40] modification to Darcy’s law below:

where v is velocity, β is the non-Darcy Forchheimer coefficient, which can be determined using the correlation proposed by Evans and Civan [41] as given below:

where k_{f} is the fracture permeability.

### 3D Semianalytical Methods

Using the method of images and considering the geometry, we can generate the solution for a unit-strength, instantaneous point source in an infinite-slab reservoir with impermeable boundaries at z = 0 and h, which is given below [42]:

where ηx , ηy , ηz is the hydraulic conductivity defined as

, and M is the interest point and M0 is the source point. If we assume the hydraulic conductivity is isotropic, then Equation (11) could be rewritten as follows:

Using the superposition theory, we can get the pressure solution for rectangles with no strike angles (Figure 1a) as follows:

Rewriting Equation (13) we get:

where xw , yw , zw is the en fracture center. Calculating the integration with respect to Z gives:

**Figure 1.** Illustration of three fracture geometries: (a) Rectangle fracture plane which is parallel to the y‐z plane and x‐z plane; (b) Rectangle fracture plane with strike angle. (c) Rectangle fracture plane y-z plane and x-z plane; (b) Rectangle fracture plane with strike angle θ; (c) Rectangle fracture plane with dip angle θ.

Thus Equation (14) could be simplified to double integral expression, as follows:

### Slanted Rectangles Fracture with Strike Angle

Sometimes the fracture plane is not orthogonal to the x-y plane, but rather has a strike angle θ as plotted in Figure in Figure 1b. In this case, the pressure solution corresponding to Equation (14) will be derived as follows:

It is necessary to point out that “csc θ” term in Equation (18) is obtained when the plane integration is performed for each fracture plane. Applying the same technique as in Equations (15) and (16), we can simplify this to the following expression:

The corresponding bounded domain formula will become:

Equation (20) will be the primary equation to be used hereafter for rectangle fractures.

### Slanted Rectangles Fracture with Dip Angle

Next we study the cases in which the fracture plane is not necessarily orthogonal to the x-y plane, and instead the angle between the fracture plane and the x-y plane is set to be θ.

In this case, the pressure drop solution corresponding to Equation (14) will be as follows:

If θ=90, then Equation (21) will be simplified to Equation (20). The corresponding bounded domain formula will then become:

### Pressure Solutions for Rectangle, Ellipse, and Disk Fractures

To show the versatility of our method, and for the purpose of faster computation, the pressure solutions for three other shapes of fractures are constructed in a slab reservoir. Following the same logic as for the rectangle fracture in a bounded reservoir, the pressure solution of other shapes could also be derived. The pressure solution for a disk fracture in a slab reservoir is given as below:

Similarly, the pressure solution for an ellipse fracture in a slab reservoir is given as:

Three cases are designed with three different shapes of fractures. In all three cases, the depth of the slab reservoir is 30 m, and the center of all three fractures is at (0, 15). The dimensionless pressure results after 1000 days are plotted in Figure 2. The dimensionless pressure is defined as the actual minimum and maximum pressure being translated to [0, 1] range by the linear translation. Figure 2a is a fracture with rectangular shape, where the length of the rectangle equals 20 m and the height equals 10 m.

Figure 2b has a disk fracture, with radius equal to 10 m. Figure 2c has an elliptical fracture, with a long radius equal to 10 m and short radius equal to 5 m. Using Equation (17) and Equations (23) and (24), the pressure solution could be calculated. To concentrate on different pressure distribution contours caused by different shapes of the fractures, only the dimensionless pressure values are plotted. The dimensionless pressure solutions are shown in Figure 2.

Comparing Figure 2a–c, there are two observations to make. First, it shows that our pressure solutions capture the pressure contours accurately, with the boundary effect clearly seen from the top and bottom of the reservoir. Second, the fracture shape impact the pressure distribution thus the shale gas production significantly, so that it is worthwhile obtaining better measurement about the fracture shapes.

**Figure 2.** Illustrations of dimensionless pressure distributions for various shapes of fractures.(a) Rectangle fracture located at (−10, 10) × (10, 20); (b) Disk fracture, centered at (0, 15), and radius equals to 10 m; (c) Ellipse fracture, centered at (0, 15), with long radius a= 10, and short radius b = 5.

### Semianalytical Model for 3D Fractures in Shale Gas Reservoir

For the non-planar fracture geometry with varying fracture width and fracture permeability along fracture length, a semi-analytical model is developed to efﬁciently simulate the shale gas production.The semi-analytical model discretizes the non-planar fracture geometry into a number of small fracture segments in order to approximately describe the non-planar fracture geometry. Each fracture segment can have different fracture properties such as fracture length, width, and permeability to approximately represent the non-planar fracture geometry.

The semi-analytical model mainly consists of two parts for simulating shale gas production: gas ﬂow from shale matrix into fracture and gas ﬂow from fracture to wellbore. First, we used an analytical solution to approximately solve the revised diffusivity equation for gas transport from shale matrix into fracture segments. Second, gas ﬂow in each fracture segment is described by the non-Darcy ﬂow equation. For the ﬁrst part, an analytical solution is used to solve the gas diffusivity equation with the following assumptions:

(1) The reservoir is bounded;

(2) The reservoir is isotropic and homogeneous;

(3) It only allows gas flow from fractures to wellbore.

(4) There is no gravity effect.