Abstract

Shale reservoirs have some natural fractures with a certain density and connectivity. The basic percolation model of shale gas reservoir: the black oil model of gas-water phase is used as the basic model, and the dissolved gas is used to simulate adsorbed gas. Accurate description of natural fractures: random distributed discrete fracture model is used as the basic model to describe natural fractures. By comparing the calculation results of single medium (including random distributed discrete fracture model) and double medium model, the model for predicting shale gas productivity is optimized.

1. Introduction

China’s shale gas resource potential is enormous, but the reservoir matrix is very tight [1–4]. With the innovation of fracturing technology [5], volumetric fracturing horizontal well technology has become a key technology for the development of unconventional oil and gas reservoirs in China [6].

Microseismic monitoring technology shows a complex fracture network formed in the reservoir after volume fracturing. However, the microseismic data can only obtain the general shape of fracture reconstruction area and cannot obtain accurate fracture network information. Well test is an essential means to obtain reservoir information. Through well test data interpretation, key parameters affecting gas well production can be inversed and obtained, which is of great significance for reservoir evaluation, fracturing evaluation, and dynamic prediction [7–9].

At present, a lot of well test models have been put forward based on the assumption of linear flow because of the long linear flow in bulk fractured horizontal wells. Wattenbarger et al. [10] developed a linear flow model based on the assumption that matrix fluids flow linearly toward fractures, and many improved models have been based on this assumption. Bello et al. considered the matrix as a double medium model and improved the linear flow model [11]. Nobakht et al. [12] and Wu et al. [13–15] used this model to carry out productivity analysis. Based on Aiahmadi et al. [16] model and Bello Model, a three-hole model of horizontal well in volume fracturing was established by considering the secondary fractures formed after fracturing modification. Brohi et al. [17] established a three-linear flow model considering the effect of the unreformed zone. Assuming that the fluid flow in each zone is linear, Stalgorova and Mattar (2013) further considered the incomplete reconstruction between the stages of fracturing and developed a five-zone model. Jia et al. [18] have established the mathematical models of gas-water two-phase and oil-gas two-phase flow, respectively, considering the multiphase flow in shale reservoir. The abovementioned foreign models only analyze the productivity performance of horizontal wells with volume fracturing but not the well test analysis. Many oil and gas fields in China use well test to evaluate reservoirs. Zhang et al. [19] have studied the method of well test for straight fractured wells. Li et al. [20] and Wang et al. [21] have put forward the model of flat fracture interpretation by well test for multisection fractured horizontal wells. Wang et al. [22] further proposed a new well test interpretation deconvolution method. Zhu et al. [23] and Gao et al. [24] used the flat fracture model to analyze the well test interpretation while considering the factors of shale gas adsorption and desorption, negative skin effect, and the outer zone of modification. Given the heterogeneity of large-scale fractures, Zhang et al. established the numerical simulation method of discrete fracture seepage flow [25].

The gas seepage law of shale gas reservoir is very complicated, which is different from that of conventional low permeability reservoir [26–43]. In the reservoir, the gas flowing medium includes matrix pore, natural fracture, and artificial fracture. Artificial fractures are generally a few millimetres in size, and natural fractures are smaller than the size of artificial fractures, while matrix pores are smaller and can reach the nanometer scale, so conventional Duthie’s law is difficult to describe their flow patterns accurately. Besides, it is difficult to describe the gas desorption and matrix-to-fracture diffusion in shale gas reservoirs, which raises a difficult problem for predicting the productivity of shale gas reservoirs.

Considering that there is no percolation model for shale gas reservoir, the simulation model for shale gas reservoir is based on black oil model or dual-medium model. The black oil model is the most classical model in the numerical simulation of oil and gas reservoirs. However, it cannot involve gas desorption, so it is not suitable for shale gas reservoir. In this paper, the method of simulating adsorbed gas with dissolved gas was used to simulate the Langmuir isotherm adsorption curve with dissolved gas-oil ratio curve, and the oil phase was regarded as the nonflowing phase. The gas desorption velocity on the shale surface was assumed much faster than the gas flow velocity in the natural microfractures, thus ignoring the gas diffusion. Based on the double medium model of double porosity and single permeability, the desorption and diffusion of gas were considered, and the simulation results were compared with that of a single medium.

2. Development of the Black Oil Model

We first present the governing equations of mass conservation for each of the three components in the black-oil problem under saturated condition (i.e., for three phases). Following the introduction of black-oil terminology, we derived a pressure equation; then, the analogy between the developed black oil model and a compositional model formulation was established. We also discussed our approach to updating phase saturation in the black oil model. Next, we performed the analysis for the undersaturated condition. Finally, the numerical methods used to solve the governing equations were briefly discussed. Throughout this article, lowercase letter subscripts , , and denotes the water, oil, and gas phases, respectively, while water, oil, and gas components are, respectively, indicated by . The superscript refers to standard (surface) conditions.

2.1. Single-Medium Model

Black oil model is the most perfect and mature model for nonvolatile crude oil, and it is also the most widely used model (Lu et al., 2001). The black oil model considers the three-phase flow of oil, gas, and water, in which the gas phase exists in the form of free gas and dissolved gas.

2.1.1. Black Oil Model

The basic assumptions of the black oil model are as follows: (1) There are three phases of oil, gas, and water at most, and the seepage law follows Darcy’s law. (2) The oil-gas two-phase reaches the equilibrium state in an instant, the oil and water are not mutually soluble, and the gas is not soluble in water. (3) If there is dissolved gas, only consider the gas component dissolved in the oil component. (4) The percolation in the reservoir is isothermal. (5) Reservoirs and fluids are slightly compressible.

Then, the oil-gas-water three-phase black oil model can be expressed as:

The auxiliary equations are as follows:

where is the formation permeability, ; , , and are the relative permeability of oil, gas, and water; , , , and are the density of oil, dissolved gas, gas, and water, ; , , and are the pressure of oil, gas, and water phase, MPa; , , , and are the weight of oil, oil, gas, and water phase containing dissolved gas, , where ,, and ; is the depth from a certain datum, m; is the formation porosity; , , and are the oil, gas, and water saturation, respectively; is the time, s; is the oil-water capillary pressure, MPa; and is the oil-gas capillary pressure, MPa. Equations (1)–(6) is a complete oil-gas-water three-phase black oil model; it can be seen that the system is closed due to six unknowns—, , , , , and . So, by adding some boundary conditions and initial conditions, the equations can be solved.

Assuming that the boundary of the model is closed, the boundary conditions are as follows:

where , , and are the lengths of the , , and directions, respectively, m.

The initial conditions are as follows:

where is the original formation pressure, MPa; is the original water saturation, and is the original gas saturation. With these boundary and initial conditions (7)–(12), we need to linearize the (1)–(3) differential equations. At present, the most widely used linearization method is a difference discrete method. First, the partial derivative is replaced by the difference quotient; then, the partial differential equation is recurrence relation, and the corresponding algebraic equations are formed on the finite points in the solution domain; finally, the equations are solved by the fully implicit or semi-implicit method. The exact difference method is described in Section 2.3 of this chapter.

2.1.2. Adsorption Gas Treatment

The process of producing shale gas is mainly divided into two stages. The first stage is the extraction of free gas in natural fractures. With the decrease of formation pressure, the desorption of adsorbed gas from the matrix makes the development of shale gas enter into the second stage. Generally speaking, the gas output in the first stage decreases rapidly, and the gas output is relatively stable after entering the second stage, so the adsorbed gas is the most crucial factor to influence the gas production in the later stage of shale gas development. If the desorption rate of shale gas from the matrix to fracture is much faster than the flow rate in fracture, the desorption is relatively unimportant in the process of simulating shale gas production, so it can be assumed that the desorption process is instantaneous. With this hypothesis, shale gas adsorption can be simulated by gas dissolution in immobile oil. That is to say, at a given pressure, the amount of gas adsorbed by shale gas in the matrix is similar to the amount of gas dissolved in crude oil at the corresponding pressure, and the dissolved gas-oil ratio curve can simulate the Langmuir curve. The introduction of oil component needs to increase the porosity and modify the saturation of the reservoir and modify the relative permeability curve of gas and water.

When using dissolved gas to simulate the adsorption of gas, the total amount of water and gas in the black oil model should be equal to the total amount of water and gas in shale:

In the formula, and are the gas and water saturation, is the model porosity, and is the actual porosity. The above two types of the merger are:

If , then:

In the formula, is the oil saturation of the model. Equation (16) relates the actual porosity of the shale to that of the black oil model. Because the oil saturation value is chosen randomly, the oil phase saturation must be constant. When the formula (16) is combined with the formulas (13) and (14), the following are:

The relative permeability curves of gas and water can be adjusted by the equations (17) and (18): the relative permeability values corresponding to the actual and of the shale reservoir should be assigned to the corresponding and in the model, and it is considered that oil does not flow, which can be used to set the relative permeability of 0 or superviscosity to express. If =0.45, the real and modified infiltration curves are shown in Table 1.

According to the Langmuir isothermal adsorption equation, the amount of gas adsorbed on a unit volume of shale is:

In the formula, is the gas adsorption, ; is the Langmuir volume, ; is the pressure, MPa; and is the Langmuir pressure, MPa. Since the amount of adsorbed gas depends on the mass rather than the volume of the shale, another Langmuir equation is more suitable for describing the adsorption characteristics of the shale:

In the formula, is the gas adsorption, ; is the Langmuir constant, ; and is the Langmuir pressure constant, MPa^-1. In actual shale reservoirs, the amount of adsorbed gas per unit volume is:

In the formula, is the volume density of the shale, . In the model, the amount of gas dissolved in the same unit volume can be expressed as:

In the formula, is the dissolved gas-oil ratio, , and is the volume coefficient of crude oil formation. In order to ensure material balance, the volume factor must be constant. The combination of (21) and (22) is:

Equation (23) can be used to convert adsorbed gas in the shale matrix into dissolved gas in stagnant crude oil. The productivity of shale gas reservoir can be evaluated by the oil-gas-water three-phase black oil model.

Table 2 shows a comparison of adsorbed gas and dissolved gas-oil in a shale reservoir. During the calculation, the other parameters in equation (23) are as follows: =1500 kg/m³, =0.0244 m³/kg, =1.2, =0.3 MPa^-1, =0.36, and =0.45.

2.1.3. Discrete Fracture Model

Because the seepage law of fractured reservoir is quite different from that of the porous reservoir, the description of fractures in the fractured reservoir is very important to the simulation of development effect, and fracture permeability-dominated reservoirs are often more complicated than matrix permeability-dominated reservoirs. At present, Warren-Root model and discrete fracture network model (DFN model) are the main models to describe fractured reservoirs, as shown in Figure 1.

(a)

(b)

(a)
(b)

Figure 1 

Warren-Root model (a) and DFN model (b).

The Warren-Root model assumes that fractures are evenly distributed in the matrix, and the oil and gas in the matrix flow into the fracture and are recovered. In the numerical model, the dual-medium model, including the dual-porosity and dual-permeability model and the dual-porosity and single-permeability model (detailed in Section 2.2), is used to deal with this kind of model. As detailed in the introduction, the DFN model has the following two advantages over the Warren-Root model: (1)More representative of the actual fracture network shape(2)The actual data (such as fracture monitoring, production, and well test data) can be used to simulate reservoir connectivity and continuity more realistically

In this paper, the Fracman software is used to generate a discrete fracture network. In order to describe the fracture better by the mathematical model, Fracman assumes that the fracture is polygonal. It is well known that an ideal fracture in a homogeneous formation should have an elliptical shape. However, since the strata cannot be homogeneous, the shape of the fractures cannot be elliptical. Dershowitz’s results show that the observed cracks are closer to polygons than ovals because of the interweaving of cracks. Both the Veneziano model and Dershowitz model treat cracks as polygons. The Fracman software provides three models for generating discrete fracture networks: Enhanced Baecher model (EB model), Levy-Lee Fractal model (LLF model), and Nearest Neighbor model (NN model). EB model assumes that all fractures are evenly distributed in the reservoir, and the distribution law follows Poisson process, and the distance between the center points of all fractures follows an exponential distribution. Each one may intersect with other cracks, and the Fracman software will modify the cross-cracks according to specific conditions. LLF model is based on the process of “Levitation Flight.” It is a kind of random movement problem in mathematics. The NN model assumes that the fracture density decreases exponentially with the distance from the main fracture. By comparing the three models, Mou Songru found that the calculation points generated by EB model are most widely distributed in the XY plane, and the overall area has the highest coverage, and the length distribution is closer to the power-law distribution. Therefore, it is more reasonable to use EB model to generate cracks. Figure 2 shows the discrete fracture model generated by Fracman (a) and the coarsening results in the numerical model (b).

(a)

(b)

(a)
(b)

Figure 2 

Complex fracture pattern generated by Fracman and coarsening result in reservoir simulation model.

2.2. Dual-Medium Model

The dual-medium model considers the reservoir as matrix pore and fracture system, and the flow laws of the two can be treated separately, and the flow between them can be realized by the flow coefficient. The double-porosity and double-porosity and double-permeability models can be used to extract oil and gas from both fractured and matrix systems. The latter assumes that oil and gas can only be extracted from the fracture system. The typical dual-media model is the Warren-Root model. For a shale gas reservoir, the gas does not flow in the matrix, and the gas enters the fracture to participate in the flow if there is a particular pressure difference between the matrix and the fracture, and the adsorbed gas is desorbed; in this paper, a two-hole single-permeability model is used to simulate shale gas flow.

2.2.1. Adsorption and Desorption of Shale Gas

Adsorption is the phenomenon of gas adhering to the surface of a solid. The fundamental reason is that when the gas and the solid come into contact, they do not meet each other in thermodynamic equilibrium. The adsorption can be divided into physical adsorption and chemical adsorption.

Shale matrix contains a large number of micropores so that the shale matrix inner surface area is very large. When the shale gas molecules come into contact with the shale matrix, both the shale matrix and the shale gas molecules in the outer space attract the shale gas molecules on the solid surface. Therefore, an adsorption field is created on the solid surface, and the shale gas molecules are continuously adsorbed on the solid surface until saturation, or thermal equilibrium, is reached since the shale matrix adsorbs the shale gas molecules by the van der Waals forces. Therefore, shale gas shale adsorption belongs to the physical adsorption. The adsorption of the shale matrix is faster in the initial stage, and the adsorption rate becomes slower and slower until the adsorption equilibrium is reached.

At present, the Langmuir equation is generally used to describe gas adsorption, as shown in Figure 3 and (20). Langmuir volume is the maximum adsorption capacity of shale; Langmuir pressure is the pressure corresponding to half Langmuir volume. The diagram shows that the adsorption gas volume increases gradually with the increase of pressure until the maximum adsorption volume, that is, the Langmuir volume. For example, coal, shale also has an excellent adsorption capacity, especially for the middle-low pressure range of shale reservoirs; its initial capacity is more excellent than conventional sandstone.

Figure 3 

Langmuir isothermal adsorption curve.

Because the adsorption process of shale gas is physical adsorption, the desorption process is basically opposite to the adsorption process, that is, when the pressure is reduced, the gas adsorbed in the pores of shale matrix will be desorbed and return to the pores as free gas. It can be concluded that the desorption process also conforms to the Langmuir equation. This can be based on the decline in shale reservoir pressure to deduce the desorption volume.

Of course, the actual desorption process is not only pressure-related but also shale gas saturation. As shown in Figure 3, when the shale gas content at a given pressure is on a curve, the gas is gradually desorbed from the shale surface as the pressure decreases. And when the shale gas content at a given pressure is below the curve, the gas is not saturated, that is, the amount of adsorbed gas is less than the maximum amount of adsorbed gas under this pressure. Although the pressure drops, the shale gas will not be desorbed immediately but will not be desorbed until the pressure further drops to the reservoir pressure corresponding to the gas content on the isothermal adsorption curve, the shale gas will gradually be sucked out.

2.2.2. Diffusion of Gas

As mentioned earlier, the pore size of shale matrix is very small, and the permeability is very poor, so the typical Duthie’s law cannot accurately describe the migration of shale gas in the matrix; it is generally believed that the migration law of shale gas in matrix pores follows the first diffusion law, Fick’s first law, as shown in formula (24) .

In the formula, is the shale gas volume diffused from shale matrix, m³/d; is diffusion coefficient, m²/d; is the shape factor of shale matrix, m^-2; is the shale matrix volume, m³; is the average gas concentration in shale matrix, m³/m³; and is the gas concentration at the shale fracture interface, m³/m³.

Diffusion is the process by which molecules flow under the driving force of concentration difference. For the development of shale gas, the gas in the fracture is first extracted, which leads to the decrease of the gas concentration in the fracture, while the gas concentration in the matrix pores is higher. Fick’s first law assumes that the concentration at both the fracture and the shale surface is an average, and since the concentration gradient from the matrix pore to the fracture is difficult to describe accurately, this kind of diffusion actually belongs to quasi-stable diffusion.

2.2.3. Double Porosity and Single Permeability Model

It is assumed that in the dual porosity and single permeability model the reservoir includes matrix porosity and fracture, but the gas can only flow from the fracture to the bottom of the well, and the gas flow from the matrix to fracture is characterized by the flow coefficient. Therefore, the pore system and fracture system should be modelled, respectively. (1)Pore system

In the pore system, most of the gas is in the adsorption state. In the matrix pores near the fracture surface, the gas desorption rate is fast and in equilibrium with the free gas, which can be simulated by Langmuir equation, such as formula (19) or formula (20); gas and free gas in matrix pores far from the fracture surface are nonequilibrium. Gas migration in matrix pores is characterized by Fick’s first diffusion law, as in equation (24). (2)Rift system

The flow in the fracture system follows Duthie’s law, which differs from the black oil model in that only gas and water phases are considered, so the differential equation of the fracture system is:

The auxiliary equations are as follows:

In the formula, is the oil-water capillary pressure, MPa. Equations (25)–(28) is a complete gas-water two-phase model. It can be seen from the four equations that , , , and are four unknowns, so the equations are closed.

Assuming that the boundary of the model is closed, the boundary conditions are as follows:

The initial conditions are as follows:

2.3. Finite Difference Discrete Method

This section briefly describes the difference discretization method for block-center grids. Block-center grid is a common grid. The properties of the center point represent the properties of the whole grid. There is an interface between the adjacent grids, and the fluid flows from one grid to another grid. Block-center grid difference method is a simple and mature method, which is widely used in numerical simulation.

Figure 4 shows a block-center grid in the -direction of a one-dimensional reservoir. The -direction is composed of grids, the sum of which is the length of the reservoir in the -direction, representing the point of each grid in the center of the grid block. Figure 5 shows grid in the -direction and the grid adjacent to it. The grid block size (, , ), the grid block boundary (, ), the distance between the center points of different grid blocks (, ), and so on are shown.

Figure 4 

Block-centered network for 1D in the -direction.

Figure 5 

Grid in the -direction and its neighboring grid.

Add a grid block in the -direction to a one-dimensional model, and it becomes a two-dimensional model. As you can imagine, the size of the grid blocks in the 2D model, the boundaries of the grid blocks, and the distance between the centers of the different grid blocks include the sizes in both the - and -directions, such as and , and , and and .

Since the mathematical problem described by the partial differential equation is a continuity problem, it is necessary to discretize it, that is, to treat the -direction continuity problem as a separate problem of . The discrete method is to replace the partial derivative with the difference quotient, that is, to turn the partial differential equation into a difference equation. If the reservoir thickness is small and the dip angle is small, the two-dimensional plane model can be used and ignore the influence of the gravity term, then take equation (25) for example:

The equations (33) are discretized as follows:

In the formula, , , and are the ordinal numbers of the grid in -, - and -directions, respectively; is the time step. Multiply the formula (34) by , and you will get:

Further summarized as follows:

When the IMPES (implicit pressure vivid saturation) method is used to solve the pressure difference from time to time, the pressure difference between time and time is not solved directly, and the results are as follows:

The formula (37) is brought into the formula (36) as shown in formula (38). The gas-phase pressure has been linearized in the left side of equation (38), but for microcompressible fluids, the relative permeability, density, and viscosity of the gas phase are all related to the gas phase pressure or saturation, and these parameters need to be linearized because of their nonlinearity in the equation. Because the absolute permeability is only related to the spatial position but not to the time, the average harmonic method is usually used to deal with the absolute permeability values on the boundary.

For gas-phase relative permeability, gas density and viscosity adopt the “upstream weight” principle, that is, if the fluid flows from the point to the point, it takes the value at the upstream point. If the flow is from six adjacent points, it takes the value at the upstream adjacent point.

From the above, in the equation at the left of (39), all the time-related parameters except the pressure are taken as the value of the last time step. So far, all terms in the equation at the left end of equation (38) have been linearized. The linearization steps of the equation at the right end of the equation (38) are discussed below. Since IMPES transforms the solution of step unknowns into the solution of the difference between step and step , the equation on the right side of (39) also needs to be treated as follows:

The difference for the product of a variable or function can be made as follows:

Using the above rule to treat formula (41), to simplify the omission of , the following are:

In the above equation, is a function of the unknown, and and are the functions of the pressure. The nonlinearity is very weak, so it can be approximately equal to the value of the last time step, namely:

After the above treatment, only and parameters in equation (43) cause nonlinearity. For both gas density and matrix porosity, there is desorption, namely:

The differential of the above formula gives:

Bring in the top type (44) to get:

Combining the equations (38) and (47), the final form of difference dispersion is obtained. This equation contains only two unknowns ( and ). In addition, the water phase equation is just a closed system of equations. The final form of difference dispersion is as follows: