A Transient Productivity Model of Fractured Wells in Shale Reservoirs Based on the Succession Pseudo-Steady State Method

Abstract

After volume fracturing, shale reservoirs can be divided into nonlinear seepage areas controlled by micro- or nanoporous media and Darcy seepage areas controlled by complex fracture networks. In this paper, firstly, on the basis of calculating complex fracture network permeability in a stimulated zone, the steady-state productivity model is established by comprehensively considering the multi-scale flowing states, shale gas desorption and diffusion after shale fracturing coupling flows in matrix and stimulated region. Then, according to the principle of material balance, a transient productivity calculation model is established with the succession pseudo-steady state (SPSS) method, which considers the unstable propagation of pressure waves, and the factors affecting the transient productivity of fractured wells in shale gas areas are analyzed. The numerical model simulation results verify the reliability of the transient productivity model. The results show that: (1) the productivity prediction model based on the SPSS method provides a theoretical basis for the transient productivity calculation of shale fractured horizontal well, and it has the characteristics of simple solution process, fast computation speed and good agreement with numerical simulation results; (2) the pressure wave propagates from the bottom of the well to the outer boundary of the volume fracturing zone, and then propagates from the outer boundary of the fracturing zone to the reservoir boundary; (3) with the increase of fracturing zone radius, the initial average aperture of fractures, maximum fracture length, the productivity of shale gas increases, and the increase rate gradually decreases. When the fracturing zone radius is 150 m, the daily output is approximately twice as much as that of 75 m. If the initial average aperture of fractures is 50 μm, the daily output is about half of that when the initial average aperture is 100 μm. When the maximum fracture length increases from 50 m to 100 m, the daily output only increases about by 25%. (4) When the Langmuir volume is relatively large, the daily outputs of different Langmuir volumes are almost identical, and the effect of Langmuir volume on the desorption output can almost be ignored.

1. Introduction

Shale gas reservoirs have the characteristics of tiny pores and throats, extremely low permeability, abundant natural fractures and diverse gas storage modes [1,2,3,4,5,6,7,8,9]. After the reservoir is fractured, a complex fracture network is generated, composed of activated natural fractures and induced hydraulic fractures, and the gas flow channels are converted from nanoscale pores to microscale pores so that the gas flow state in the stimulated zone is changed [10,11,12,13,14,15,16,17,18]. The Darcy flow in the stimulated region near the wellbore is coupled with the nano-/microscale flow in the matrix, resulting in the complex multiscale flow in shale gas reservoirs [19]. Because interactions between desorption, diffusion and seepage exist in reservoir shale gas, the original linear seepage theory is no longer applicable as a result of the change of flow patterns [20]. Therefore, it is necessary to establish a new transient productivity calculation model of fractured wells considering of multi-scale flow and nonlinear seepage.

Currently, there are few studies on the transient productivity of fractured wells in shale gas reservoirs, and the characteristics of multi-scale flow and non-linear flow are neglected in the majority of productivity models [21,22,23]. The current analytical models are based on Laplace transform, superposition principle, source function, complex function theory and other mathematical methods to analyze the transient pressure response for shale reservoir, while the Laplace transform and source function are too complicated to solve. Meanwhile, the assumption that the pressure has already reached the boundary of the reservoir under initial condition is inconsistent with the actual situation. The new model established in this paper can avoid these disadvantages.

For the complex fracture network intermingled with natural fractures and artificial fractures near the wellbore [24,25,26], Min et al. [27] used the cubic law to calculate the permeability of complex fracture networks, and fracture closure and shear failure of complex fracture networks under different stresses are considered. In a study of the unsteady productivity prediction of shale gas, Stalgorova et al. [28] divided the seepage area of shale reservoir into stimulated region and un-stimulated region, and a pressure transient model based on tri-line flow model was established. However, this model neglected shale gas desorption and diffusion effects so that it was inaccurate for shale gas reservoir productivity evaluation. On the basis of the tri-line flow model, Swami et al. [29] established a dual-porosity seepage model of shale gas reservoir considering adsorption and desorption in matrix. However, this model still assumed that shale gas is linear seepage and the characteristic of multi-scale flow is neglected. Based on the previous studies, Deng [12] corrected the B-K model with different slip coefficients, and established a multi-scale flow model for shale gas reservoir considering diffusion, slippage, desorption and adsorption, while the assumption of homogeneous medium does not match with volume fracturing and not take into account the unsteady seepage of shale reservoir. Zhang et al. [30] and Su et al. [31] extended the dual-region composite flow model proposed by Zhao et al. [32] to the shale gas productivity model. The transient pressure response model of horizontal well was established by point-source function and Laplace transform, while the inappropriate assumption that the initial pressure extends to the reservoir boundary at the beginning was made. In addition, some researchers have studied other methods for calculating fracturing well properties and productivity [33,34,35,36,37]. With the production of shale gas, the effective stress of the reservoir increases, and a stress-sensitive phenomena occur. Meanwhile, adsorbed shale gas will undergo surface diffusion in addition to desorption [38,39,40].

In order to overcome the deficiencies of the source function and other mathematical methods, Shahamat et al. [41] analyzed the pressure response under the condition of transition flow and boundary flow with the succession pseudo-steady state (SPSS) method, which considers the correlation of pressure wave propagation with reservoir properties, fluid properties, and production time. The SPSS method is simple and easy to implement programming, which avoids the unreasonable assumption that the pressure has already reached the boundary of the reservoir under the initial condition. However, Shahamat did not take into consideration the multi-scale flow and non-linear flow after fracturing in the shale gas reservoir.

In this paper, first, based on the Beskok-Karniadakis apparent permeability model in regard of different slip coefficients, the seepage region in shale gas reservoir is divided into a stimulated zone and a matrix zone. The radial Darcy flow is considered in the stimulated zone and non-linear flow characterized by Knusen Number is considered in the matrix zone. Second, the mathematical model of steady-state productivity of shale horizontal wells with volume fracturing is derived. Then, a transient productivity calculation model of fractured wells in shale combined the material balance equation is established with the SPSS method, and the numerical model simulation results verify the reliability of the transient productivity model. Finally, the horizontal well productivity prediction and the analysis of influencing factors are carried out.

2. Physical Model and Basic Assumptions

As is shown in Figure 1, based on the fracture network after the volume fracturing in shale reservoir, a composite flow model for fractured horizontal well is established, and the seepage regions in shale reservoir are divided into inside flow region and outside flow region [41].

The basic assumptions are made for multi-stage fractured horizontal wells:

(1)

The model is for isothermal single-phase shale gas flow, and vertical flow is neglected.

(2)

The reservoir is composite with the matrix zone and stimulated zone, and the reservoir has a constant and uniformed thickness with the upper and lower boundaries closed.

(3)

The gas seepage is characterized by Kundsen number in matrix zone with the radius of $r_{\mathrm{e}}$ while the gas flow is consistent with Darcy law in stimulated zone with the radius of $r_{\mathrm{f}}$.

3. Mathematical Model

3.1. Steady-State Productivity Model

The shale reservoir after volume fracturing is divided into stimulated zone and matrix seepage zone, and the gas flow equation at different zone is established respectively. On this basis, the mathematical model of steady-state productivity of shale horizontal well with volume fracturing is derived.

3.1.1. Shale Matrix Gas Seepage Model

Knudsen dimensionless number K_n is defined as:

$$K_{\mathrm{n}}=\frac{\overline{\lambda }}{r_{\mathrm{p}}}$$

(1)

where $\overline{\lambda }$ is the gas molecules’ mean free path, nm; r_p is the pore radius, nm.

With regard to the shale matrix gas seepage model, Beskok and Karniadakis [42] established an ideal gas flow equation that is universally applicable to continuous flow, slip flow, transition flow, and molecular flow:

$$v=-\frac{K_{\mathrm{m}}}{\mu }(1+\alpha K_{\mathrm{n}})(1+\frac{4K_{\mathrm{n}}}{1-\mathrm{b}{K}_{\mathrm{n}}})\frac{\mathrm{d}p}{\mathrm{d}x}$$

(2)

where v is the gas seepage velocity, m/s; K_m is the permeability of matrix, 10⁻³μm²; $\mu$ is the viscosity of shale gas, mPa·s; α is the rarefaction effect factor, dimensionless; K_n is the Knudsen number, dimensionless; b is the gas slippage constant, dimensionless.

Knudsen diffusion coefficient is given by Faruk [43]:

$$D_{\mathrm{k}}=\frac{4r_{\mathrm{p}}}{3}\sqrt{\frac{2Z{\mathrm{R}}_{\mathrm{g}}T}{\pi M_{\mathrm{w}}}}\times {10}^{-3}$$

(3)

where D_k is the diffusion coefficient, mm²/s; Z is the gas deviation factor, dimensionless; R_g is the universal gas constant, J·mol⁻¹·K⁻¹; T is the shale formation temperature, K; M_w is the molar mass of the gas, kg·mol⁻¹.

The permeability in shale matrix is defined as [44]:

$$K_{\mathrm{m}}=\frac{r_{\mathrm{p}}{}^2{\varphi }_{\mathrm{m}}}{8{\tau }^2}\times {10}^{-3}$$

(4)

where ${\varphi }_{\mathrm{m}}$ is the matrix porosity of shale, dimensionless; $\tau$ is the tortuosity, dimensionless (The value is 1).

The gas mean free path of a molecule is defined as [28]:

$$\overline{\lambda }=\sqrt{\frac{\pi Z{\mathrm{R}}_{\mathrm{g}}T}{2M_{\mathrm{w}}}}\frac{\mu }{p}$$

(5)

where p is the formation pressure, MPa.

By combining Equation (1), Equation (3), Equation (4), and Equation (5), K_n can be expressed as follows:

$$K_{\mathrm{n}}=\frac{3\pi \mu {\varphi }_{\mathrm{m}}{D}_{\mathrm{k}}}{64K_{\mathrm{m}}p}$$

(6)

If α = 0, b = 1, combining Equations (2) and (6), differential equation of gas seepage in shale matrix pores can be expressed as follows:

$$v=-\frac{K_{\mathrm{m}}}{\mu }(1+\frac{3\pi \mu {\varphi }_{\mathrm{m}}{D}_{\mathrm{k}}}{16K_{\mathrm{m}}p})\frac{\mathrm{d}p}{\mathrm{d}x}$$

(7)

The seepage mathematical model in the external matrix area can be expressed as follows:

$$v=-\frac{K_{\mathrm{m}}}{\mu }(1+\frac{3\pi \mu {\varphi }_{\mathrm{m}}{D}_{\mathrm{k}}}{16K_{\mathrm{m}}p})\frac{\mathrm{d}p}{\mathrm{d}r}$$

(8)

3.1.2. Stimulated Region Gas Seepage Model

After volume fracturing, the inner area will form a complex fracture network containing natural fractures and artificial fractures [45]. Assuming the fracture fluid flows in a smooth plane, the permeability is calculated by cubic law according to the projection length of each fracture on the x and y axis. By superimposing the permeability of each fracture, the permeability of the whole fracture network of x and y directions without any stress is calculated.

$$k_x={\displaystyle \sum _{i=1}^N\frac{b_i{}^3{l}_i}{12A}\cdot \cos {\theta }_i}$$

(9)

$$k_y={\displaystyle \sum _{i=1}^N\frac{b_i{}^3{l}_i}{12A}\cdot \sin {\theta }_i}$$

(10)

where k_x and k_y are the permeability in x direction and y direction respectively, 10⁻³ μm²; N is the number of fractures, dimensionless; A is the area of fractured region, mm²; b_i is the aperture of different fracture, μm; l_i is the length of different fracture, mm; θ_i is the angle of the fracture and the direction of X, °.

The permeability of fractures under different stress is composed of normal open fracture and shear type fracture. In the calculation of permeability, the normal stress and shear stress are calculated, respectively, and the ultimate permeability is the superposition of the permeability caused by the normal stress and shear stress [27]:

$$k_x=k_{nx}+k_{dx}$$

(11)

$$k_y=k_{ny}+k_{dy}$$

(12)

where k_nx and k_ny are the permeability in the x and y direction due to normal closure of fractures respectively, 10⁻³ μm²; k_dx and k_dy are the permeability in the x and y direction due to shear failure of fractures, respectively, in 10⁻³ μm².

In the fracture network, single fractures do not reflect fractures with different lengths and apertures. Therefore, the concept of equivalent fracture frequency (f_x, f_y) and equivalent aperture (b_x, b_y) are used to express the complex fracture deformation in fracture network. The permeability of the open or closed fractures produced by the positive stress in the direction of X and Y (k_nx, k_ny) can be calculated by the following equations [46]:

$$k_{nx}=\frac{f_x}{12}{b}_x^3$$

(13)

$$k_{ny}=\frac{f_y}{12}{b}_y^3$$

(14)

where f_x and f_y are the equivalent frequency caused by the positive stress in the x direction and y direction, respectively, 1/m; b_x and b_y are the equivalent aperture caused by the positive stress in the x direction and y direction respectively, in μm.

The equivalent frequency of fractures can be obtained by numerical experiments and initial aperture inversion.

$$f_x=\frac{12k_{nx}}{b_{ei}^3}$$

(15)

$$f_y=\frac{12k_{ny}}{b_{ei}^3}$$

(16)

where b_ei is the initial average aperture of fractures, μm.

The size of the fracture aperture in the x and y direction can be expressed as follows respectively [47]:

$$b_x=b_r+b_{\max }=b_r+b_{\max }\exp \{-({\alpha }_x{\sigma }_x+{\alpha }_y{\sigma }_y)\}$$

(17)

$$b_y=b_r+b_{\max }=b_r+b_{\max }\exp \{-({\beta }_x{\sigma }_x+{\beta }_y{\sigma }_y)\}$$

(18)

where b_r is the residual aperture, μm; b_max is the maximum deformable mechanical aperture caused by positive stress, μm; α_x and α_y are the x aperture stress coefficient in the x and y direction, respectively, dimensionless; β_x and β_y are the y aperture stress coefficient in the x and y direction, respectively, dimensionless; ${\sigma }_x$ and ${\sigma }_y$ are the normal stress in the x direction and y direction, respectively, MPa.

Similar equations to Equation (13) and Equation (14), fracture permeability under shear stress can be expressed as follows:

$$k_{dx}=\frac{f_{dx}}{12}{d}_x^3$$

(19)

$$k_{dy}=\frac{f_{dy}}{12}{d}_y^3$$

(20)

where f_dx and f_dy are the equivalent frequency caused by the shear stress in the x direction and y direction, respectively, 1/m; d_x and d_y are the equivalent apertures caused by the shear stress in the x direction and y direction respectively, μm.

Only when the pressure reaches the critical pressure, the fracture will be shear failure, so only some fractures will appear shear deformation, and other fractures will not appear shear deformation. The aperture produced by shear deformation can be calculated by the following formulae:

$$k\prec k_c,d_x=0,d_y=0$$

(21)

$$k\succ k_c,\{\begin{array}{l}{d}_x=d_{\max }\{1-\exp \left[-{\gamma }_x(k-k_c)\right]\}\\ d_y=d_{\max }\{1-\exp \left[-{\gamma }_y(k-k_c)\right]\}\end{array}$$

(22)

where k is the stress ratio of x and y, dimensionless; k_c is the critical stress ratio, dimensionless; d_max is the maximum aperture of fracture after shear failure, μm; γ_x and γ_y are the shear stress coefficient in x direction and y direction respectively, dimensionless;

Assuming that the fracture is not cohesive, the critical pressure ratio and critical direction of fracture failure can be calculated by the following Coulomp failure criterion [48]:

$$k_c=\frac{1+\sin \phi }{1-\sin \phi }$$

(23)

$${\phi }_{\mathrm{f}}=90-\left(45+\frac{\phi }{2}\right)$$

(24)

where φ is the internal friction angle of rock, °; φ_f is the critical failure angle of rock, °.

Finally, the permeability (k_x, k_y) of the x and y directions in the stimulated area can be obtained by superimposing the fracture permeability of normal closure and shear deformation [27]:

$$k_x=\frac{f_x}{12}{\left[b_r+b_{\max }\exp \{-({\alpha }_x{\sigma }_x+{\alpha }_y{\sigma }_y)\}\right]}^3+\frac{f_{dx}}{12}{\left[d_{\max }\left[1-\exp \{-{\gamma }_x(k-k_c)\}\right]\right]}^3$$

(25)

$$k_y=\frac{f_y}{12}{\left[b_r+b_{\max }\exp \{-({\beta }_x{\sigma }_x+{\beta }_y{\sigma }_y)\}\right]}^3+\frac{f_{dy}}{12}{\left[d_{\max }\left[1-\exp \{-{\gamma }_y(k-k_c)\}\right]\right]}^3$$

(26)

The comprehensive permeability of stimulated fracture network area is obtained through the root mean square of x and y permeability [49]:

$$K_{\mathrm{f}}=\frac{\sqrt{k_x^2+k_y^2}}{2}$$

(27)

where K_f is the permeability in volume fracturing zone, 10⁻³ μm².

According to the hypothesis, the Darcy seepage equation in the internal volume fracturing area can be expressed as follows:

$$\frac{\mathrm{d}p}{\mathrm{d}r}=-\frac{\mu }{K_{\mathrm{f}}}v$$

(28)

3.1.3. Steady-State Productivity Model

Considering the flow is the product of gas seepage velocity and gas seepage area, the flow of volume fracturing area can be expressed as follows:

$$q_{1\mathrm{sc}}=\frac{\pi K_{\mathrm{f}}h{T}_{\mathrm{sc}}(p_{\mathrm{f}}{}^2-p_{\mathrm{wf}}{}^2)}{\overline{\mu }\overline{Z}{p}_{\mathrm{sc}}T\ln \frac{r_{\mathrm{f}}}{r_{\mathrm{w}}}}$$

(29)

where $q_{1\mathrm{sc}}$ is the volume fracturing zone flow, m³/d; h is the reservoir thickness, m; T_sc is the temperature under standard condition, K; $p_{\mathrm{f}}$ is the outer boundary pressure of volume fracturing zone, MPa; $p_{\mathrm{wf}}$ is the bottom hole flow pressure, MPa; $r_{\mathrm{f}}$ is the radius of volume fracturing zone, m; $r_{\mathrm{w}}$ is the borehole radius, m; p_sc is the pressure under standard conditions, MPa; $\overline{\mu }$ is the gas viscosity under average formation pressure, m Pa·s; $\overline{Z}$ is the compression factor under average formation pressure, dimensionless.

The flow of external matrix seepage area is obtained from Equation (8):

$$q_{2\mathrm{sc}}=\frac{\pi K_{\mathrm{m}}h{T}_{\mathrm{sc}}(p_{\mathrm{e}}{}^2-p_{\mathrm{f}}{}^2)}{\overline{\mu }\overline{Z}{p}_{\mathrm{sc}}T\ln \frac{r_{\mathrm{e}}}{r_{\mathrm{f}}}}+\frac{\pi K_{\mathrm{m}}{T}_{\mathrm{sc}}h(3\pi \overline{\mu }{\varphi }_{\mathrm{m}}{D}_{\mathrm{k}})(p_{\mathrm{e}}-p_{\mathrm{f}})}{8K_{\mathrm{m}}{p}_{\mathrm{sc}}T\overline{\mu }\overline{Z}\ln \frac{r_{\mathrm{e}}}{r_{\mathrm{f}}}}$$

(30)

where $q_{2\mathrm{sc}}$ is the matrix seepage area flow, m³/d; $r_{\mathrm{e}}$ is the supply radius, m.

According to the law of mass conservation [50], the gas volume flow of the matrix zone is equal to that of volume fracturing zone under the standard conditions:

$$q_{\mathrm{sc}}=q_{2\mathrm{sc}}=q_{1\mathrm{sc}}$$

(31)

where $q_{\mathrm{sc}}$ is the production of shale gas under the standard condition, m³/d.

By combining Equation (29), Equation (30), and Equation (31), the volume flow for steady seepage in fractured horizontal well in shale gas reservoir can be expressed as follows:

$$q_{\mathrm{sc}}=\frac{p_{\mathrm{e}}{}^2-p_{\mathrm{wf}}{}^2}{A+BC}$$

(32)

where $A=\frac{\overline{\mu }\overline{Z}{p}_{\mathrm{sc}}T\ln \frac{r_{\mathrm{f}}}{r_{\mathrm{w}}}}{\pi K_{\mathrm{f}}h{T}_{\mathrm{sc}}}$; $B=\frac{\overline{\mu }\overline{Z}{p}_{\mathrm{sc}}T\ln \frac{r_{\mathrm{e}}}{r_{\mathrm{f}}}}{\pi K_{\mathrm{m}}h{T}_{\mathrm{sc}}}$; $C=\frac{p_{\mathrm{e}}{}^2-p_{\mathrm{f}}{}^2}{p_{\mathrm{e}}{}^2-p_{\mathrm{f}}{}^2+\frac{3\pi \mu {\varphi }_{\mathrm{m}}{D}_{\mathrm{k}}(p_{\mathrm{e}}-p_{\mathrm{f}})}{8K_{\mathrm{m}}}}$.

3.2. Unsteady-State Productivity Model

The seepage of shale gas is an unstable process, and the initial production decline is very obvious. Only when the pressure wave reaches to the boundary in the later stage of production, the production tends to be stable. In shale reservoir, the amount of dissolved gas is fairly small, so the effect of dissolved gas is not considered in the process of deriving the unstable mathematical model, and the expansion of the adsorbed gas volume with pressure drop is ignored.

Considering the effects of shale gas adsorption and desorption, the Langmuir isotherm equation [51] represents the desorbed volume of shale gas from the matrix as pressure changes:

$$Q=V_{\mathrm{L}}(\frac{p_{\mathrm{i}}}{p_{\mathrm{L}}+p_{\mathrm{i}}}-\frac{p}{p_{\mathrm{L}}+p})$$

(33)

where Q is the desorption gas volume of per kilogram mass shale matrix, m³/kg; $V_{\mathrm{L}}$ is the Langmuir volume, m³/kg; $p_{\mathrm{i}}$ is the original formation pressure, MPa; $p_{\mathrm{L}}$ is the Langmuir pressure, MPa.

From the initial stage of production to the end of production, the propagation of pressure waves in the reservoir can be divided into two stages. In the first stage, the pressure wave propagates from the bottom of the well to the boundary of the reservoir, and the pressure wave radius gradually increases to $r_{\mathrm{e}}$. In the second stage, when the pressure wave reaches the boundary, the pressure wave radius no longer changes, and the boundary pressure gradually decreases. Meanwhile, the gas production decreases, and finally, the pressure on the boundary tends to the bottom hole flow pressure. In this paper, the first stage of pressure wave propagation is divided into two periods, that is, the pressure wave firstly propagates from the bottom of the well to the outer boundary of the volume fracturing zone, and then it propagates from the outer boundary of the fracturing zone to the reservoir boundary. The distance of pressure wave propagation is only related to time, the physical properties of reservoir and the fluid.

In order to derive a transient productivity model for shale reservoir fractured well, this paper adopts a succession pseudo-steady state method (SPSS) proposed by Shahamat [41]. Specifically, on the basis of the pressure wave radius formula, steady state shale gas production formula, and material balance equation, assuming a time step, it is considered that the seepage within this time step is a steady flow, and the cumulative output is calculated. The material balance equation can calculate the average formation pressure in the area affected by the pressure wave during this period of time. The average formation pressure is used as the boundary pressure for the next time step to calculate the output at the next time step. By analogy, the relationship between output and time is obtained, that is, the output under non-steady state. The specific steps are as follows:

3.2.1. The Solution of Initial Production

The whole production stages of fractured horizontal well in shale reservoir are divided into several time steps. For the first time step Δt, Hsieh et al. [52] proposed the formula for calculating the propagation radius of pressure wave at the time of Δt. The value of the Δt is as small as possible in order to avoid R₁ from exceeding the outer boundary of the volume fracturing area:

$$R_1=0.5879\sqrt{\frac{K_{\mathrm{f}}\mathsf{\Delta }t}{{\varphi }_{\mathrm{f}}\mu C_t{}_{\mathrm{f}}}}$$

(34)

where R₁ is the propagation radius of pressure wave at the time of Δt, m; Δt is the time of a production time step, d; ${\varphi }_{\mathrm{f}}$ is the fracture porosity in volume fracturing area, dimensionless; C_tf is the comprehensive compression coefficient of fractured zone, MPa⁻¹.

According to the formula of steady-state productivity, the production at the time of Δt can be expressed as follows:

$$q_1=\frac{\pi K_{\mathrm{f}}h{T}_{\mathrm{sc}}(p_i{}^2-p_{\mathrm{wf}}{}^2)}{\overline{\mu }\overline{Z}{p}_{\mathrm{sc}}T\ln \frac{R_1}{r_{\mathrm{w}}}}$$

(35)

where q₁ is the production at the time of Δt, m³.

3.2.2. The Solution of Production at the Next Production Time Step

Using q₁ as initial output, the pressure wave propagation distance R₂ at the time of Δt + Δt_f is obtained:

$$R_2=0.5879\sqrt{\frac{K_{\mathrm{f}}(\mathsf{\Delta }t+\mathsf{\Delta }{t}_{\mathrm{f}})}{{\varphi }_{\mathrm{f}}\mu C_t{}_{\mathrm{f}}}}$$

(36)

where R₂ is the propagation radius of pressure wave at the time of Δt + Δt_f, m; Δt_f is the propagation time of pressure waves in fractured zone, d.

In the time of the Δt_f, the seepage of shale gas is stable seepage, and the cumulative output Gp₂ can be expressed as follows:

$$G_{\mathrm{p}2}=q_1\times \mathsf{\Delta }{t}_{\mathrm{f}}$$

(37)

where Gp₂ is the cumulative yield in the time of the Δt_f, m³.

The free gas geological reserves in the pressure wave propagation radius is calculated by the volume method:

$$G_{\mathrm{m}2}=\pi R_2{}^{2}h{\varphi }_{\mathrm{m}}(1-S_{\mathrm{w}})/B_{\mathrm{gi}}$$

(38)

$$G_{\mathrm{f}2}=\pi R_2{}^{2}h{\varphi }_{\mathrm{f}}(1-S_{\mathrm{w}})/B_{\mathrm{gi}}$$

(39)

where G_m2 is the free gas volume in the shale matrix in the pressure wave radius at the time of Δt + Δt_f, m³; S_w is the irreducible water saturation, dimensionless; B_gi is the volume coefficient of shale gas under the condition of original formation, dimensionless; G_f2 is the free gas volume in the fractures within pressure wave radius at the time of Δt + Δt_f, m³.

According to the principle of material balance, under the ground standard condition, material balance equation can be expressed as follows [53]:

$$G_{\mathrm{m}}+G_{\mathrm{f}}+G_{\mathrm{a}}=G_{\mathrm{p}}+G_{\mathrm{m}}^{\prime }+G_{\mathrm{f}}^{\prime }+G_{\mathrm{a}}^{\prime }$$

(40)

where $G_{\mathrm{m}}$ is the free gas volume in the matrix under the original formation pressure, m³; $G_{\mathrm{f}}$ is the free gas volume in the fractures under the original formation pressure, m³; $G_{\mathrm{a}}$ is the adsorptive gas volume under the original formation pressure, m³; $G_{\mathrm{p}}$ is the produced gas volume, m³; $G_{\mathrm{m}}^{\prime }$ is the free gas volume in the matrix under the current formation pressure, m³; $G_{\mathrm{f}}^{\prime }$ is the free gas volume in the fractures under the current formation pressure, m³; $G_{\mathrm{a}}^{\prime }$ is the adsorptive gas volume under the current formation pressure, m³.

Equation (40) can be rewritten as follows:

$$\begin{array}{l}{G}_{\mathrm{m}}+G_{\mathrm{f}}+\frac{G_{\mathrm{m}}{B}_{\mathrm{gi}}}{{\varphi }_{\mathrm{m}}(1-S_{\mathrm{w}})}{\rho }_{\mathrm{s}}\frac{V_{\mathrm{L}}{p}_{\mathrm{i}}}{p_{\mathrm{L}}+p_{\mathrm{i}}}\\ =G_{\mathrm{p}}+G_{\mathrm{m}}\frac{B_{\mathrm{gi}}}{B_{\mathrm{g}}}\left(1-\frac{C_{\mathrm{m}}+C_{\mathrm{w}}{S}_{\mathrm{w}}}{1-S_{\mathrm{w}}}\mathsf{\Delta }p\right)+G_{\mathrm{f}}\frac{B_{\mathrm{gi}}}{B_{\mathrm{g}}}+\frac{G_{\mathrm{m}}{B}_{\mathrm{gi}}}{{\varphi }_{\mathrm{m}}(1-S_{\mathrm{w}})}{\rho }_{\mathrm{s}}\frac{V_{\mathrm{L}}p}{p_{\mathrm{L}}+p}\end{array}$$

(41)

The equation about the average formation pressure p₂ in the range of pressure propagation radius (R₂) can be expressed as follows:

$$\begin{array}{l}\frac{p_2}{Z}[G_{\mathrm{m}2}\left(1-\frac{C_{\mathrm{m}}+C_{\mathrm{w}}{S}_{\mathrm{w}}}{1-S_{\mathrm{w}}}\left(p_{\mathrm{i}}-p_2\right)\right)+G_{\mathrm{f}2}]=\\ \frac{p_{\mathrm{i}}}{Z_{\mathrm{i}}}[G_{\mathrm{m}2}+G_{\mathrm{f}2}-G_{\mathrm{p}2}+\frac{G_{\mathrm{m}2}{B}_{\mathrm{gi}}{\rho }_{\mathrm{s}}}{{\varphi }_{\mathrm{m}}(1-S_{\mathrm{w}})}(\frac{V_{\mathrm{L}}{p}_{\mathrm{i}}}{p_{\mathrm{L}}+p_{\mathrm{i}}}-\frac{V_{\mathrm{L}}p}{p_{\mathrm{L}}+p})]\end{array}$$

(42)

where $p_2$ is the average formation pressure in the pressure wave radius at the time of Δt + Δt_f, MPa; C_m is the compression coefficient of shale matrix, MPa⁻¹; C_w is the compression coefficient of water, MPa⁻¹; Z_i is the gas compression factor in the original state, dimensionless; ρ_s is the shale matrix density, kg/m³.

According to the formula of steady-state productivity, the output at the time of Δt + Δt_f can be rewritten as follows:

$${\{\begin{array}{l}{q}_2=\frac{p_2{}^2-p_{\mathrm{wf}}{}^2}{A}\\ A=\frac{\overline{\mu }\overline{Z}{p}_{\mathrm{sc}}T\ln \frac{R_2}{r_{\mathrm{w}}}}{\pi K_{\mathrm{f}}h{T}_{\mathrm{sc}}}\end{array}}_{}$$

(43)

Then we can repeat step (2) to solve for the production at the next production time, and calculate the average formation pressure and the sweep radius under different production times. Furthermore, we can combine the steady-state output calculation formula to obtain the output at different production times. Repeating the above steps, we can obtain the production of shale gas reservoir fracturing wells throughout the production phase. It can be seen that using the SPSS method to calculate the non-steady state production is very convenient, and it is also extremely easy to program. At the same time, it avoids the unreasonable assumption that the pressure has already reached to the boundary of the reservoir under the initial conditions when other mathematical methods solve the productivity model. The different seepage law of gas in volume fracturing zone and shale matrix zone, the adsorption and desorption effect of shale gas are also considered, which is more consistent with the actual situation.

3.3. Model Validation

This model divides the seepage area into volume fracturing area and shale matrix area, taking into account the multi-scale flow and non-linear seepage characteristics of shale gas in horizontal fracturing well. In order to verify the reliability of the SPSS method for solving multi-scale flow problems, a two-region radial seepage numerical model considering shale gas desorption and diffusion is established by using Eclipse [54]. The outer zone of the model is the matrix seepage zone based on the dual media model, and the inner zone is the volume fracturing zone.

By inputting the same model parameters, the cumulative gas production of model in this paper and the Eclipse numerical model are obtained. The two curves have a high degree of conformity, and the output change trend is consistent (Figure 2). The cumulative output of volume fracturing horizontal well composite flow model is slightly higher than that of the Eclipse numerical model. The reason for this is that the SPSS method adopted in this paper avoids the false assumption that the pressure wave spreads to the reservoir boundary under the initial condition. Therefore, it has more desorption gas and describes the actual productivity more accurately.

Figure 3 shows the reservoir pressure distribution of the Eclipse model at 300 days of production. It can be seen that in the shale gas production process, the reservoir pressure gradually decreases from near the wellbore to the reservoir boundary. It reflects that the pressure wave propagates from the bottom of the well to the outer boundary of the volume fracturing zone, and then propagates from the outer boundary of the fracturing zone to the reservoir boundary.

4. Results and Discussion

Due to the complexity of geological features and seepage laws, shale gas productivity is affected by many factors after horizontal well volume fracturing. Based on the combined flow model of fracturing horizontal wells and the parameters of fracturing well in a certain shale gas reservoir (

Table 1. Model calculation basic parameters table.

Parameters (Unit)	Value
Radius of volume fracturing zone (m)	110
Permeability in volume fracturing zone (10−3μm2)	200
Comprehensive compression coefficient of fractured zone (MPa−1)	0.035
Fracture porosity in volume fracturing area (-)	0.1
Supply radius (m)	400
Diffusion coefficient (mm2/s)	281
Comprehensive compression coefficient of shale matrix (MPa−1)	0.019
Compression coefficient of shale matrix (MPa−1)	0.0001
Matrix porosity (-)	0.044
Shale matrix density (kg/m3)	2500
Bottom hole flow pressure (MPa)	5
Viscosity of shale gas (mPa·s)	0.02
Langmuir pressure (MPa)	10
Langmuir volume (m3/kg)	0.05
Original formation pressure (MPa)	30
Formation temperature (K)	360
Borehole radius (m)	0.1
Compressibility of water(MPa−1)	0.0004
Irreducible water saturation (-)	0.1
Reservoir thickness (m)	30
Permeability of matrix (10−3 μm2)	0.005
Pore radius of matrix (nm)	500
Gas slippage constant (-)	−1
Rarefaction effect factor (-)	-

), according to the SPSS method, the effect of different factors such as fracturing zone radius (r_f), Langmuir volume (V_L), the initial average aperture of fractures (b_ei), maximum fracture length (l_max) on gas productivity are studied, and the results are shown in Figure 4, Figure 5, Figure 6 and Figure 7.

Figure 4 shows that with the increase of the radius of complex fracture network in the stimulated area, the daily output will also increase significantly. However, with the increase of production time, the difference of daily output under each fracturing radius is also decreasing.

Figure 5 shows that with the increase of Langmuir volume, daily gas production gradually increased, but the growth rate gradually decreased. This is because the larger Langmuir volume indicates that the reservoir has more adsorbed gas, so in the process of production, more adsorbed gas is desorbed and the free gas is taken out when the pressure is reduced. However, when the Langmuir volume reaches 0.075 m³/kg, the daily output of different Langmuir volumes are almost identical, at this time, the effect of Langmuir volume on the desorption output can almost be ignored.

Figure 6 shows that under the condition of small initial fracture aperture, the influence of initial fracture aperture on daily gas production is significant. The flow of gas in open fractures is the main channel for shale gas transport. The larger the aperture of fractures, the more the flow of gas provided. With the increase of initial fracture aperture, daily gas production gradually increased, but the growth rate gradually decreased.

The fracture length of shale gas well after fracturing is an important index affecting shale gas production. The longer the fracture length, the wider the area affected by volume fracturing. Figure 7 shows that maximum fracture length has the most obvious effect on the initial production of shale gas well. As production continues, the daily output is getting closer.

5. Conclusions

(1)

The productivity prediction model based on the SPSS method provides a theoretical basis for the transient productivity calculation of shale fractured horizontal wells, and it has the characteristics of simple solution process, fast computation speed and high agreement with numerical simulation results.

(2)

The pressure wave propagates from the bottom of the well to the outer boundary of the volume fracturing zone, and then propagates from the outer boundary of the fracturing zone to the reservoir boundary.

(3)

With the increase of fracturing zone radius, the initial average aperture of fractures, maximum fracture length, the productivity of shale gas increases, and the increase rate gradually decreases. When the fracturing zone radius is 150 m, the daily output is approximately twice as much as that of 75 m. If the initial average aperture of fractures is 50 μm, the daily output is about half of that when the initial average aperture is 100 μm. When the maximum fracture length increases from 50 m to 100 m, the daily output only increases about by 25%.

(4)

When the Langmuir volume is relatively large, the daily outputs of different Langmuir volumes are almost identical, and the effect of Langmuir volume on the desorption output can almost be ignored.