## Abstract

The shale gas experiences many different spatial scales during its flow in the reservoir, which will engender different flow mechanisms. In order to accurately simulate the production performance of shale gas well, it is essential to establish a multi-continuum model for shale gas reservoir. Based on the geometrical scenario of multistage horizontal well fracturing, this paper builds up a triple-continuum model incorporating three systems: matrix with extremely low permeability, less permeable natural fractures and highly permeable hydraulic fractures.

#### Authors:

Wei Zhang • Jianchun Xu • Ruizhong Jiang

China University of Petroleum (East China), Qingdao 266580, Shandong Province, China

Received: 21 November 2015 / Accepted: 27 August 2016

© The Author(s) 2016.

This numerical model employs Langmuir adsorption equation to present the influence of desorption gas in matrix and considers the Klinkenberg effect in matrix and natural fractures by adjusting the apparent permeability. The solution of this model is achieved using implicit scheme. Eventually, this model is applied on the single well production situation in a synthetic reservoir, production decline curves and cumulative production curves are obtained, then the sensitivity analysis is made on various kinds of parameters; thus, the influences of these parameters on production rate are obtained: The gas rate will rise with the increase in hydraulic fracture half-length, meshing size, Langmuir volume and Langmuir pressure, but with the decrease in hydraulic fracture spacing.

## Introduction

Nowadays, shale gas production has been making a significant contribution to the world’s gross energy supply (Li 2009) and accounts for over 50 % of American natural gas production (Montgomery et al. 2005). The exploration shows that China is also abundant in terms of the total recoverable shale gas reserves (Hu et al. 2010). Shale gas reservoir distinguishes itself by its extremely low permeability and low porosity, which engenders great difficulty on the exploitation and production of shale gas (Curtis 2002; Li et al. 2015). In addition, the pore size in shale matrix is only about 2–50 nm and 0.1–5 μm for the natural fractures (Sondergeld et al. 2010). Instead of the conventional method, the effective and economic development of shale gas reservoir requires the formation to be hydraulically fractured, and the multistage hydraulic fracturing method has been prevalently employed in shale gas reservoir (Ozcan et al. 2014).

Naturally fractured reservoir (NFR) can be defined as a reservoir contains a connected network of fractures created by natural processes and proved to have an effect on fluid flow (Javadpour 2009; Sun et al. 2016). NFR contains more than 20 % of the world’s hydrocarbon reserves (Sarma and Aziz 2006). Shale gas reservoir generally contains natural fractures; thus, shale gas reservoir can be classified into NFR. In shales, natural fractures provide permeability and the matrix provides storage for most of the gas (Bello and Wattenbarger 2010).

The gas molecules are stored by a combination of compression in the pores and adsorption on the surface of the solid shale matter (Bello and Wattenbarger 2008). It is believed that compared to the induced fractures, the permeability of natural fractures is too small to be considered. Yet this is proved to be wrong because the natural fracture network can actually enhance the productivity of the reservoir greatly (Cipolla et al. 2010); thus, the local natural fractures shall not be ignored.

From the nanoscale matrix to the hydraulic fractures whose width is several millimeters, shale gas experiences many different spatial scales during its flow in the reservoir (Sheng et al. 2012; Guo et al. 2015). According to different spatial scales, the flow of gas in the shale formation will result in many different mechanisms (Jiang and Wang 2014; Xu et al. 2015). As the natural fractures and the hydraulic fractures present a big difference in terms of fracture conductivity and connectivity, it is more realistic to assume fractures having different properties (Hassan and Wattenbarger 2011).

**Fig. 1** Illustration of gas flow from nanoscale to macroscale.

Commonly, the dual-porosity model is often used for the simulation of oil and gas reservoir, but only two continuums are unable to fully describe the multi-mechanism and multi-scale gas flow in shale formation (Cheng and Dong 2012). In order to accurately simulate the complicated shale gas flow and production, multi-continuum models shall be built up (Zhu and Zhang 2013). Furthermore, some microscopic mechanisms, such as desorption and diffusion effect, should be considered in the model as well (Ozkan et al. 2010); thus, triple-continuum model begins to emerge (Fig. 1).

The first triple porosity was introduced by Abdassah and Ershaghi (1986), and they divided the matrix to have different properties with single fracture. Then Al-Ahmadi and Ershaghi (1996) first assume the fractures to have different properties, their model was presented using a radial system, and the natural fractures and the hydraulic fractures can both feed the well. Drier (2004) improved the triple-continuum model originally proposed by Al-Ahmadi and Ershaghi (1996) by considering transient flow condition between natural fractures and hydraulic fractures.

Bello and Wattenberger (2008, 2009, 2010) applied the triple-continuum model to analyze the production rate in horizontal well in tight fractured reservoirs. Ozkan et al. (2009) proposed a trilinear model comprising three continuous media: finite conductivity hydraulic fractures, dual-porosity inner reservoir between the hydraulic fractures and outer reservoir beyond the tip of the hydraulic fractures. Apaydin et al. (2012) examine the effects of matrix micro-fractures on effective matrix permeability of a dual-porosity medium. They stated that matrix micro-fractures accelerate production by providing earlier and more effective contribution of the matrix into flow rates. Therefore, the multi-continuum model for shale gas reservoir commonly comprises three media: matrix, natural fractures and hydraulic fractures, which is also called dual-fracture model (Al-Ahmadi and Ershaghi 1996).

The triple-continuum model has undergone great development recently, and models for reservoirs with sophisticated geometries and conditions have been built up in previous work. The findings show that the triple-continuum model can capture the reservoir heterogeneity very well (Al-Ahmadi 2010). Nonetheless, most previous triple-continuum models were solved using analytical method instead of numerical approach to predict the transient pressure in the reservoir. In addition, diffusion mechanism, which is proved to have significant impact on gas flow, was rarely considered in previous models. In this work, a new numerical triple-continuum model incorporating desorption and diffusion is proposed and applied for the production simulation in fractured shale gas reservoir.

## Incorporated mechanisms

#### Adsorption and desorption of shale gas

Shale gas can exist as free gas phase or as adsorbed gas on solids in pores, and the methane molecules are mainly adsorbed to the carbon-rich components, i.e., kerogen (Lu et al. 1995). As observed, the adsorbed gas accounts for a significant fraction of gas reserves in place. With the drawdown of the formation pressure, the adsorbed gas is released from the surface of solids and becomes free gas, contributing to the total production. In this work, the mass flux of the adsorbed gas in a micro-unit within a short time is described as:

Here the Langmuir adsorption model is employed to represent the adsorbed gas content V_{E}:

Where V_{L} stands for the Langmuir volume, p_{L} is the Langmuir pressure, and p is the reservoir pressure.

#### Klinkenberg effect

In low-permeability shale gas reservoirs with small pore space, the Klinkenberg effect may alter the permeability significantly, especially in low reservoir pressure conditions (Wu et al. 1998). In this work, Klinkenberg effect is incorporated into the gas flow equations by modifying the apparent permeability of gas as a function of pressure (Wu et al. 2009):

where k_{a} is the apparent permeability of the gas phase; k_{i} is the constant, equals to the gas permeability in large pores without Klinkenberg effect; b is the Klinkenberg factor. Although the Klinkenberg factor may change with gas nature and pore size, we adopt b as a constant in our simulation, and its value can be calculated as (Yao et al. 2013):

## Model construction

#### Model assumptions

Before building the triple-continuum model, some assumptions are made based on the characteristics of each medium:

- Three media are incorporated in the model: matrix, less permeable natural fractures and more permeable hydraulic fractures.
- As shown in Fig. 2, the hydraulic fractures perpendicular to the horizontal well are discretely distributed, and the gas flow in hydraulic fractures is considered as one-dimensional flow, whereas the matrix system and the natural fracture system are continuously distributed over the whole rectangular reservoir, and the flow in these two media is considered as two-dimensional flow.
- Flow is sequential from one medium to another: from matrix to natural fractures to hydraulic fractures.
- The desorption of adsorbed gas is only considered in matrix, and the Klinkenberg effect exerts influences on both matrix and natural fractures.
- The rock is incompressible, and the porosity is seen as a constant.
- The gas is regarded as compressible real gas, and the gas viscosity and gas compressibility factors are both the function of pressure.

**Fig. 2** Geometric model of the fractured shale gas reservoir.

#### The flow equation of matrix

The continuity equation of matrix is given by:

q_{m} is the inter-porosity flow term from matrix to natural fractures, which is expressed as:

where α_{mf} is the shape factor from matrix to natural fractures, A_{mf} represents the interface area between matrix and natural fractures per unit volume, and Lf represent the characteristic length of natural fracture, which can be seen as its average spacing. The natural fracture spacing in the shale formation is generally 0.05–10 m (Bustin et al. 2008).

In the continuity equation of matrix, when considering Klinkenberg effect, the apparent permeability in the kinetic equation shall be adjusted:

And the gas density can be expressed as:

Substituting Eqs. (6)–(10) into Eq. (5), then we obtain:

Taking the transformation we can obtain:

Substituting Eq. (12) into Eq. (11), we obtain:

Equation (13) is final form of the flow equation of matrix system, and its initial and boundary conditions are:

#### The flow equation of natural fractures

The continuity equation in natural fracture is:

where q_{f} is the flow term from the natural fractures to the hydraulic fractures:

α_{fF} is the shape factor from natural fractures to hydraulic fractures, and it is defined as the interface area of natural fractures and hydraulic fractures per unit volume divided by the characteristic length of hydraulic fracture.

When considering Klinkenberg effect in the natural fractures:

Substituting Eqs. (19) and (20) into Eq. (16), we obtain:

#### The flow equation of hydraulic fractures

Unlike the matrix and the natural fractures, hydraulic fractures are discretely distributed with a certain spacing over the reservoir; consequently, the flow in the hydraulic fractures is viewed as one-dimensional flow. Moreover, neither the desorption nor the Klinkenberg effect is considered in the hydraulic fractures; thus, its flow equation is:

where q_{w} is the flow term from hydraulic fractures to the horizontal well, V_{i,j} is the volume of the grid at which the hydraulic fractures intersect with the wellbore, and q is the gas production rate. Figure 3 is an illustration of flow in hydraulic fracture grids.

Note that q_{f} is the inter-porosity flow from NF to HF and this is applied to all hydraulic fracture grids, because each HF grid receives flow from NF system. q_{w} is the flow term from HF to the horizontal well, it is only applied to the grid where the HF intersects with the horizontal well, and the production rate is assumed to be constant.

**Fig. 3** Illustration of flow in hydraulic fractures.

Substituting all the related equations into Eq. (24), we obtain the flow equation of hydraulic fractures, followed by its initial and boundary conditions:

So far we have derived the differential equations for the three media [Eqs. (13), (21) and (29)] together with their initial and boundary conditions; therefore, the entire triple-continuum model is established.

#### Solution of the mathematical model

This triple-continuum model is solved using finite difference approach, and we discretized the differential equations in each medium via implicit scheme. Within the same time step, we first solve the hydraulic fracture equation to obtain the hydraulic fracture pressure p_{F} and then put it into the natural fracture equation to derive p_{f}, and eventually p m is obtained by solving the matrix equation. By substituting p_{F} into Eq. (28), the production rate of each time cycle is calculated.

Based on this solving algorithm, a program is developed for our simulator. In the program, the gas compressibility factor is calculated by iterative calculation and an empirical equation is employed to calculate the gas viscosity. The detailed calculation procedures for gas properties are shown in Appendix A and B.

### Model application and discussion

#### Model validation

To verify the correctness of this model, we compare the gas rate curve obtained from our model with the prediction from an existing model proposed by Wang (2015). Here we neglect the desorption and Klinkenberg effect and set the parameters of the reservoir according to their paper (reservoir length = 100 m, reservoir width = 100 m, HF height = 100 m, HF half length = 40 m, reservoir permeability = 10^{−3}mD, reservoir porosity = 0.1, reservoir initial pressure = 15 MPa, bottom-hole pressure = 5 MPa). Figure 4 presents the result of the comparison.

**Fig. 4** Contrast between this study and the result proposed by Wang (2015).

From the results, we can notice that the gas rate curve predicted using our model declines slower than the reference curve in the early stage of the production, it may due to the different calculation of the gas compressibility factor and gas viscosity in the program, and after 400 days of production time, the two gas rate curves can match with each other with negligible disparity. Therefore, despite minor difference from the existing model in the early stage, our model is quite reliable in general.

#### Reservoir description and model setup

A synthetic shale gas reservoir is employed to demonstrate the application of our proposed model. The reservoir is drilled by a horizontal well and followed by multistage fracturing, and the basic parameters of the shale gas reservoir and the horizontal well are listed in Table 1.

**Table 1** Basic parameters for the reservoir and horizontal well.

This hypothetical shale reservoir is nominally assumed to be 1500 m deep, and a horizontal well with six hydraulic fractures is incorporated in the model. The size of the reservoir is set to be 1300 m × 600 m × 20 m, the grid size is 75 × 33 × 1 in X, Y and Z direction, and the flow on the vertical direction is neglected. The width of the grids containing hydraulic fractures is set to be 2 m, and the mesh size near the fractures and horizontal well is smaller and finer in order to increase the accuracy. The actual width of hydraulic fracture is assumed to be 5 mm; thus, in our numerical model, the effective permeability for HF grids is calculated using the equation: k_{eff} = k_{F} W_{F}/W_{grid} = 20md (Rubin 2010). This synthetic reservoir is set for the following sensitivity analysis, and the model setup maybe moderately adjusted for the convenience of sensitivity analysis.

#### The influence of meshing size on gas rate

In order to verify the stableness of this mathematical model, the near wellbore meshing size is altered for different simulation processes, and its influence on gas rate is examined and illustrated in Fig. 5.

**Fig. 5** Production performances for different values of meshing size.

As we can see from Fig. 5, although it shows difference between the three curves at the beginning of the simulation, this difference disappears gradually with time elapsing and only small difference can be seen in the later stage of the simulation. Consequently, the gas rate is insensitive to the change of meshing size, which indicates good stableness of the proposed model as well (Wang 2015).