Bottom-hole pressure drawdown management of fractured horizontal wells in shale gas reservoirs using a semi-analytical model

Due to strong stress sensitivity resulted from unconventional tight formationsit is of practical interest to formulate a reasonable pressure drawdown plan to improve gas extraction recovery. The impact of water-shale interactions on the reservoir permeability was previously ignored in the managed pressure drawdown optimization. The controlled-pressure production dynamic analysis was mostly conducted using numerical simulation, lack of rigorous theoretical support. Hence in this paper, a theoretical production prediction model was proposed and verified with HIS RTA 2015by incorporating multiple pressure drawdown mechanisms and various non-linear gas flow process. The on-site production effects dominated by two different pressure drop methods was further compared, indicating that compared to depressurization production, the production reversion can occur in the controlled pressure production process and the EUR of single well can be increased by about 30% under the control of managed pressure drawdown approach. Finally, the pressure drawdown optimization strategy was carried out on the field test from the both production effect and economic benefits, which demonstrated that the best economic solution can generally be obtained in the early stage of production. The research results can be closely linked to the on-site production practice of shale gas wells, providing insights into designing optimized production strategy scheme.

Since 2010, three industrial production demonstration zones in the marine shale gas reservoirs have been established in China, making China the third country in the world to realize the commercial development of shale gas reservoirs 1 . Deep shale gas resources will become the main body of natural gas production in the future. By 2040, China is likely to be the largest shale gas producing area after North America 2 . Shale gas has become an important replacement resource to make up for the shortage of conventional energy 2,3 , and it has affected profoundly the world's energy production and consumption patterns.
Shale gas reservoir is typically characterized by low-porosity and extremely low-permeability. The permeability of the reservoir with poor pore connectivity is between 10 -5 ~ 10 −3 mD 4-6 . The large resistance in gas seepage of shale gas reservoir makes no productivity under natural conditions. The horizontal well and hydraulic fracturing technology are generally adopted to achieve industrialized gas flow. The rapid decline rate of initial daily gas production for shale gas wells can be attributed to the strong stress sensitivity. During the production process of a shale gas well, the reservoir conductive medium is gradually compacted with the increase of the closure pressure, and the shale gas flow channel is deformed, resulting in a decrease in seepage capacity.
In the practice of shale gas reservoir development, it is necessary to clarify a suitable bottom hole pressure drawdown strategy, which is conducive to reducing the daily gas production decline rate and increasing the ultimate productivity of shale gas single well. Generally, the production methos of shale gas wells can be divided into depressurization production and controlled pressure production. Depressurization production means that there is no flow restriction at the wellhead, and the bottom hole flowing pressure at the initial production stage of gas well is rapidly reduced to the constant pressure. The shale gas wells with depressurization production can www.nature.com/scientificreports/ channel width of the fracture is reduced, and the fracture conductivity can be impaired 29 . Pressure-controlled production can appropriately adjust the bottom hole pressure decline rate, control the increase in effective stress, and effectively avoid the reduction of fracture conductivity, to increase the long-term EUR of the gas well.
Assuming that the impact of proppant elastic embedding and deformation of proppant on the width of hydraulic fractures is linearly related to the net confining pressure 18 , the expressions (1) and (2) are expressed as follows: Width loss caused by elastic embedding of proppant: Width loss caused by elastic fracture of proppant: The effective width of the proppant-fracture is: where p avg is average pore pressure, Pa; v m is reservoir Poisson's ratio, dimensionless; D p is median value of proppant particle size, m; E p is Young's elastic modulus of proppant, Pa; w 0 is initial width of hydraulic fracture, m; v p is proppant Poisson's ratio, dimensionless.
Microfracture water-rock interaction. Compared with matrix and proppant hydraulic fractures, unsupported microfractures are more sensitive to net reservoir stress. The invasion of excess fracturing fluid into the fractures near the well hole bore can make the fractures surface hydration 30,31 and osmotic hydration 30,31 , successively causing clay minerals to expand 32 , weakening the mechanical properties of the reservoir 33 , and enhancing the stress sensitivity (Figs. 2 and 3). Even after the stress is restored, the permeability can only recover to 10% to 15% of the initial permeability 34 . In this work, the permeability is assumed to decrease exponentially with the increase of effective stress 35 . The stress sensitivity coefficient after fracturing fluid soaking is not constant, but is negatively correlated with effective stress 9 . The permeability in the water-bearing equivalent fracture network is: where:   www.nature.com/scientificreports/ where K f is facture permeability considering stress sensitivity, m 2 ; K fi is intrinsic fracture permeability, m 2 ; γ f is stress sensitivity coefficient of the fracture, Pa -1 ; p e is initial formation pressure, Pa; p f is fracture pore pressure, Pa.
Matrix stress sensitivity. The increase effective stress on the rock skeleton of the reservoir results in obvious elastoplastic deformation of the rock pore structure (Fig. 4), and thus matrix permeability, porosity, and rock physical parameters have changed. Proper pressure-controlled production can reduce the shrinkage rate of rock pore volume and improve the connectivity between effective pores. Therefore, an exponential stress sensitivity empirical model 35 was introduced in the paper to characterize the influence of matrix stress sensitivity on matrix seepage capacity: where K m is matrix permeability considering stress sensitivity, m 2 ; K mi is intrinsic matrix permeability, m 2 ; γ m is stress sensitivity coefficient of the fracture, Pa -1 ; p m is matrix pore pressure, Pa.

Methodology
Physical model description. The diversity of pore structure scales in shale reservoirs makes the shale gas flow process complicated, and the flow path covers the molecular scale to the macro scale 36 . The shale gas production period is generally attributed to multiple migration mechanisms across scales, and it is mainly divided into three stages ( Fig. 5): In the first stage, free gas in the proppant fractures and unproppant fractures moves towards the wellbore by pressure drop near the wellbore; During the second stage, free gas in the matrix flow to natural fractures, motivated by the pressure difference between the matrix and the fracture; In the third stage, the reservoir pressure change in the matrix pores can promote microscopic gas flow such as the diffusion and slippage of free gas and the surface diffusion of adsorbed gas to enhance gas flow capacity in the matrix. When the matrix pore pressure is declined to the critical desorption pressure, the adsorbed gas starts to mobilize into free gas in the matrix pores. Irregular fractures distribution and strong heterogeneity in the fractured reservoir make it difficult to acquire the specific fracture and certain reservoir parameters. Thus, the fracture and matrix can be simplified to facilitate the calculation process. It is assumed that all the hydraulic fractures can be characterized by bi wing transverse 37 , note that the increased effective stress may lead to elastic embedding and deformation of proppant in hydraulic fractures. Furthermore, the secondary fracture network and matrix are considered as equivalent medium. Considering that there are some unstimulated reservoir zones between the hydraulic fractures and outside the fracture tip, the proposed model in the paper can include the following five flow areas: hydraulic fractures, fracture network area 1, matrix area 2 between hydraulic fractures, and unstimulated matrix area 3 and matrix area 4, as is shown in the Fig. 6. The single-phase methane gas flow behavior in different areas is one-dimensional flow. The stress sensitivity, high-pressure adsorption, diffusion and viscous flow are all considered in the matrix 2, 3, and 4. Moreover, the hydration expansion of unproppant fracture network can aggravate the clay swelling, which is the focus of pressure control protection for gas wells.
Where ψ e is initial formation pseudo-pressure, Pa 2 /(Pa • s); ψ is the formation pseudo-pressure, Pa2/(Pa • s); K F is hydraulic fracture permeability, m 2 ; ϕ 2 is porosity of the matrix zone 2, ,dimensionless; ϕ 1 is porosity of fracture network, dimensionless; ϕ 3 is porosity of matrix zone 3. m 2 ; C t3 is comprehensive compressibility of matrix zone 3, Pa -1 ; C t2 is comprehensive compressibility of matrix zone 2, Pa -1 ; C t1 is comprehensive compressibility of fracture network, Pa -1 ; t a is pseudo-time.s; L f is hydraulic fracture spacing, m; A cw is wellbore crossflow area,m 2 ; ϕ F is porosity of hydraulic fracture, m; C tF is comprehensive compressibility of hydraulic fracture,Pa -1 ;  where ϕ 4 is porosity of matrix zone 4, dimensionless; C t4 is the matrix comprehensive compressibility coefficient, Pa -1 ; P 4 is pore pressure in matrix zone 4, Pa; t is production duration, s; ρ g is gas density,,kg/m 3 ; K 4a is apparent permeability of matrix zone 4, m 2 ; µ 4 is gas viscosity, Pa • s. where: the matrix comprehensive compression factor is: Gas compression factor: where C g4 is gas compressibility, Pa −1 ; C d4 is modified supercritical desorption gas compression coefficient of matrix zone 4, Pa −1 ; C f4 is formation compressibility, Pa −1 .
Considering that the conventional Langmuir equation cannot be applied to high pressure isotherm adsorption process (Fig. 7), the high pressure isotherm adsorption model 38 is adopted and desorption gas compression factor 39 is transformed as Eq. (10): where V L is Langmuir volume, m 3 /m 3 ; ϕ 4 is porosity of the matrix zone 4, dimensionless; T is formation temperature, K; Z is gas compression factor, dimensionless; P sc is standard atmospheric pressure, Pa; Z sc is standard gas compression factor, taking the value of 1, dimensionless; T sc is standard atmospheric temperature, K; P L is Langmuir pressure, Pa; ρ a is absorbed gas density, kg/m 3 .
In this paper, the stress sensitivity of the matrix and the desorption of gas in the matrix are coupled, and the apparent permeability model of gas flow in the shale matrix micro/nanopores based on molecular dynamics theory was used to superimpose the viscous flow and Knudsen diffusion of real gas molecules with weighting coefficients: where: effective Knudsen number: www.nature.com/scientificreports/ Considering stress sensitivity and desorption process the effective hydraulic flow radius: Surface coverage of adsorbed gas on matrix pore wall: where K 4a is apparent permeability of matrix zone 4, m 2 ; K 4i is intrinsic permeability of matrix zone 4, m 2 ; D is gas diffusion coefficient, m 2 /s; is molecular free path of methane molecule, m; r e is effective hydraulic radius, m; r 0 is original hydraulic radius, m; K B is Boltzmann constant, 1.38065*10-23 J/K; d is methane molecular collision diameter,m; p is average pore pressure in the circular pore, m; d CH 4 is methane molecular collision diameter, m. The pseudo pressure and pseudo-time 37 are adopted to linearize the high-pressure physical property parameters in the flowing control equation to solve the equations easily.
The pseudo-pressure is: The pseudo-time is: where µ i is gas viscosity at initial formation pressure condition, Pa • s ; C ti is comprehensive compression factor of matrix zone 4 at initial formation pressure condition, Pa -1 ; µ(p) is gas viscosity at formation pressure of p , Pa • s ; C ti is comprehensive compression factor of matrix zone 4 at formation pressure of p ,, Pa -1 ; According to the definition of pseudo-pressure and pseudo-time and the dimensionless parameters definition Table 1, Eq. (7) can be rewritten as follows: Boundary conditions can be expressed in dimensionless form are as follows: Governing equation for gas flow in the matrix of zone 3. Considering the reservoir stress sensitivity, gas desorption, diffusion and viscous flow from matrix zone 3 into zone 1 along the y direction, the dimensionless gas percolation equation in the matrix of zone 3 can be derived in Eq. (21):

Dimensionless parameter Definition Dimensionless parameter Definition
Dimensionless pseudo-pressure Dimensionless formation conductivity www.nature.com/scientificreports/ As the outer boundary condition is closed and the pressure on the inner boundary is continuous, the dimensionless initial and boundary conditions are reflected in Eqs. (22)(23)(24) as follows: Governing equation for gas flow in matrix of zone 2 40 . Considering the unsteady gas flow exchange in matrix zone 4 and matrix zone 2, the dimensionless gas seepage equation in the matrix of zone 2 is established in Eq. (25): The outer boundary is closed, the inner boundary pressure is continuous, the initial and boundary conditions are expressed in Eqs. (26)(27)(28): The initial condition: Outer boundary: Inner boundary: Governing equation for gas flow in fracture network of zone 1. Consider the influence of hydration on the seepage capacity of the equivalent fracture network in zone 1, the flow rate on the outer boundary between zone 1 and zone 2 is continuous, and the inner boundary pressure is continuous. The dimensionless gas flow equation in fracture network zone 1 can be depicted in expression (29): The boundary conditions are as followed in Eqs. (30)(31)(32): Initial condition: Outer boundary: Inner boundary: Governing equation for gas flow in inner zone hydraulic fracture. The gas mass transfer in the hydraulic fracture is mainly dominated by viscous flow, and the elastic embedding and deformation of proppant in the hydraulic fractures cannot be ignored. The dimensionless flow control Eq. (33) in the hydraulic fracture is established: where w F is hydraulic fracture width,m.Considering the constant bottom hole pressure of the gas well and the outer boundary is closed, the initial and boundary conditions are as followed: The initial condition: The Heaviside (x) function 42 was introduced to transform Eq. (37) into a continuity function (38): The dimensionless bottom hole pseudo pressure ψ wfD (t aD ) was transformed in the Laplace space as the expression (39) : The second integral term on the right side of the above formula (39) can be simplified to obtain expression (40): Discretize Eq. (40), as shown in Fig. 8: where: www.nature.com/scientificreports/ Next, the dimensionless seepage equations can be transformed in Laplace domain and the final solution was acquired. Then the semi-analytical solution of the dimensionless production in the real space is acquired with the Stehfest numerical inversion 43 . By the Newton iteration method, the production solution in the real space at constant pressure can be derived and the specific solution process was shown in the Fig. 9. As the high lights of this paper are theoretical analysis of dynamic production performance of gas wells with managed pressure drawdown and optimization of pressure drop strategies rather than model derivation, the solution of the model is directly given in this paper. The specific solution derivation details have been illustrated 44 .
The total dimensionless production rate at the bottom hole of the shale gas well in Laplace space can be derived in Eq. (42): where N is the number of hydraulic fractures, dimensionless.

Model verification and analysis
As is shown in the Fig. 10, a multi-zone composite linear flow numerical model by the IHS simulator was adopted to verify the proposed model with the same model input parameters in the Table 2. The proposed productivity model was simplified due to that the only Darcy flow is considered in the simulator, mainly ignoring the non-linear gas flow mechanism, supercritical desorption and reservoir stress sensitivity, while constant stress sensitivity of the matrix and fracture, high-pressure physical properties (see Figs. 11 and 12), and Langmuir desorption were all considered. The bottom hole pressure drop path in the numerical simulation is assumed to be stepped and the pressure control period is 3 years (Fig. 13), the production prediction results of the two models were obtained, as shown in Fig. 14 www.nature.com/scientificreports/ method to conduct actual productivity prediction and dynamic performance analysis for shale gas wells with managed pressure drawdown. Figures 15, 16, 17 show that when the non-linear gas flow mechanism is not considered in the production model, the peak daily gas production is low, and the daily gas production decline significantly. The production capacity predicted by the model can be underestimated by 30.3%, indicating diffusion and desorption is significant to the later-staged-productivity; the calculated EUR of considering variable stress sensitivity is 227.7 million cubic meters, 11.8% higher than that calculated by the constant stress sensitivity productivity model. The increase in net stress leads to the reservoir compaction, and the reduced reservoir seepage capacity; the model considering the reservoir water saturation of 45% predicts that EUR is 39.4% lower than that of dry reservoir. It can be seen that the fracturing fluid indeed caused damage to reservoir seepage capacity; when the gas well is put into production, the hydraulic fracture proppant can suffer elastic embedding and elastic deformation to result in 5.75% loss of the EUR. If the above comprehensive factors are considered in the model, the final EUR is about       www.nature.com/scientificreports/ 35.5% lower than that of the weakened model. To sum up, the comprehensive factors considered in the production model is conducive to improve the accuracy of predicting the medium and long-term productivity. It aims to provide theoretical reference and engineering value for the optimization analysis of pressure reduction strategies.

Case study
Production effect analysis between managed drawdown and high drawdown. Two shale gas wells with different production systems, HD (High Drawdown) and MD (Managed Drawdown), were selected as production effects comparison for basically same testing production, the number of fracturing sections, and geological background. HD well adopts short-term pressure control method for production, while long-term pressure control strategy is applied to MD well. The constant pressure productivity model was used to historically fit the production data of the HD well (Fig. 18) and the pressure-controlled productivity model in the paper was used to historically fit the MD gas well production data (Fig. 19) to invert the unknown geological and engineering parameters of the two wells, and then respectively forecast the EUR for two wells. As is shown in Fig. 20, the short-term-cumulative gas production of managed pressure drawdown is lower than that of high pressure drawdown. The advantages of managed drawdown production in the later stage begin to appear, and "capacity reversal" occurs. The ultimate EUR calculated by the proposed model is about 33.35% higher than that of the constant pressure productivity model. Therefore, the managed pressure drawdown production system can achieve the production goal of increasing the medium and long-term cumulative production of a single well. As can be seen in the Fig. 21, the annual average daily gas production decline rate of gas wells with the managed pressre drawdown system is generally lower than that of gas wells with the high pressre drawdown system. Starting from the second year of production, the annual production decline rate of pressure-controlled gas wells has been reduced from 36.02 to 5.12%, and has been relatively low from the 15th year, basically entering the stable production stage. The annual daily production decline rate of depressurized gas wells has been reduced from 50.80 to 10.55%, which is still in a period of obvious decline in production. It can be seen that pressure-controlled  www.nature.com/scientificreports/ production can control the degree of production decline, protect the seepage capacity of the reservoir, extend the stable production period of gas wells, and achieve high-efficiency increase in ultimate production of gas wells.
Optimization analysis of pressure drawdown management strategy. The pressure drop strategy of an example well in the Sichuan shale gas area was optimized from the production effect and economic benefit evaluation with the proposed productivity model, including the pressure control duration, the initial controlled pressure, and different pressure drop paths, respectively. The geological parameters and engineering of this example well are shown in Table 3.
The impact of pressure control duration on production effects. The stepped pressure drop path can be approximated as a linear pressure drop path during the long-term production period of a gas well, thus the limited pressure drop path is linear. The influence of different pressure control durations of the gas well on the production effect was studied, and the optimal pressure control duration was clarified for the example well. As is shown in Figs. 22,23,24,25,26, comparing the production effect of depressurization production and pressure control for 0.1 to 3 years, with the increase of pressure control duration, the annual average daily gas production of gas wells decreases gradually, and the stable production period lasts for a long time. There is an overall trend of gas well EUR first increasing and then decreasing when the pressure control duration increases. When the pressure is controlled for 3 years, the EUR will show a downward trend, indicating that the pressure control period has an optimal value of 2 years, to obtain the most ideal productivity. Reasonable pressure control production for 2 years can make single EUR be increased by 10.02%. On the contrary, too long pressure control time is not conducive to gas well production performance.  www.nature.com/scientificreports/ The influence of the initial controlled pressure on the production effect. In the early stage of gas well production, well test work is usually carried out to explain the reservoir properties and parameters, and pressure control will not be implemented from the beginning of production process. Thus, the gas well is set to control pressure from different bottom hole pressures starting point to study the effect of different initial controlled pressures on gas well production. Comparing the production performance charts of gas wells with different initial controlled pressures for 2 years of pressure control (Fig. 27, 28, 29, 30, 31), it can be seen that when the initial controlled Table 3. Geological parameters and engineering parameters of an example well.   www.nature.com/scientificreports/ pressure is closer to the original reservoir pressure, that is, the pressure is controlled at the beginning of production process for shale gas wells, low annual average daily gas production decline rate and long stable production period can be obtained. Specifically, the EUR of initial reservoir pressure relative to that of 1/5 of initial reservoir pressure be increased by 8.54%. Conversely, the initial controlled pressure is close to delivery pressure, and it is mainly dominated by depressurization production. It can illustrate that the fast pressure drop in the reservoir near the wellbore will lead to high initial gas production of the gas well, and slow reservoir production rate in   www.nature.com/scientificreports/ the later stage of production, which is attributed to permeability loss and thus the low long-term production. Therefore, the formation pressure drops slowly when the pressure is controlled in the early stage of production, and the long-term cumulative gas production of gas wells is greater, on the contrary, the cumulative gas production of depressurization is the lowest.
Influence of pressure drop path type on production effect. Different reservoir stress sensitivity at different pressure drop paths can lead to obvious distinction in reservoir seepage capacity, which affects the final dynamic pro-   www.nature.com/scientificreports/ duction performance. Thus, the production effect of gas well under different pressure drop paths was evaluated (Fig. 32, 33,34,35): The production effect of a single well is the most ideal when bottom hole flowing pressure decreases linearly, that is, a robust drawdown. When the pressure drop occurs as a conservative drawdown, too low pressure drop rate is not conducive to increasing the long-term production of gas reserves in the matrix. On the contrary, the prominent concave effect of pressure drop path means great production pressure difference near the well, and can cause serious damage to the fractured reservoir conductivity and the boundary production utilization of reservoir reserves.   www.nature.com/scientificreports/ The influence of managed pressure drawdown parameters on economic benefits. The final objective function for optimizing the economic benefits of shale pressure controlled gas wells is determined by maximizing cumulative gas production and minimizing production costs, that is, the 20-year net present value (NPV), and an economic evaluation model is introduced 46 :   www.nature.com/scientificreports/ where i r is annual interest rate,%; m is gas price, yuan/m 3 . Assuming that the investment cost of a single horizontal well in this example is 68 million yuan, the gas price is 1.275 yuan/m 3 , and the annual interest rate is 10%, study the impact of different pressure control durations, pressure drop paths, and initial controlled pressures on net present value to seek the optimal combination of pressure drop parameters and obtain the greatest economic return, and formulate the best economic pressure control plan.
Comparing the NPV chart of 0.1~3-year-pressure controlled production and depressurization production (Fig. 36), the ultimate NPV value of 1-year-managed pressure production is the highest, 17.56% higher than that of depressurization production; Fig. 37 reveals the 20-year NPV value is positively correlated with the initial controlled pressure. The max NPV value can be up to 11.12%; As is shown in Fig. 38, the simulated deep gas well is produced with different pressure drop paths from commissioning to abandonment (20 years of production). The 20-year-NPV of linear pressure drop path is better than the conservative path and the aggressive path. , Which can be increased by 5.46% and 1.24% respectively.

Summary and conclusions
(1) The importance of water-rock interaction of the unproppant fractures to reservoir seepage capacity and long-term production was emphasized in pressure control mechanism. A five-zone composite pressure control productivity model was proposed with comprehensively considering the multiple pressure control mechanisms, and then verified to be reliable with numerical simulator. Then it was proved that more comprehensive gas flow mechanisms and pressure control mechanisms considered in the proposed model can effectively promote higher prediction accuracy in the medium and long-term productivity of the gas well, which aims to provide theoretical basis for the optimization analysis of pressure drop strategy. (2) The production effect between the managed pressure drawdown and high pressure drawdown was compared and analyzed by the proposed model. It was found that the annual average daily gas production decline rate by managed drawdown production was generally lower than that of high drawdown production.Theshale gas well controlled by managed drawdown scheme has a relatively earlier stable-production stage; Compared to depressurization production, the production reversion can occur in the controlled pressure production process and the EUR of single well can be increased by about 30%. Thus, for strongsensitive-shale formations, pressure-controlled production can alleviate the production decline rate, and finally achieve the production goal of increasing the medium and long-term cumulative production. (3) The optimization analysis of pressure drawdown strategy on productivity was carried out.There is an optimal value for the pressure control time of a specific gas well. In this work, two-year-pressure control production can increase 10.02% of the single EUR. On the contrary, too long pressure control time is not conducive to gas well production; when pressure control begins as soon as a gas well is put into production, the productivity can reach the peak and thus the EUR can be increased by 8.54% compared with the depressurization production. The aggressive or conservative bottom hole pressure drop path is not conducive to the long-term production of a single well, while the steady pressure drop path has the least conductivity loss of gas wells to obtain the most ideal reservoir utilization. (4) The economic evaluation model was introduced to formulate the optimal economic pressure drawdown strategy. Generally speaking, the best economic solution can generally be obtained in the initial stage of production.The optimal economic pressure control duration is at an inflection point. The economic benefits first increase and then decrease with the increase in the pressure control duration; Furthermore, the 20-year-NPV value of gas wells is positively correlated with the initial controlled pressure; The NPV of the gas well with the linear pressure drop path is better than the conservative path and the aggressive path. To sum up, the research results can be closely linked to the on-site production practice of shale gas wells, and is of theoretical reference significance and engineering application value for formulating the optimized production strategy plan to improve the ultimate recovery rate of a single well.

Data availability
The datasets generated and/or analyzed during the current study are not publicly available due to the fact that the data forms part of current ongoing project but are available from the corresponding author on reasonable request.