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

After volume fracturing, shale reservoirs can be divided into nonlinear seepage areas controlled by microor 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.


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.

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

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].(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 e r , while the gas flow is consistent with Darcy law in stimulated zone with the radius of f r .

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.

Shale Matrix Gas Seepage Model
Knudsen dimensionless number Kn is defined as: where  is the gas molecules' mean free path, nm; rp is the pore radius, nm.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 e while the gas flow is consistent with Darcy law in stimulated zone with the radius of r f .

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.

Shale Matrix Gas Seepage Model
Knudsen dimensionless number K n is defined as: where λ 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: where v is the gas seepage velocity, m/s; K m is the permeability of matrix, 10 −3 µm 2 ; µ 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]: where D k is the diffusion coefficient, mm 2 /s; Z is the gas deviation factor, dimensionless; R g is the universal gas constant, J•mol −1 •K −1 ; T is the shale formation temperature, K; M w is the molar mass of the gas, kg•mol −1 .The permeability in shale matrix is defined as [44]: where φ m is the matrix porosity of shale, dimensionless; τ is the tortuosity, dimensionless (The value is 1).The gas mean free path of a molecule is defined as [28]: 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: If α = 0, b = 1, combining Equations ( 2) and ( 6), differential equation of gas seepage in shale matrix pores can be expressed as follows: The seepage mathematical model in the external matrix area can be expressed as follows:

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.
where k x and k y are the permeability in x direction and y direction respectively, 10 −3 µm 2 ; N is the number of fractures, dimensionless; A is the area of fractured region, mm 2 ; 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]: where k nx and k ny are the permeability in the x and y direction due to normal closure of fractures respectively, 10 −3 µm 2 ; k dx and k dy are the permeability in the x and y direction due to shear failure of fractures, respectively, in 10 −3 µm 2 .
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]: 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.
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]: 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; σ x and σ 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: 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: 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]: 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]: The comprehensive permeability of stimulated fracture network area is obtained through the root mean square of x and y permeability [49]: where K f is the permeability in volume fracturing zone, 10 −3 µm 2 .
According to the hypothesis, the Darcy seepage equation in the internal volume fracturing area can be expressed as follows:

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: where q 1sc is the volume fracturing zone flow, m 3 /d; h is the reservoir thickness, m; T sc is the temperature under standard condition, K; p f is the outer boundary pressure of volume fracturing zone, MPa; p wf is the bottom hole flow pressure, MPa; r f is the radius of volume fracturing zone, m; r w is the borehole radius, m; p sc is the pressure under standard conditions, MPa; µ is the gas viscosity under average formation pressure, m Pa•s; Z is the compression factor under average formation pressure, dimensionless.The flow of external matrix seepage area is obtained from Equation ( 8): where q 2sc is the matrix seepage area flow, m 3 /d; r 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: where q sc is the production of shale gas under the standard condition, m 3 /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: where A =

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: where Q is the desorption gas volume of per kilogram mass shale matrix, m 3 /kg; V L is the Langmuir volume, m 3 /kg; p i is the original formation pressure, MPa; p 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 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:

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 1 from exceeding the outer boundary of the volume fracturing area: where R 1 is the propagation radius of pressure wave at the time of ∆t, m; ∆t is the time of a production time step, d; φ f is the fracture porosity in volume fracturing area, dimensionless; C tf is the comprehensive compression coefficient of fractured zone, MPa −1 .
According to the formula of steady-state productivity, the production at the time of ∆t can be expressed as follows: where q 1 is the production at the time of ∆t, m 3 .

The Solution of Production at the Next Production Time Step
Using q 1 as initial output, the pressure wave propagation distance R 2 at the time of ∆t + ∆t f is obtained: where R 2 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 2 can be expressed as follows: where Gp 2 is the cumulative yield in the time of the ∆t f , m 3 .The free gas geological reserves in the pressure wave propagation radius is calculated by the volume method: 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 3 ; 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 3 .
According to the principle of material balance, under the ground standard condition, material balance equation can be expressed as follows [53]: where G m is the free gas volume in the matrix under the original formation pressure, m 3 ; G f is the free gas volume in the fractures under the original formation pressure, m 3 ; G a is the adsorptive gas volume under the original formation pressure, m 3 ; G p is the produced gas volume, m 3 ; G m is the free gas volume in the matrix under the current formation pressure, m 3 ; G f is the free gas volume in the fractures under the current formation pressure, m 3 ; G a is the adsorptive gas volume under the current formation pressure, m 3 .Equation ( 40) can be rewritten as follows: The equation about the average formation pressure p 2 in the range of pressure propagation radius (R 2 ) can be expressed as follows: 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 −1 ; C w is the compression coefficient of water, MPa −1 ; Z i is the gas compression factor in the original state, dimensionless; ρ s is the shale matrix density, kg/m 3 .According to the formula of steady-state productivity, the output at the time of ∆t + ∆t f can be rewritten as follows: 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.

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.

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), according to the SPSS method, the effect of different factors such as fracturing zone radius (rf), Langmuir volume (VL), the initial average aperture of fractures (bei), maximum fracture length (lmax) on gas productivity are studied, and the results are shown in Figures 4-7.    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.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 3 /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 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 3 /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 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 3 /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.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.

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

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.

Figure 1
Figure 1 Physical model of fractured horizontal well in shale reservoir.

Figure 1 .
Figure 1.Physical model of fractured horizontal well in shale reservoir.

Figure 2 .
Figure 2. Comparison of cumulative gas production of model in this paper and the Eclipse numerical model.

Figure 3 .
Figure 3. Reservoir pressure simulation of the Eclipse model.

Figure 2 .
Figure 2. Comparison of cumulative gas production of model in this paper and the Eclipse numerical model.

Figure 3
Figure3shows 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.

Figure 2 .
Figure 2. Comparison of cumulative gas production of model in this paper and the Eclipse numerical model.

Figure 3 .
Figure 3. Reservoir pressure simulation of the Eclipse model.

Figure 3 .
Figure 3. Reservoir pressure simulation of the Eclipse model.

Energies 2018 , 17 Figure 5 .
Figure 5.The influence of Langmuir volume on the daily gas production.

Figure 6 .
Figure 6.The influence of initial average aperture of fractures on the daily gas production.

Figure 5 .
Figure 5.The influence of Langmuir volume on the daily gas production.

Energies 2018 , 17 Figure 5 .
Figure 5.The influence of Langmuir volume on the daily gas production.

Figure 6 .
Figure 6.The influence of initial average aperture of fractures on the daily gas production.

Figure 6 .
Figure 6.The influence of initial average aperture of fractures on the daily gas production.

Figure 7 .
Figure 7.The influence of maximum fracture length on the daily gas production.

Figure 7 .
Figure 7.The influence of maximum fracture length on the daily gas production.