# Model Optimization of Shale Gas Reservoir Volume Fracturing Dissolved Gas Simulation Adsorbed Gas

Authors: Hao Dong; Yi Zhang; Zongwu Li; Chao Jiang; Jiaze Li; Tao Wu; Liting Wang; Chuangjiang Wang; Hao Wang; Fujing Li; Qian Ru
Journal: Geofluids
 (2021)

---

## 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.


---

## Body

## 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 subscriptsw, o, and g denotes the water, oil, and gas phases, respectively, while water, oil, and gas components are, respectively, indicated by i≡W,O,G. The superscript s 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:(1)∇⋅kkroρoμo∇po−γog∇D=∂φρoSo∂t,(2)∇⋅kkroρgdμo∇po−γog∇D+kkrgρgμg∇pg−γo∇D=∂φρgdSo+ρgSg∂t,(3)∇⋅kkrwρwμw∇pw−γw∇D=∂φρwSw∂t.The auxiliary equations are as follows:(4)So+Sg+Sw=1,(5)pw=po−pcow,(6)pg=po−pcog,wherek is the formation permeability, µm2; kro, krg, and krw are the relative permeability of oil, gas, and water; ρo, ρgd, ρg, and ρw are the density of oil, dissolved gas, gas, and water, kg/m3; po, pg, and pw are the pressure of oil, gas, and water phase, MPa; γog, γo, γg, and γw are the weight of oil, oil, gas, and water phase containing dissolved gas, kg⋅m−2⋅s−2, where γog=γo+γgd,γo=ρog, and γw=ρwg; D is the depth from a certain datum, m; φ is the formation porosity; So, Sg, and Sw are the oil, gas, and water saturation, respectively; t is the time, s; pcow is the oil-water capillary pressure, MPa; and pcog 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—po, pg, pw, So, Sg, and Sw. 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:(7)∂po∂tx=0,Lx=0,∂pg∂tx=0,Lx=0,∂pw∂tx=0,Lx=0,(8)∂po∂ty=0,Ly=0,∂pg∂ty=0,Ly=0,∂pw∂ty=0,Ly=0,(9)∂po∂tz=0,Lz=0,∂pg∂tz=0,Lz=0,∂pw∂tz=0,Lz=0,whereLx, Ly, and Lz are the lengths of the x, y, and z directions, respectively, m.The initial conditions are as follows:(10)po,w,gx,y,z,tt=0=pi,(11)Swx,y,z,tt=0=Swi,(12)Sgx,y,z,tt=0=Sgi,wherepi is the original formation pressure, MPa; Swi is the original water saturation, and Sgi 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:(13)Sgmφm=Sgφ,(14)Swmφm=Swφ.In the formula,Sgm and Swm are the gas and water saturation, φm is the model porosity, and φ is the actual porosity. The above two types of the merger are:
(15)Sgm+Swmφm=Sg+Swφ.IfSgm+Swm+Som=1, then:
(16)φm=φ1−Som.In the formula,Som 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:
(17)Sgm=1−SomSg,(18)Swm=1−SomSw.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 Sg and Sw of the shale reservoir should be assigned to the corresponding Sgm and Swm 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 Som=0.45, the real and modified infiltration curves are shown in Table 1.Table 1
Relative permeability curves for reality and modification.
SgSwkrwkrgSwmSgm010.3100.5500.050.950.230.020.52250.02750.10.90.120.060.4950.0550.150.850.050.190.46750.08250.20.80.010.340.440.110.250.75010.41250.1375According to the Langmuir isothermal adsorption equation, the amount of gas adsorbed on a unit volume of shale is:(19)C=VLpp+pL.In the formula,C is the gas adsorption, m3/m3; VL is the Langmuir volume, m3/m3; p is the pressure, MPa; and PL 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:
(20)V=Vmbp1+bp.In the formula,V is the gas adsorption, m3/kg; Vm is the Langmuir constant, m3/kg; and b is the Langmuir pressure constant, MPa-1. In actual shale reservoirs, the amount of adsorbed gas per unit volume is:
(21)ρBVmbp1+bp.In the formula,ρB is the volume density of the shale, kg/m3. In the model, the amount of gas dissolved in the same unit volume can be expressed as:
(22)RsφmSomBo.In the formula,Rs is the dissolved gas-oil ratio, m3/m3, and Bo 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:
(23)Rs=ρBVmBobpφmSom1+bp.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.Table2 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: ρB=1500 kg/m3, Vm=0.0244 m3/kg, Bo=1.2, b=0.3 MPa-1, φm=0.36, and Som=0.45.Table 2
Adsorbed gas volume and dissolved gas-oil ratio in shale.
p, MpaV, m3/tRs, m3/m30.000.000.001.3714.0679.092.7422.33122.514.1127.78149.945.4831.63168.846.8534.51182.668.2236.74193.209.5938.52201.51


#### 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 Figure1.Figure 1
Warren-Root model (a) and DFN model (b).
(a)(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 Section2.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 realisticallyIn 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. Figure2 shows the discrete fracture model generated by Fracman (a) and the coarsening results in the numerical model (b).Figure 2
Complex fracture pattern generated by Fracman and coarsening result in reservoir simulation model.
(a)(b)


### 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 Figure3 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 Figure3, 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) .
(24)qs=DσVmC¯−Cp.In the formula,qs is the shale gas volume diffused from shale matrix, m3/d; Dis diffusion coefficient, m2/d; σ is the shape factor of shale matrix, m-2; Vm is the shale matrix volume, m3; C− is the average gas concentration in shale matrix, m3/m3; and Cp is the gas concentration at the shale fracture interface, m3/m3.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 systemIn 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 systemThe 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:(25)∇⋅kkrgρgμg∇pg−γg∇D=∂φρgSg∂t,(26)∇⋅kkrwρwμw∇pw−γw∇D=∂φρwSw∂t.The auxiliary equations are as follows:(27)Sg+Sw=1,(28)pw=pg−pcgw.In the formula,pcgw 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 pg, pw, Sg, and Sw are four unknowns, so the equations are closed.Assuming that the boundary of the model is closed, the boundary conditions are as follows:(29)∂pg∂tx=0,Lx=0,∂pw∂tx=0,Lx=0,∂pg∂ty=0,Ly=0,∂pw∂ty=0,Ly=0,∂pg∂tz=0,Lz=0,∂pw∂tz=0,Lz=0.The initial conditions are as follows:(30)pw,gx,y,z,tt=0=pi,(31)Swx,y,z,tt=0=Swi,(32)Sgx,y,z,tt=0=Sgi.


### 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.Figure4 shows a block-center grid in the x-direction of a one-dimensional reservoir. The x-direction is composed of Nx grids, the sum of which is the length of the reservoir in the x-direction, representing the point of each grid in the center of the grid block. Figure 5 shows grid I in the x-direction and the grid adjacent to it. The grid block size (Δxi−1, Δxi, Δxi+1), the grid block boundary (xi−1/2, xi+1/2), the distance between the center points of different grid blocks (Δxi−1/2, Δxi+1/2), and so on are shown.Figure 4
Block-centered network for 1D in thex-direction.Figure 5
Gridi in the x-direction and its neighboring grid.Add a grid block in they-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 x- and y-directions, such as Δxi and Δyj, xi−1/2 and yj−1/2, and Δxi−1/2 and Δyj−1/2.Since the mathematical problem described by the partial differential equation is a continuity problem, it is necessary to discretize it, that is, to treat thex-direction continuity problem as a separate problem of Nx. 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:
(33)∇⋅kkrgρgμg∇pg=∂φρgSg∂t.The equations (33) are discretized as follows:
(34)1Δxikkrgρgμgi+1/2pgi+1n+1−pgin+1Δxi+1/2+kkrgρgμgi−1/2pgi−1n+1−pgin+1Δxi−1/2+1Δyjkkrgρgμgj+1/2pgj+1n+1−pgjn+1Δyj+1/2+kkrgρgμgj−1/2pgj−1n+1−pgjn+1Δyj−1/2+1Δzkkkrgρgμgk+1/2pgk+1n+1−pgkn+1Δzk+1/2+kkrgρgμgk−1/2pgk−1n+1−pgkn+1Δzk−1/2=φρgSgn+1−φρgSgnΔt.In the formula,i, j, and k are the ordinal numbers of the grid in x-, y- and z-directions, respectively; n is the n time step. Multiply the formula (34) by ΔxiΔyjΔzk, and you will get:
(35)ΔyjΔzkkkrgρgμgi+1/2pgi+1n+1−pgin+1Δxi+1/2+kkrgρgμgi−1/2pgi−1n+1−pgin+1Δxi−1/2+ΔxiΔzkkkrgρgμgj+1/2pgj+1n+1−pgjn+1Δyj+1/2+kkrgρgμgj−1/2pgj‐1n+1−pgjn+1Δyj−1/2+ΔxiΔyjkkrgρgμgk+1/2pgk+1n+1−pgkn+1Δzk+1/2+kkrgρgμgk−1/2pgk−1n+1−pgkn+1Δzk−1/2=ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt.Further summarized as follows:(36)ΔyjΔzkΔxi+1/2kkrgρgμgi+1/2pgi+1n+1−pgin+1+ΔyjΔzkΔxi−1/2kkrgρgμgi−1/2pgi−1n+1−pgin+1+ΔxiΔzkΔyj+1/2kkrgρgμgj+1/2pgj+1n+1−pgjn+1+ΔxiΔzkΔyj−1/2kkrgρgμgj−1/2pgj−1n+1−pgjn+1+ΔxiΔyjΔzk+1/2kkrgρgμgk+1/2pgk+1n+1−pgkn+1+ΔxiΔyjΔzk−1/2kkrgρgμgk−1/2pgk−1n+1−pgkn+1=ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt.When the IMPES (implicit pressure vivid saturation) method is used to solve the pressure difference fromn+1 time to n+1 time, the pressure difference between n time and n+1 time is not solved directly, and the results are as follows:
(37)δp=pn+1−pn.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 k 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.
(38)ΔyjΔzkΔxi+1/2kkrgρgμgi+1/2pgi+1n−pgin+δpgi+1−δpgi+ΔyjΔzkΔxi−1/2kkrgρgμgi−1/2pgi−1n−pgin+δpgi−1−δpgi+ΔxiΔzkΔyj+1/2kkrgρgμgj+1/2pgj+1n−pgjn+δpgj+1−δpgj+ΔxiΔzkΔyj−1/2kkrgρgμgj−1/2pgj−1n−pgjn+δpgj−1−δpgj+ΔxiΔyjΔzk+1/2kkrgρgμgk+1/2pgk+1n−pgkn+δpgk+1−δpgk+ΔxiΔyjΔzk−1/2kkrgρgμgk−1/2pgk−1n−pgkn+δpgk−1−δpgk=ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt,(39)ki+1/2=2kiki+1ki+ki+1.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.(40)krg,ρg,μgi+1/2=krg,ρg,μgi,formitoi+1,krg,ρg,μgi+1,fromi+1toi.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 n+1 unknowns into the solution of the difference between step n+1 and step n, the equation on the right side of (39) also needs to be treated as follows:
(41)ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt=ΔxiΔyjΔzkδφρgSgΔt.The difference for the product of a variable or function can be made as follows:(42)δab=an+1δb+bnδa.Using the above rule to treat formula (41), to simplify the omission of ΔxiΔyjΔzk, the following are:
(43)δφρgSgΔt=φn+1ρgn+1δSg+Sgnφnδρg+ρgn+1δφΔt.In the above equation,δSg is a function of the unknown, and φn+1 and ρn+1 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:
(44)φn+1=φn,ρgn+1=ρgn.After the above treatment, onlyδρg and δφ parameters in equation (43) cause nonlinearity. For both gas density and matrix porosity, there is desorption, namely:
(45)φ=φ01+Cφpg−pg0,ρg=ρg01+Cppg−pg0.The differential of the above formula gives:(46)δφ=φ0Cφδpg,δρg=ρg0Cpδpg.Bring in the top type (44) to get:
(47)δφρgSgΔt=φnρgnδSg+Sgnφnρg0Cρδpg+ρgnφ0CφδpgΔt.Combining the equations (38) and (47), the final form of difference dispersion is obtained. This equation contains only two unknowns (δSg and δp). In addition, the water phase equation is just a closed system of equations. The final form of difference dispersion is as follows:
(48)ΔyjΔzkΔxi+1/2kkrgρgμgi+1/2pgi+1n−pgin+δpgi+1−δpgi+ΔyjΔzkΔxi−1/2kkrgρgμgi−1/2pgi−1n−pgin+δpgi−1−δpgi+ΔxiΔzkΔyj+1/2kkrgρgμgj+1/2pgj+1n−pgjn+δpgj+1−δpgj+ΔxiΔzkΔyj−1/2kkrgρgμgj−1/2pgj−1n−pgjn+δpgj−1−δpgj+ΔxiΔyjΔzk+1/2kkrgρgμgk+1/2pgk+1n−pgkn+δpgk+1−δpgk+ΔxiΔyjΔzk−1/2kkrgρgμgk−1/2pgk−1n−pgkn+δpgk−1−δpgk=ΔxiΔyjΔzkφnρgnδSg+Sgnφnρg0Cρδpg+ρgnφ0CφδpgΔt.In the numerical simulation of shale gas, since the setting of natural fractures and artificial fractures are involved, the setting of the two often requires the method of grid densification, so the method of grid encryption is critical, and reasonable grid encryption can improve the simulation accuracy and reduce the possibility of nonconvergence.Local mesh encryption is the use of a dense mesh only in the most desired part of the solution area and a coarse mesh elsewhere. That is, the encrypted grid lines do not necessarily extend to the boundary of the solution area but can be distributed only where encryption is required, as shown in Figure6. Local grid encryption can subdivide a grid without affecting the shape of its neighbors. Due to the variabilities of the contact relationship between grids in space and time in local grid encryption, it is necessary to implement this method in the program, such as discretization, mesh generation, and mesh numbering; in addition, a convenient and flexible data management system is also needed to solve the coefficient matrix.Figure 6
Grid refinement in Cartesian coordinate.In conventional grid systems, the fluid exchange between grids is relatively simple because no one side of the grid is adjacent to only one grid. But for the local grid encryption system, because some of the adjacent grids may have more than one on one side, it is necessary to consider the particularity when considering the discretization of traffic between adjacent grids. If there is a locally encrypted grid of contacts as shown in Figure7, assuming that the absolute permeability, relative permeability, and fluid viscosity of the grids remain constant, then the sum of the flows Q which is from grid 2 to grid 6 to grid 1 can be defined as:
(49)Q=kkrΔzμΔxp2−p1Δy2+p3−p1Δy3+⋯+p6−p1Δy6.Figure 7
Grid refinement in Cartesian coordinate.In the formula,Δx is the distance from grid 2 to grid 6 from the center of the grid to grid 1; Δyi is the length of grid i in the y-direction; Δz is the thickness of the reservoir; and pi is the pressure of grid i, i=1,2,3,⋯,6.In the formula (49) knows the flow rate between the grids, it is possible to derive a set of differential equation suitable for local grid encryption. The remaining steps are not covered. This method is reliable and straightforward, and the error is within the engineering calculation range.


## 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:(1)∇⋅kkroρoμo∇po−γog∇D=∂φρoSo∂t,(2)∇⋅kkroρgdμo∇po−γog∇D+kkrgρgμg∇pg−γo∇D=∂φρgdSo+ρgSg∂t,(3)∇⋅kkrwρwμw∇pw−γw∇D=∂φρwSw∂t.The auxiliary equations are as follows:(4)So+Sg+Sw=1,(5)pw=po−pcow,(6)pg=po−pcog,wherek is the formation permeability, µm2; kro, krg, and krw are the relative permeability of oil, gas, and water; ρo, ρgd, ρg, and ρw are the density of oil, dissolved gas, gas, and water, kg/m3; po, pg, and pw are the pressure of oil, gas, and water phase, MPa; γog, γo, γg, and γw are the weight of oil, oil, gas, and water phase containing dissolved gas, kg⋅m−2⋅s−2, where γog=γo+γgd,γo=ρog, and γw=ρwg; D is the depth from a certain datum, m; φ is the formation porosity; So, Sg, and Sw are the oil, gas, and water saturation, respectively; t is the time, s; pcow is the oil-water capillary pressure, MPa; and pcog 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—po, pg, pw, So, Sg, and Sw. 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:(7)∂po∂tx=0,Lx=0,∂pg∂tx=0,Lx=0,∂pw∂tx=0,Lx=0,(8)∂po∂ty=0,Ly=0,∂pg∂ty=0,Ly=0,∂pw∂ty=0,Ly=0,(9)∂po∂tz=0,Lz=0,∂pg∂tz=0,Lz=0,∂pw∂tz=0,Lz=0,whereLx, Ly, and Lz are the lengths of the x, y, and z directions, respectively, m.The initial conditions are as follows:(10)po,w,gx,y,z,tt=0=pi,(11)Swx,y,z,tt=0=Swi,(12)Sgx,y,z,tt=0=Sgi,wherepi is the original formation pressure, MPa; Swi is the original water saturation, and Sgi 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:(13)Sgmφm=Sgφ,(14)Swmφm=Swφ.In the formula,Sgm and Swm are the gas and water saturation, φm is the model porosity, and φ is the actual porosity. The above two types of the merger are:
(15)Sgm+Swmφm=Sg+Swφ.IfSgm+Swm+Som=1, then:
(16)φm=φ1−Som.In the formula,Som 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:
(17)Sgm=1−SomSg,(18)Swm=1−SomSw.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 Sg and Sw of the shale reservoir should be assigned to the corresponding Sgm and Swm 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 Som=0.45, the real and modified infiltration curves are shown in Table 1.Table 1
Relative permeability curves for reality and modification.
SgSwkrwkrgSwmSgm010.3100.5500.050.950.230.020.52250.02750.10.90.120.060.4950.0550.150.850.050.190.46750.08250.20.80.010.340.440.110.250.75010.41250.1375According to the Langmuir isothermal adsorption equation, the amount of gas adsorbed on a unit volume of shale is:(19)C=VLpp+pL.In the formula,C is the gas adsorption, m3/m3; VL is the Langmuir volume, m3/m3; p is the pressure, MPa; and PL 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:
(20)V=Vmbp1+bp.In the formula,V is the gas adsorption, m3/kg; Vm is the Langmuir constant, m3/kg; and b is the Langmuir pressure constant, MPa-1. In actual shale reservoirs, the amount of adsorbed gas per unit volume is:
(21)ρBVmbp1+bp.In the formula,ρB is the volume density of the shale, kg/m3. In the model, the amount of gas dissolved in the same unit volume can be expressed as:
(22)RsφmSomBo.In the formula,Rs is the dissolved gas-oil ratio, m3/m3, and Bo 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:
(23)Rs=ρBVmBobpφmSom1+bp.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.Table2 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: ρB=1500 kg/m3, Vm=0.0244 m3/kg, Bo=1.2, b=0.3 MPa-1, φm=0.36, and Som=0.45.Table 2
Adsorbed gas volume and dissolved gas-oil ratio in shale.
p, MpaV, m3/tRs, m3/m30.000.000.001.3714.0679.092.7422.33122.514.1127.78149.945.4831.63168.846.8534.51182.668.2236.74193.209.5938.52201.51


### 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 Figure1.Figure 1
Warren-Root model (a) and DFN model (b).
(a)(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 Section2.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 realisticallyIn 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. Figure2 shows the discrete fracture model generated by Fracman (a) and the coarsening results in the numerical model (b).Figure 2
Complex fracture pattern generated by Fracman and coarsening result in reservoir simulation model.
(a)(b)


## 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:(1)∇⋅kkroρoμo∇po−γog∇D=∂φρoSo∂t,(2)∇⋅kkroρgdμo∇po−γog∇D+kkrgρgμg∇pg−γo∇D=∂φρgdSo+ρgSg∂t,(3)∇⋅kkrwρwμw∇pw−γw∇D=∂φρwSw∂t.The auxiliary equations are as follows:(4)So+Sg+Sw=1,(5)pw=po−pcow,(6)pg=po−pcog,wherek is the formation permeability, µm2; kro, krg, and krw are the relative permeability of oil, gas, and water; ρo, ρgd, ρg, and ρw are the density of oil, dissolved gas, gas, and water, kg/m3; po, pg, and pw are the pressure of oil, gas, and water phase, MPa; γog, γo, γg, and γw are the weight of oil, oil, gas, and water phase containing dissolved gas, kg⋅m−2⋅s−2, where γog=γo+γgd,γo=ρog, and γw=ρwg; D is the depth from a certain datum, m; φ is the formation porosity; So, Sg, and Sw are the oil, gas, and water saturation, respectively; t is the time, s; pcow is the oil-water capillary pressure, MPa; and pcog 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—po, pg, pw, So, Sg, and Sw. 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:(7)∂po∂tx=0,Lx=0,∂pg∂tx=0,Lx=0,∂pw∂tx=0,Lx=0,(8)∂po∂ty=0,Ly=0,∂pg∂ty=0,Ly=0,∂pw∂ty=0,Ly=0,(9)∂po∂tz=0,Lz=0,∂pg∂tz=0,Lz=0,∂pw∂tz=0,Lz=0,whereLx, Ly, and Lz are the lengths of the x, y, and z directions, respectively, m.The initial conditions are as follows:(10)po,w,gx,y,z,tt=0=pi,(11)Swx,y,z,tt=0=Swi,(12)Sgx,y,z,tt=0=Sgi,wherepi is the original formation pressure, MPa; Swi is the original water saturation, and Sgi 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:(13)Sgmφm=Sgφ,(14)Swmφm=Swφ.In the formula,Sgm and Swm are the gas and water saturation, φm is the model porosity, and φ is the actual porosity. The above two types of the merger are:
(15)Sgm+Swmφm=Sg+Swφ.IfSgm+Swm+Som=1, then:
(16)φm=φ1−Som.In the formula,Som 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:
(17)Sgm=1−SomSg,(18)Swm=1−SomSw.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 Sg and Sw of the shale reservoir should be assigned to the corresponding Sgm and Swm 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 Som=0.45, the real and modified infiltration curves are shown in Table 1.Table 1
Relative permeability curves for reality and modification.
SgSwkrwkrgSwmSgm010.3100.5500.050.950.230.020.52250.02750.10.90.120.060.4950.0550.150.850.050.190.46750.08250.20.80.010.340.440.110.250.75010.41250.1375According to the Langmuir isothermal adsorption equation, the amount of gas adsorbed on a unit volume of shale is:(19)C=VLpp+pL.In the formula,C is the gas adsorption, m3/m3; VL is the Langmuir volume, m3/m3; p is the pressure, MPa; and PL 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:
(20)V=Vmbp1+bp.In the formula,V is the gas adsorption, m3/kg; Vm is the Langmuir constant, m3/kg; and b is the Langmuir pressure constant, MPa-1. In actual shale reservoirs, the amount of adsorbed gas per unit volume is:
(21)ρBVmbp1+bp.In the formula,ρB is the volume density of the shale, kg/m3. In the model, the amount of gas dissolved in the same unit volume can be expressed as:
(22)RsφmSomBo.In the formula,Rs is the dissolved gas-oil ratio, m3/m3, and Bo 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:
(23)Rs=ρBVmBobpφmSom1+bp.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.Table2 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: ρB=1500 kg/m3, Vm=0.0244 m3/kg, Bo=1.2, b=0.3 MPa-1, φm=0.36, and Som=0.45.Table 2
Adsorbed gas volume and dissolved gas-oil ratio in shale.
p, MpaV, m3/tRs, m3/m30.000.000.001.3714.0679.092.7422.33122.514.1127.78149.945.4831.63168.846.8534.51182.668.2236.74193.209.5938.52201.51


## 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 Figure1.Figure 1
Warren-Root model (a) and DFN model (b).
(a)(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 Section2.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 realisticallyIn 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. Figure2 shows the discrete fracture model generated by Fracman (a) and the coarsening results in the numerical model (b).Figure 2
Complex fracture pattern generated by Fracman and coarsening result in reservoir simulation model.
(a)(b)


## 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 Figure3 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 Figure3, 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) .
(24)qs=DσVmC¯−Cp.In the formula,qs is the shale gas volume diffused from shale matrix, m3/d; Dis diffusion coefficient, m2/d; σ is the shape factor of shale matrix, m-2; Vm is the shale matrix volume, m3; C− is the average gas concentration in shale matrix, m3/m3; and Cp is the gas concentration at the shale fracture interface, m3/m3.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 systemIn 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 systemThe 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:(25)∇⋅kkrgρgμg∇pg−γg∇D=∂φρgSg∂t,(26)∇⋅kkrwρwμw∇pw−γw∇D=∂φρwSw∂t.The auxiliary equations are as follows:(27)Sg+Sw=1,(28)pw=pg−pcgw.In the formula,pcgw 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 pg, pw, Sg, and Sw are four unknowns, so the equations are closed.Assuming that the boundary of the model is closed, the boundary conditions are as follows:(29)∂pg∂tx=0,Lx=0,∂pw∂tx=0,Lx=0,∂pg∂ty=0,Ly=0,∂pw∂ty=0,Ly=0,∂pg∂tz=0,Lz=0,∂pw∂tz=0,Lz=0.The initial conditions are as follows:(30)pw,gx,y,z,tt=0=pi,(31)Swx,y,z,tt=0=Swi,(32)Sgx,y,z,tt=0=Sgi.


## 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 Figure3 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 Figure3, 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) .
(24)qs=DσVmC¯−Cp.In the formula,qs is the shale gas volume diffused from shale matrix, m3/d; Dis diffusion coefficient, m2/d; σ is the shape factor of shale matrix, m-2; Vm is the shale matrix volume, m3; C− is the average gas concentration in shale matrix, m3/m3; and Cp is the gas concentration at the shale fracture interface, m3/m3.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 systemIn 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 systemThe 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:(25)∇⋅kkrgρgμg∇pg−γg∇D=∂φρgSg∂t,(26)∇⋅kkrwρwμw∇pw−γw∇D=∂φρwSw∂t.The auxiliary equations are as follows:(27)Sg+Sw=1,(28)pw=pg−pcgw.In the formula,pcgw 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 pg, pw, Sg, and Sw are four unknowns, so the equations are closed.Assuming that the boundary of the model is closed, the boundary conditions are as follows:(29)∂pg∂tx=0,Lx=0,∂pw∂tx=0,Lx=0,∂pg∂ty=0,Ly=0,∂pw∂ty=0,Ly=0,∂pg∂tz=0,Lz=0,∂pw∂tz=0,Lz=0.The initial conditions are as follows:(30)pw,gx,y,z,tt=0=pi,(31)Swx,y,z,tt=0=Swi,(32)Sgx,y,z,tt=0=Sgi.


## 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.Figure4 shows a block-center grid in the x-direction of a one-dimensional reservoir. The x-direction is composed of Nx grids, the sum of which is the length of the reservoir in the x-direction, representing the point of each grid in the center of the grid block. Figure 5 shows grid I in the x-direction and the grid adjacent to it. The grid block size (Δxi−1, Δxi, Δxi+1), the grid block boundary (xi−1/2, xi+1/2), the distance between the center points of different grid blocks (Δxi−1/2, Δxi+1/2), and so on are shown.Figure 4
Block-centered network for 1D in thex-direction.Figure 5
Gridi in the x-direction and its neighboring grid.Add a grid block in they-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 x- and y-directions, such as Δxi and Δyj, xi−1/2 and yj−1/2, and Δxi−1/2 and Δyj−1/2.Since the mathematical problem described by the partial differential equation is a continuity problem, it is necessary to discretize it, that is, to treat thex-direction continuity problem as a separate problem of Nx. 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:
(33)∇⋅kkrgρgμg∇pg=∂φρgSg∂t.The equations (33) are discretized as follows:
(34)1Δxikkrgρgμgi+1/2pgi+1n+1−pgin+1Δxi+1/2+kkrgρgμgi−1/2pgi−1n+1−pgin+1Δxi−1/2+1Δyjkkrgρgμgj+1/2pgj+1n+1−pgjn+1Δyj+1/2+kkrgρgμgj−1/2pgj−1n+1−pgjn+1Δyj−1/2+1Δzkkkrgρgμgk+1/2pgk+1n+1−pgkn+1Δzk+1/2+kkrgρgμgk−1/2pgk−1n+1−pgkn+1Δzk−1/2=φρgSgn+1−φρgSgnΔt.In the formula,i, j, and k are the ordinal numbers of the grid in x-, y- and z-directions, respectively; n is the n time step. Multiply the formula (34) by ΔxiΔyjΔzk, and you will get:
(35)ΔyjΔzkkkrgρgμgi+1/2pgi+1n+1−pgin+1Δxi+1/2+kkrgρgμgi−1/2pgi−1n+1−pgin+1Δxi−1/2+ΔxiΔzkkkrgρgμgj+1/2pgj+1n+1−pgjn+1Δyj+1/2+kkrgρgμgj−1/2pgj‐1n+1−pgjn+1Δyj−1/2+ΔxiΔyjkkrgρgμgk+1/2pgk+1n+1−pgkn+1Δzk+1/2+kkrgρgμgk−1/2pgk−1n+1−pgkn+1Δzk−1/2=ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt.Further summarized as follows:(36)ΔyjΔzkΔxi+1/2kkrgρgμgi+1/2pgi+1n+1−pgin+1+ΔyjΔzkΔxi−1/2kkrgρgμgi−1/2pgi−1n+1−pgin+1+ΔxiΔzkΔyj+1/2kkrgρgμgj+1/2pgj+1n+1−pgjn+1+ΔxiΔzkΔyj−1/2kkrgρgμgj−1/2pgj−1n+1−pgjn+1+ΔxiΔyjΔzk+1/2kkrgρgμgk+1/2pgk+1n+1−pgkn+1+ΔxiΔyjΔzk−1/2kkrgρgμgk−1/2pgk−1n+1−pgkn+1=ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt.When the IMPES (implicit pressure vivid saturation) method is used to solve the pressure difference fromn+1 time to n+1 time, the pressure difference between n time and n+1 time is not solved directly, and the results are as follows:
(37)δp=pn+1−pn.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 k 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.
(38)ΔyjΔzkΔxi+1/2kkrgρgμgi+1/2pgi+1n−pgin+δpgi+1−δpgi+ΔyjΔzkΔxi−1/2kkrgρgμgi−1/2pgi−1n−pgin+δpgi−1−δpgi+ΔxiΔzkΔyj+1/2kkrgρgμgj+1/2pgj+1n−pgjn+δpgj+1−δpgj+ΔxiΔzkΔyj−1/2kkrgρgμgj−1/2pgj−1n−pgjn+δpgj−1−δpgj+ΔxiΔyjΔzk+1/2kkrgρgμgk+1/2pgk+1n−pgkn+δpgk+1−δpgk+ΔxiΔyjΔzk−1/2kkrgρgμgk−1/2pgk−1n−pgkn+δpgk−1−δpgk=ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt,(39)ki+1/2=2kiki+1ki+ki+1.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.(40)krg,ρg,μgi+1/2=krg,ρg,μgi,formitoi+1,krg,ρg,μgi+1,fromi+1toi.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 n+1 unknowns into the solution of the difference between step n+1 and step n, the equation on the right side of (39) also needs to be treated as follows:
(41)ΔxiΔyjΔzkφρgSgn+1−φρgSgnΔt=ΔxiΔyjΔzkδφρgSgΔt.The difference for the product of a variable or function can be made as follows:(42)δab=an+1δb+bnδa.Using the above rule to treat formula (41), to simplify the omission of ΔxiΔyjΔzk, the following are:
(43)δφρgSgΔt=φn+1ρgn+1δSg+Sgnφnδρg+ρgn+1δφΔt.In the above equation,δSg is a function of the unknown, and φn+1 and ρn+1 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:
(44)φn+1=φn,ρgn+1=ρgn.After the above treatment, onlyδρg and δφ parameters in equation (43) cause nonlinearity. For both gas density and matrix porosity, there is desorption, namely:
(45)φ=φ01+Cφpg−pg0,ρg=ρg01+Cppg−pg0.The differential of the above formula gives:(46)δφ=φ0Cφδpg,δρg=ρg0Cpδpg.Bring in the top type (44) to get:
(47)δφρgSgΔt=φnρgnδSg+Sgnφnρg0Cρδpg+ρgnφ0CφδpgΔt.Combining the equations (38) and (47), the final form of difference dispersion is obtained. This equation contains only two unknowns (δSg and δp). In addition, the water phase equation is just a closed system of equations. The final form of difference dispersion is as follows:
(48)ΔyjΔzkΔxi+1/2kkrgρgμgi+1/2pgi+1n−pgin+δpgi+1−δpgi+ΔyjΔzkΔxi−1/2kkrgρgμgi−1/2pgi−1n−pgin+δpgi−1−δpgi+ΔxiΔzkΔyj+1/2kkrgρgμgj+1/2pgj+1n−pgjn+δpgj+1−δpgj+ΔxiΔzkΔyj−1/2kkrgρgμgj−1/2pgj−1n−pgjn+δpgj−1−δpgj+ΔxiΔyjΔzk+1/2kkrgρgμgk+1/2pgk+1n−pgkn+δpgk+1−δpgk+ΔxiΔyjΔzk−1/2kkrgρgμgk−1/2pgk−1n−pgkn+δpgk−1−δpgk=ΔxiΔyjΔzkφnρgnδSg+Sgnφnρg0Cρδpg+ρgnφ0CφδpgΔt.In the numerical simulation of shale gas, since the setting of natural fractures and artificial fractures are involved, the setting of the two often requires the method of grid densification, so the method of grid encryption is critical, and reasonable grid encryption can improve the simulation accuracy and reduce the possibility of nonconvergence.Local mesh encryption is the use of a dense mesh only in the most desired part of the solution area and a coarse mesh elsewhere. That is, the encrypted grid lines do not necessarily extend to the boundary of the solution area but can be distributed only where encryption is required, as shown in Figure6. Local grid encryption can subdivide a grid without affecting the shape of its neighbors. Due to the variabilities of the contact relationship between grids in space and time in local grid encryption, it is necessary to implement this method in the program, such as discretization, mesh generation, and mesh numbering; in addition, a convenient and flexible data management system is also needed to solve the coefficient matrix.Figure 6
Grid refinement in Cartesian coordinate.In conventional grid systems, the fluid exchange between grids is relatively simple because no one side of the grid is adjacent to only one grid. But for the local grid encryption system, because some of the adjacent grids may have more than one on one side, it is necessary to consider the particularity when considering the discretization of traffic between adjacent grids. If there is a locally encrypted grid of contacts as shown in Figure7, assuming that the absolute permeability, relative permeability, and fluid viscosity of the grids remain constant, then the sum of the flows Q which is from grid 2 to grid 6 to grid 1 can be defined as:
(49)Q=kkrΔzμΔxp2−p1Δy2+p3−p1Δy3+⋯+p6−p1Δy6.Figure 7
Grid refinement in Cartesian coordinate.In the formula,Δx is the distance from grid 2 to grid 6 from the center of the grid to grid 1; Δyi is the length of grid i in the y-direction; Δz is the thickness of the reservoir; and pi is the pressure of grid i, i=1,2,3,⋯,6.In the formula (49) knows the flow rate between the grids, it is possible to derive a set of differential equation suitable for local grid encryption. The remaining steps are not covered. This method is reliable and straightforward, and the error is within the engineering calculation range.


## 3. Contrast and Analysis-Comparison between Single-Medium Model and Double-Medium Model

The conceptual models are established for the single-medium model and the dual-medium model, respectively, in which the adsorbed gas in the single medium is represented by the dissolved gas in the oil phase, and the natural fracture system in the dual medium is represented by the flow coefficient. Both models adopt the grid division mode of80×49×5; each grid is square in the plane, the side length is 25 m, the vertical height of each grid is 10 m, the horizontal well length is 1000 m, the artificial fracture is 6, and the fracture length is 200 m; the artificial fracture and natural fracture are simulated by local mesh encryption method.The simulated formation depth is 3000 m, the formation pressure is 30 MPa, the average permeability is0.001×10−3μm2, the average porosity is 0.05, the surface density of crude oil is 850 kg/m3, the compressibility of water is 1.0×10−4MPa−1, the viscosity of water is 0.2 mPa·S, and the compressibility of rock is 1.0×10−5MPa−1. The curves of volume coefficient and viscosity of free gas with pressure are shown in Figure 8. The curves of dissolved gas-oil ratio and volume coefficient with bubble point pressure are shown in Figure 9, and the curves of relative permeability of gas to water are shown in Figure 10. If Som=0.4 and Bo=1.2, the adsorption gas volume and the corresponding gas-water relative permeability curves can be calculated by the method described in Section 2.1.2, as shown in Figures 11 and 12, respectively; the results are shown in Tables 3 and 4.Figure 8
Volume factor and viscosity changing with pressure for free gas.Figure 9
Dissolved gas-oil ratio and volume factor changing with bubble pressure.Figure 10
Gas-water relative permeability curve in single medium.Figure 11
Langmuir isothermal adsorption curve.Figure 12
Gas-water relative permeability curve in dual medium.Table 3
Relative permeability curves for reality and modification.
SgSwkrwkrgSwmSgm0.001.00100.600.120.880.5690.000890.5290.0710.170.840.4450.00130.5010.0990.220.780.320.00240.4660.1340.280.720.230.00380.4310.1690.350.650.1530.00650.3880.2120.400.600.1210.00980.360.240.470.530.0870.0230.3160.2840.560.440.0560.0480.2660.3340.630.370.0240.0880.2220.3780.670.330.0120.1240.20.4Table 4
Adsorbed gas volume and dissolved gas-oil ratio in shale.
p, MPaV, m3/m3Rs, m3/m32.760.4929.395.520.9959.678.271.4889.0511.031.97118.4413.792.46147.4716.552.92175.4419.313.35201.2622.063.77226.2024.824.13247.5727.584.45267.1630.344.75284.9733.094.7535.854.7538.614.75According to the above parameters, the single medium and the double medium are simulated, respectively. As shown in Figure13, the pressure curves calculated by the two models are changed with time. As can be seen from the diagram, the pressure changes simulated by the two models are similar. For a single medium, the natural fracture generated by Fracman is orthogonal to the artificial fracture, so it can be seen from the pressure change chart of one-month production that the pressure propagation is distributed in a network with the propagation of the artificial fracture and the natural fracture. For dual media, because the simulation of natural fractures is realized by the flow of fracture system and matrix system, it can be seen from the pressure change chart of 1-month production that the pressure propagation only extends along with the artificial fractures; the propagation of pressure in natural fractures cannot be shown graphically. With the development of production, the pressure deficit gradually appears around the horizontal wellbore and the artificial fracture, which makes the pressure continuously spread outwards. After 1 year of production, the pressure drop around the artificial fracture is severe. For single medium, because the oil phase does not flow, with the pressure drop, the gas is continuously released from the oil phase, and for double medium, after the earlier period of free gas extraction, with the pressure drop, the adsorbed gas is continuously desorbed, become the main factor that assures yield. It can also be seen that for conventional low permeability gas reservoirs, the area of oil leakage around the fracture after artificial fracturing should be elliptical, while for shale gas reservoirs, the area of oil leakage should be closer to the rectangle (blue area). This is mainly because the permeability of a shale gas reservoir is very low, the percolation conditions are deplorable, and the pressure expands very slowly. When the pressure difference between two points is less than the resistance required for the gas to flow, the pressure stops expanding forward. Because the area of oil leaking from each fracture is very small, the areas of multiple fractures add up to be more rectangular.Figure 13
Pressure distribution for single medium and dual medium.Figures14 and 15 show a comparison of the daily and cumulative gas yields of the two models. As can be seen from the diagram, the gas production in 3 years simulated by dual media is slightly higher than that in single media, which is mainly related to the way the natural fractures are treated by the two models. For the single medium, except for the artificial fracture, which is the main seepage passage, it is limited by the computer memory, and the natural fracture is mainly disposed around the wellbore. However, the natural fractures in the dual medium are distributed in the whole simulation area, so when the pressure expands beyond the natural fractures, the gas in the single medium can hardly flow, and the dual medium is not affected by this factor. Taken together, however, both methods are reasonable for simulating shale gas production.Figure 14
Comparison of gas production rate for a single medium and dual medium.Figure 15
Comparison of cumulative gas production for a single medium and dual medium.


## 4. Concluding Remarks

In this paper, the black oil model is used to simulate the shale gas reservoir, and the dissolved gas is used to replace the adsorbed gas. Although the treatment of natural fractures is different, the simulation results of the two models are very close, which proves that the method of replacing adsorbed gas with dissolved gas is feasible.Based on the black oil model of a single medium and the principle of equal total gas content, the adsorption gas in shale gas reservoir is replaced by dissolved gas of formation crude oil, and the phase-permeability curve is transformed by changing the porosity and oil-gas-water saturation of the model. The relation curve between the dissolved gas-oil ratio and the pressure is found by the relation between the dissolved gas-oil ratio and the adsorbed gas quantity, as well as the natural fracture network is generated by the discrete fracture model provided by Fracman. The percolation characteristics of shale gas are simulated by adding Langmuir isothermal adsorption curve and gas diffusion principle. Comparing the results of the two models, it is found that the results of the two models are close to each other, and the black oil model can be used to simulate the seepage process of shale gas.



---

*Source: 6631994-2021-02-19.xml*