Shale reservoir has a complex pore structure, and gas seepage in the reservoir has multiple flow mechanisms such as Knudsen diffusion, surface adsorption diffusion, slip flow and viscous flow (Geng et al., 2016; Wu et al., 2016; Zhang et al., 2018) (as shown in Figure 2).

In the process of shale gas development, the thickness of adsorption layer decreases, and the effective flow pore size increases further. The effective pore radius of removing the adsorption layer can be expressed as: r a b s = r − d M θ = r − d M p / Z p / Z + p L (1) Where, d _M is the molecular diameter of the gas, m; r is the actual radius of the matrix pores, m; θ is the true gas coverage, dimensionless; p is the current pore pressure, Pa; p _L is Langmuir pressure, Pa; Z is the true gas compression factor, dimensionless.

Considering the stress-sensitive effect, the porosity and permeability of shale matrix related to pore pressure can be expressed as: φ p = φ r + φ i − φ r e − η p i − p k p = k i e − Ψ p i − p (2) Where, φ(p) is the porosity of the matrix considering stress sensitivity; φ _i is the initial porosity, dimensionless; φ _r is the porosity under maximum stress, dimensionless. k(p) is the matrix permeability after considering stress sensitivity, nD; k _i is the original permeability of matrix, nD; η is the stress sensitivity coefficient of porosity, Pa⁻¹; ψ is the stress-sensitive coefficient of permeability, Pa⁻¹; p _i is the original gas reservoir pressure, Pa.

The effective pore radius of the matrix considering the thickness of the adsorption layer and the stress-sensitive effect is as follows: r e = 8 τ k o n , a b s e − Ψ p i − p φ r + φ i − φ r e − η p i − p (3) Where, τ is the tortuosity, which is the ratio of the actual length of the flow path to the represented length.

In the circular organic nanopores, the apparent permeability contributed by the bulk phase gas based on the Beskok model can be expressed as: k o n , b u l k = − r e 2 φ p 8 τ 1 + α K n 1 + 4 K n 1 + K n (4) Where, Kn is Knudsen diffusion constant, dimensionless; α is the rarefied gas effect coefficient, dimensionless, which can be calculated by the following equation: α = α 0 2 π tan − 1 α 1 K n β (5) Where, α ₀ is the rarefied gas coefficient under K _n→∞ condition; α ₁ and β are fitting constants.

The apparent permeability considering the real gas effect can be expressed as: k o n , s u r f a c e = D s μ g r R T C s c φ p M g p 2 τ (6) Where, D _s is the surface diffusion coefficient, m²/s; μ _gr is the gas viscosity after considering the real gas effect, mPa·s; R is the universal gas constant, Pa/(mol•K); T is the absolute temperature, K; M _g is the molecular weight of gas, g/mol; C _sc is the concentration of adsorbed gas in Langmuir monolayer on the pore surface, kg/m³. The calculation equation is as follows: C s c = 4 θ M g π d M 3 N A (7) Where, N _A is Avogadro’s constant, 6.022×10²³mol^-1.

The surface diffusion coefficient can be calculated as follows: D s = D s 0 1 − θ + κ 2 θ 2 − θ + H 1 − κ 1 − κ κ 2 θ 2 1 − θ + κ 2 θ 2 (8) H 1 − κ = 0 , κ ≥ 1 1,0 ≤ κ ≤ 1 , κ = κ b κ m (9) Where, H(1-κ) is Heaviside function; κ _b is the plugging rate, m/s; κ _m is migration rate, m/s; D s 0 represents the surface diffusion coefficient when the gas coverage is 0, m²/s, which is calculated by the following equation: D S 0 = 8.29 × 10 − 7 T 0.5 e − Δ H 0.8 R T (10) Where, ΔH is the isothermal adsorption heat of gas, J/mol.

The total gas transport flux includes free gas phase transport and surface diffusion of adsorbed gas in circular tube pores. After taking into account the cross-sectional area of the circular tube, the apparent permeability of the matrix pore is as follows: k a p p = ς b u l k k o n , b u l k + ς s u r f a c e k o n , s u r f a c e = r e 2 r e 2 + θ d M k o n , b u l k + 1 − r e 2 r e 2 + θ d M k o n , s u r f a c e (11) Where ς _bulk is body phase gas weight coefficient, its value is equal to the ratio of the cross-sectional area occupied by the body phase gas and the total cross-sectional area, dimensionless; ς _surface is surface diffusion weight coefficient, its value is equal to the ratio of the cross-sectional area occupied by the surface diffusion gas to the total cross-sectional area, dimensionless.

Further, based on the dual media-discrete fracture model, the comprehensive seepage model of fractured horizontal well is established as follows:

①Seepage control equations of matrix system: ∇ ⋅ k a p p μ g B g ∇ p g m l − α l k a p p μ g B g p g m l − p g f l + q d e s = ∂ ϕ m l / B g ∂ t (12) Where, q d e s = − ρ s 1 − ϕ m l − ϕ f l ∂ V E ∂ t (13) Where, B _g is the volume coefficient of gas, dimensionless; p _gm is the matrix system gas pressure, MPa; p _gf is the fracture system gas pressure, MPa; q _des is the desorption production of matrix, m³/s; φ _m is the matrix system porosity, dimensionless; φ _f is the porosity of fracture system, dimensionless; α is the shape factor, dimensionless; ρ _s is the density of shale rock, g/cm³; The subscripts l=1, 2, 1 denotes the inner zone, 2 denotes the outer zone.

②Seepage control equations of fracture system:

For inner zone: ∇ ⋅ k f k r u 1 μ u B u ∇ p f u 1 + α 1 k a p p k r u 1 μ u B u p m u 1 − p f u 1 = ∂ ∂ t S f u 1 ϕ f 1 B u (14)

For outer zone: ∇ ⋅ k f μ g B g ∇ p f g 2 + α 2 k a p p μ g B g p m g 2 − p f g 2 = ∂ ∂ t ϕ f 2 B g (15) Where, p _f is the pressure of the natural fracture system, MPa; k _f is the permeability of natural fractures in shale reservoirs, mD; k _r is the relative permeability, dimensionless; Subscript u=g, w, where g represents the gas phase and w represents the water phase.

③Control equations of hydraulic fracture seepage: ∇ ⋅ k F k r u l μ u B u ∇ p F u + q u s c t = ∂ ∂ t ϕ F S F u B u (16) Where, k _F is hydraulic fracture permeability, mD; φ _F is the hydraulic fracture porosity, dimensionless; q _gsct is horizontal well production, m³/s; p _F is the hydraulic fracture pressure, MPa.

First, the element characteristic matrix of Eq. (12) is constructed based on triangular mesh and control volume finite element method (as shown in Figure 3), and the average mesh pressure p _gm and average saturation S _fu of the matrix system are approximated by triangle vertex interpolation: p g m x , y = ∑ v = i k N v x , y p g m v , p f u l x , y = ∑ v = i k N v x , y p f u l v S f u l x , y = ∑ v = i k N v x , y S f u l v (17) Where N _v is the interpolation function.

For the convection term in the matrix system seepage equations, Galerkin’s weighted margin method and Gauss’s law are used to obtain: ∬ Ω i ∇ ⋅ ( λ – ¯ ⋅ ∇ p g m l ) d V = ∑ v = j , k T i v p i g m l − p v g m l (18) Where T _iv represents the conductivity matrix given by the following equation: T i v = A λ x ∂ N i ∂ x ∂ N v ∂ x + λ y ∂ N i ∂ y ∂ N v ∂ y = 1 4 A λ x a i a l + λ y b i b l (19)

The above equations establish the net flow of the control volume of the inflow node i of triangle ijk. Similarly, the net inflow flow equation of nodes j and k can be established. Then the flow matrix of the virtual finite control volume unit can be obtained as: T i j + T i k − T i j − T i k − T j i T j i + T j k − T j k − T k i − T k j T k i + T k j p i g m l p j g m l p k g m l = f g m i f g m j f g m k (20)

Eq. (20) is the processing format of the convection term of the controlled volume finite element method based on triangular mesh. In Eq. (12), the cross-flow term between matrix system and fracture system, as well as the adsorption and desorption term and the time derivative term of matrix system are numerically calculated using the traditional finite element method.

Taking gas phase as an example, the finite element integral numerical calculation format for constructing the cross-flow term is as follows: ∬ Ω m N v α k a p p μ g B g p m − p f d Ω m = ∬ Ω m N v α k a p p μ g B g ∑ v = i k N v p m v − ∑ v = i k N v p f v d Ω m = α A k a p p 3 μ g B g i 0 0 0 α A k a p p 3 μ g B g j 0 0 0 α A k a p p 3 μ g B g k p m i − p f i p m j − p f j p m k − p f k (21)

For the adsorption and desorption term and the time derivative term, the finite element integral numerical calculation format is constructed by combining Euler backward difference as follows: A 3 Δ t ϕ m B g i n + 1 − ϕ m B g i n ϕ m B g j n + 1 − ϕ m B g j n ϕ m B g k n + 1 − ϕ m B g k n + A 3 Δ t 1 − ϕ m − ϕ f V E i n + 1 − 1 − ϕ m − ϕ f V E i n 1 − ϕ m − ϕ f V E j n + 1 − 1 − ϕ m − ϕ f V E j n 1 − ϕ m − ϕ f V E k n + 1 − 1 − ϕ m − ϕ f V E k n (22)

Eqs 17–22 construct the element characteristic matrix of shale matrix gas phase seepage equation based on triangular mesh and controlled volume finite element method. Similarly, the element characteristic matrix of shale fracture system seepage equation can be constructed.

For hydraulic fracturing fractures, the discrete fracture model is used to reduce their dimensionality to a one-dimensional line segment with a certain width. The interpolation function of the gas phase equation can be expressed as: p g F = x j − x Δ l p g F i + x − x i Δ l p g F j = N i p g F i + N j p g F j = N i   N j p g F i p g F j (23)

Combined with the interpolation function, Eq. (16) of the hydraulic fracture seepage control equation is constructed in the following numerical calculation format: ∬ Ω F N v ∂ ∂ l F λ g F ∂ p g F ∂ l F d Ω F = w F × Δ l F × λ g F × ∂ N v ∂ l F ∂ N i ∂ l F   ∂ N j ∂ l F p g F i p g F j = w F × λ g F × 1 Δ l F − 1 Δ l F − 1 Δ l F   1 Δ l F p g F i p g F j (24)

For the time term: ∬ Ω F N l ∂ S g F ϕ F / B g F ∂ t d Ω F = w F × Δ l F 2 Δ t × S g F ϕ F B g F i n + 1 − S g F ϕ F B g F i n S g F ϕ F B g F j n + 1 − S g F ϕ F B g F j n (25)

Similarly, the characteristic matrix of water phase unit in fracturing fracture system can be further constructed.

The pressure value of each system in the next time step is characterized by the fully implicit calculation format. The “dimensionality reduction” treatment of fractured fractures is carried out and embedded into the reservoir flow model. By assembling the characteristic matrix of each unit, the fully implicit iterative solution matrix for the whole region is constructed as follows: T g f k + δ T g f k P g f k − W g m f k − δ N g f k Δ t   δ W g m f k P g m k − P g f k + W g m f k   δ T g f k P g f k + δ N g f − S w f k Δ t W g m f k   T g m k + δ T g m k − δ W g m f k P g m k − P g f k − W g m f k − δ N g m k Δ t − δ V k Δ t T w f k   δ T g f k P g f k + δ N w f k Δ t × δ P g f k δ P g m k δ S w f k = R g f k R g m k R w f k (26)

Finally, the numerical well model of fractured horizontal well is established and put into Eq. (26): q u s c t n + 1 = ∑ i = 1 N F P I u i n + 1 p b h − p a v e i n + 1 (27) Where, q _usct is the daily gas and water production ground conditions, m³/d; p _bh is bottom-hole flow pressure, MPa; p _avei is the average pressure of the grid block where the well source and junction are located, MPa; N _F is the total number of intersection points between fractured fractures and horizontal wellbore. PI _ui is the production index, m³/d/MPa. Through iterative computations, when the accuracy requirements of δP and δS_w are met, the stable value under current time step is obtained and the calculation of the next time step is started.

For a fractured horizontal well with a complex fracture network (Figure 1) in the south Sichuan shale gas reservoir block, the triangular grid was used to discrete the spatial grid of the whole area, and a fully implicit numerical simulation program was developed to simulate and predict the production dynamic characteristics of fractured horizontal wells. The basic parameters of the model are shown in Table 1. The calculation equation of the fracture relative permeability curve is: k _rgf=1-S _wf, k _rwf=S _wf.

Figure 4 shows the change of gas well pressure of fractured horizontal wells with complex fracture networks of shale gas for 100 days, 1 year and 10 years, respectively. As can be seen from the figure, the pressure propagation presents the characteristics of non-uniform distribution in the complex fracture network area, and the pressure propagation is faster in the well section with high reconstruction degree. At the same time, the gas well pressure propagation gradually expands outward with time, and the pressure at the boundary of the inner and outer regions is also non-uniform due to the influence of the complex fracture network boundary shape. When the pressure propagates to the original reservoir area, its pressure propagation speed is significantly slower than that of the reconstruction area, indicating that the reconstruction area of complex fracture network is the main contribution area to the pressure drop and production of shale gas wells.

In addition, Figure 5 shows the production decline curves of fractured horizontal wells with and without SRV region. By comparison, it can be seen that the SRV region formed by complex fracture networks will significantly increase the daily and cumulative gas production of fractured wells, indicating that compared to traditional double wing hydraulic fracturing method, fracture network fracturing is more conducive to achieving higher single well EUR and production efficiency for horizontal wells in shale gas reservoirs.

Figure 6 shows the influence of different fracture water saturation on gas well productivity. Due to the low flowback rate of fracturing fluid in shale reservoirs, a large amount of fracturing fluid remains in the reservoir, resulting in gas-water two-phase flow characteristics in fractured wells and near wellbores. The simulation results show that the higher the fracture water saturation, the lower the daily gas production. When the fracture water saturation is 0.2, 0.3, and 0.4, the cumulative gas production after 10 years of depletion production is 194 million m³, 170 million m³and 145 million m³, respectively. Therefore, the retained fracturing fluid has a significant impact on gas well productivity.

Figure 7 shows the influence of different hydraulic fracture conductivity on gas well productivity. As can be seen from the figure, the higher the diversion capacity, the higher the early daily gas output, but the difference degree gradually decreases with the increase of production time. When the diversion capacity CFD is 0.25D cm, 0.5D cm and 1.0D cm, the cumulative gas production after 10 years of depletion is 133 million m³, 180 million m³ and 194 million m³, respectively. Therefore, in the hydraulic fracturing process of shale gas reservoir, it is necessary to improve the conductivity of fracturing fractures as much as possible to increase the production.

Figure 8 shows the influence of permeability in the stimulation area on gas well productivity. It can be seen that the higher the permeability of the stimulation area, the larger the daily gas production of the gas well in the early stage, and the larger the cumulative production of the gas well, indicating that the fracturing stimulation area around the horizontal well is more fully stimulated. When the fracture permeability in the stimulation area is 0.001mD, 0.005mD and 0.01mD, the cumulative gas production after 10 years of depletion production is 194 million m³, 261 million m³ and 308 million m³, respectively. It is the key to improve the production of shale gas wells by optimizing the construction parameters, increasing the stimulation degree of the stimulation area and improving the seepage capacity of shale reservoir.

Figure 9 shows the influence of different Langmuir volume V _L on gas well production curves. As can be seen from the figure, the larger the V _L value, the larger the daily gas production and cumulative gas production of the gas well, and the slower the decline rate of gas well production. With the progress of gas well production, the free gas in fracture and matrix system is preferentially produced. The continuous reduction of reservoir pressure makes the gas adsorbed on the surface of organic matter begin to desorption, and enter the fracture system as a reservoir supplementary gas source, delaying the decline of gas well production. When Langmuir volume V _L is 0.01 m³/kg, 0.02 m³/kg and 0.03 m³/kg, the cumulative gas production after 10 years of depletion production is 193 million m³, 260 million m³ and 313 million m³, respectively. Adsorption and desorption is a unique seepage mechanism of shale gas reservoir, which is different from conventional gas reservoir. With the progress of production, the reservoir pressure decreases and the adsorption gas is desorbed in large quantities, which will contribute to the production replenishment and stable production of fractured horizontal wells in shale gas reservoir in the middle and later stages.