Transient Pressure Analysis for Multifractured Horizontal Well with the Use of Multilinear Flow Model in Shale Gas Reservoir

Shale gas reservoirs (SGR) have been a central supply of carbon hydrogen energy consumption and hence widely produced with the assistance of advanced hydraulic fracturing technologies. On the one hand, due to the inherent ultralow permeability and porosity, there is stress sensitivity in the reservoirs generally. On the other hand, hydraulic fractures and the stimulated reservoir volume (SRV) generated by the massive hydraulic fracturing operation have contrast properties with the original reservoirs. These two phenomena pose huge challenges in SGR transient pressure analysis. Limited works have been done to take the stress sensitivity and spatially varying permeability of the SRV zone into consideration simultaneously. This paper ﬁ rst idealizes the SGR to be four linear composite regions. What is more, the SRV zone is further divided into subsections on the basis of nonuniform distribution of proppant within the SRV zone which easily yields spatially varying permeability away from the main hydraulic fracture. By means of perturbation transformation and Laplace transformation, an analytical multilinear ﬂ ow model (MLFM) is obtained and validated as a comparison with the previous models. The ﬂ ow regimes are identi ﬁ ed, and the sensitivity analysis of critical parameters is conducted to further understand the transient pressure behaviors. The research results provided by this work are of signi ﬁ cance for an e ﬀ ective recovery of SGR resources.


Introduction
Due to extremely low permeability and porosity of shale gas reservoir (SGR), multistage hydraulic fracturing has become an integral tool to improve the gas recovery.The economic feasibility of shale gas reservoirs has a strong relationship with the fracture system permeability near the wellbore.Considered to be the most effective way to produce gas resources, a multistage fractured horizontal well can create several highconductivity hydraulic fractures as flow paths and, at the same time, activate and connect existing natural fractures so as to develop a large fracture network system [1].The zone containing the main high-conductivity hydraulic fractures and large spatial network system which can effectively improve well performance is defined as SRV (stimulated res-ervoir volume), and the remaining zone which is hardly influenced by the treatment of hydraulic fracturing is similarly defined as USRV (unstimulated reservoir volume) (Mayerhofer et al. 2006, [2,3]).
The presence of a complex fracture network in the SRV has a significant impact on the pressure transient analysis of unconventional reservoirs.Analytical and semianalytical approaches have been used to model the transient flow behavior in such systems.Zhao et al. [4,5] and Wang [6] have established semianalytical solutions with the use of Laplace transformation.The point source function or line source function, coupled with superposition principle, was utilized to mathematically incorporate the interference among hydraulic fractures.Alternatively, multilinear flow modes have been extensively developed to simulate the gas flow in SGR.The SGR is generally divided into some coupled zones and linear flow patterns are assumed.These models also assumed that continuous pressure drops along the hydraulic fractures exist to push the hydrogen gas to the wellbore.El-Banbi [7], Al-Ahmadi and Wattenbarger [8], Xu et al. [9], Stalgorova and Mattar [10] simplified SGR as linear composition reservoirs, and the governing equations of each zone can be derived.
As we all know, after the hydraulic fracturing stimulation of shale gas reservoirs, the induced fracture and frac formation can be subdivided into several categories based on the fracturing pattern and proppant distribution.Therefore, it is too idealistic to simply assume that the induced hydraulic fracture is bi-wing.In order to concisely describe the fracture network (natural or induced) around the hydraulic fractures, some aforementioned methods (including analytical methods, semianalytical methods, and even time-consumption numerical methods) are proposed.While assuming uniform distribution of identical hydraulic fractures along the length of the horizontal well, Ozkan and Koseler [11] and Ozkan et al. [2] utilize the concept of a trilinear model with a naturally fractured inner reservoir to represent the MFH well performance in unconventional reservoirs.Brown et al. [12] presented an analytical trilinear flow model that incorporates transient fluid transfer from the matrix to the fracture to simulate the pressure transient and production behavior of fractured horizontal wells in unconventional reservoirs.However, those proposed analytical models lack the ability to explicitly consider the medium conductivity secondary fracture, which absolutely induces certain errors to some degree.And then, Zeng et al. [13], Stalgorova and Mattar [10], and Wang et al. [14] used an unstructured-grid simulator to analyze the type curves of an MFHW with fracture complexity; the generation of a complex unstructured grid is time consuming; moreover, the adjacency of the high-conductivity fractures is refined with fine grids making it necessary to increase the complication and economical consumption of computation, while this latest technology can make the simulation of fracture complexity more accurate.
The aforementioned analytical models significantly simplify the geometry of a complex fracture system, e.g., regular orthogonal hydraulic fractures, and stress sensitivity phenomenon.To resolve these shortages, plenty of semianalytical models without the above simplifications have been established.Zhao et al. [5] proposed a radial compound model which treated the SRV as a circular area with high permeability.This model was used to simulate the performance of multifractured horizontal wells in SGR.According to source functions, Jiang et al. [15] analyzed pressure and gas rate transient laws of a multistage fractured horizontal well for tight oil reservoirs while considering SRV, although the fractures were still confined to be vertical to the wellbore and the stress sensitivity was also ignored in their study.Zongxiao et al. [16,17] amalgamated the perturbation technique with a linear source function method to consider the effect of stress sensitivity when they researched transient pressure behaviors of horizontal wells.However, the com-plexity of fracture networks was still neglected.Jia et al. [18] proposed a new semianalytical model by combining finite-difference and line source functions to study the flow behaviors of a horizontal well after hydraulic fracturing; unfortunately, the stress sensitivity effect was still ignored.Given the stress sensitivity effect of fractures and reservoirs, Wang et al. [19] established a semianalytical model suitable to study transient flow behaviors of a multistage fractured horizontal well with complex fracture networks without considering the property contrast between the stimulated zone and the unstimulated zone.
Recently, some works tend to consider the well interference using semianalytical models, such as Xiao et al. [20] and Jia et al. [21].
Although a semianalytical radial flow model could be employed to characterize hydraulic fractures with complex topology, the multilinear flow models are easy to be used in the real-field application.By means of perturbation transformation and Laplace transformation, we proposed a new analytical multilinear flow model to systematically investigate the effects of stress sensitivity in SGR.In addition, the SRV zone is further divided into subsections on the basis of nonuniform distribution of proppant within the SRV zone which easily yields spatially varying permeability away from the main hydraulic fracture.This paper also discusses the influence of relevant parameters on the of fractured horizontal wells in stress-sensitive SGR, including stress sensitivity, mobility ratio of the SRV and the outer region, SRV size, and coefficient of permeability variation.Corresponding solutions can be useful for fracture design and well test interpretation in field practice.

Mathematical Model and Analytical Solution
2.1.Model Descriptions.As illustrated in Figure 1(a), the microseismic data maps show that hydraulic fracturing treatments create irregular fracture geometry.A simplified physical model of a multifractured horizontal well in SGR is illustrated in Figure 1(b).The entire reservoir after hydraulic fracturing can be conceptually divided into four coupled linear flow zones with different properties, including the main hydraulic fracture Zone I, SRV Zone II, SRV Zone III, and the other outer zone without stimulation Zone IV.The fracture half-length is x f and the width and length of the entire zone are x e and y e , respectively.The length of SRV Zone II is l.All these variables have been depicted in Figure 1(b).The initial permeability of these four zones are separately denoted as k 1 , k 2 , k 3 , and k 4 .
As far as we know, the nonuniform distribution of proppant within the SRV zone easily yields spatially varying permeability away from the main hydraulic fracture.In this paper, we further divide Zone III into N small zones as illustrated in Figure 1(c).All zones are idealized as dual-porosity media with different types of permeability.The stress sensitivity of permeability is taken into consideration.Singlephase and microcompressible gas is assumed to flow in SGR, which follows Darcy's Law.  3 Geofluids 2.2.Mathematical Formula of Multilinear Model.To begin with, we define the pseudopressure as follows: and then, According to the theory proposed by Pedrosa Jr. [22], stress-sensitive permeability could be described as follows: where the reference permeability k 0 is the initial permeability in SGR, m 2 , and γ i is the permeability modulus, which is generally obtained by means of indoor experiments.We also assume that four zones have different permeability moduli.The mechanism of adsorption can be classified into instant adsorption and time-dependent adsorption, and the former is selected to describe adsorption phenomenon in this paper.At present, the Langmuir isotherm equation is used to describe the instant adsorption process of shale gas [23], and its expression is where V E is the adsorption equilibrium concentration, sm 3 /m 3 ; V L is the Langmuir adsorption concentration, sm 3 /m 3 ; P L is the Langmuir pressure, MPa; and m L is the Langmuir pressure, MPa 2 /(mPa•s).
For the convenience and simplicity of deducing formulas, some dimensionless parameters are introduced first: A multilinear flow model (MLFM) will be used to derive the mathematical equations.In the following parts, the governing equations are separately established for each zone.

Zone I.
In this zone, the gas is supplied from the adjacent reservoir Zone II to the main hydraulic fractures and then flows to the wellbore with a linear flow pattern.Gas pseudopressure is used to consider the effects of gas compressibility.When the stress sensitivity is considered, the governing function can be presented as follows: with the inner and outer boundary conditions: Substituting predefined dimensionless variables to Equations ( 6) and (7), their dimensionless formula is as follows: Equation ( 8) is a strongly nonlinear partial differential equation, which is not convenient to be solved analytically.

Geofluids
A perturbation transformation proposed by Pedrosa Jr. [22] can be used to eliminate the nonlinearity.The new dimensionless variables η jD related to the dimensionless pressure are introduced as follows: According to the simplified method proposed by [6], due to the dimensionless permeability modulus γ D which is usually a small value, η jD can be expanded as a power series in the parameter γ D .
Substituting Equations (A.6) and (A.7) into Equation ( 8), we can get a sequence of linear problem that can be solved for η jD0 , η jD1 , and so on.According to Liu et al. [24] the zero-order approximation η jD0 was accurate enough for pressure analysis.
Finally, Laplace transformation can be used to transform these equations from time domain to Laplace domain so as to eliminate the effects of the time domain.

Zone II.
Zone II is a highly permeable SVR zone which is adjacent to the main hydraulic fracture as illustrated in Figure 1(b).This area is significantly stimulated by the massive hydraulic fracturing operation; as a result, it can be assumed to single porous media with high permeability due to the support by the proppant.When the stress sensitivity is considered, the governing function can be presented as follows: The left condition is adjacent to Zone I, while the right condition is connected to Zone III.Specifically speaking, these two conditions can be presented as follows: Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows: with the following boundary conditions: Similarly, Laplace transformation can be used to transform these equations from the time domain to the Laplace domain.
with the following boundary conditions: 2.2.3.Zone III.In this work, Zone III is also assumed to be a stimulated area induced by massive hydraulic fracturing operation.Unlike Zone II, this zone is partially supported by the proppant.In addition, the preexisting natural fractures are stimulated as well.On condition to these facts, Zone III is idealized as a dual-porosity media where the fracture and matrix are coupled.As we have mentioned above, to systematically investigate the influences of spatially varying permeability away from the main hydraulic fracture, we have divided Zone III into N small zones as illustrated in Figure 1(b).We need to establish flow equations for each small zones and couple them together on the basis of continuity condition at the interface of neighboring zones.Here, we take the j-th small zone as an example to derive a set of generic governing equations.We also should separately derive the governing equations for the fracture system and matrix system.The effects of adsorption and pseudo-steady state interporosity flow from the matrix to the micro fracture system are considered simultaneously.When the stress sensitivity is 5 Geofluids considered, the governing function can be presented as follows: Fracture system: The boundary conditions are derived from the interface of neighboring two sections, where the pressure and gas flux should be strictly equal between those two sections.
In particular, the first small zone will directly connect to Zone II, and the right side of the N-th small zone is considered to be impermeable.The boundary conditions could be explicitly presented as follows: Obtained from Fick's first law of diffusion, the diffusion rate of shale gas can be expressed as Substituting the Langmuir isotherm Equation (4) into Equation ( 23), we can obtain Substituting predefined dimensionless variables to Equations ( 19)-( 24) could obtain their dimensionless formulas, and then the aforementioned perturbation transformation proposed by Pedrosa Jr. [22] is used to eliminate the nonlinearity.Finally, their dimensionless formulas with the consideration of stress sensitivity effects can be as follows: Fracture system: with the following boundary conditions: Similarly, Laplace transformation can be used to transform these equations as follows: Fracture system: with the following boundary conditions: Zone IV is assumed to be an USVR area without any fracturing stimulation as illustrated in Figure 2(a).As we have shown in the previous process of model derivations, Zone IV simultaneously supplies gas to Zone II and Zone III.When the stress sensitivity is considered, the governing function can be presented as follows: The upper condition is assumed to be impermeable: The lower boundary is adjacent to Zone II and Zone III along with the y direction.Specifically, this lower boundary condition can be presented as follows: Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows: with the following upper and lower boundary conditions:

Geofluids
Laplace transformation is used to transform these equations to the Laplace domain: 2.3.Solution of Multilinear Model.On the basis of solution methods proposed by Ozkan et al. [2] and Stalgorova and Mattar [10], these coupled multilinear models from Equations ( 6) to (36) can be solved starting with Zone IV.The general solutions of partial differential Equations ( 35) and ( 36) can be presented as follows: The upper boundary condition can derive the relationship between A 4 and B 4 : Two lower boundary conditions as shown in Equation (31) could get different solutions: After obtaining the solutions in Zone IV, the equations for Zone II and Zone III should be coupled together to simultaneously obtain the analytical solutions.After substituting Equation (39) to Equations ( 17) and ( 27), the effects of Zone IV can be incorporated into Zone II and Zone III: Zone III : To simplify the notation, Similarly, the general solutions for Equations ( 38) and (39) can be described as follows: Zone III: There are in total 2ðN + 1Þ coefficients that need to be determined.That is to say, we should write 2ðN + 1Þ equations to obtain unique solutions.Specifically, there are N interfaces in the coupled area of Zone II and Zone III which allows us to write 2N equations.
At the interface of neighboring small subsections in Zone III, γ 3,jD A 3,j e α j l j,D + B 3,j e −α j l j,D = γ 3,j+1D A 3,j+1 e α j+1 l j,D + B 3,j+1 e −α j+1 l j,D , k 3,jD α j A 3,j e α j l j,D − α j B 3,j e −α j l j,D = k 3,j+1D α j+1 A 3,j+1 e α j+1 l j,D − α j+1 B 3,j+1 e −α j+1 l j,D : ð43Þ At the interface of Zone II and Zone III, Two other equations can be further obtained at the inner boundary of Zone II and outer boundary of Zone III: All these equations can be used to determine the required coefficients.After obtaining the solution at Zone II, we substitute this equation into Equation (13) as follows: The general solutions should be Conditioning to the inner and outer boundary conditions, these two coefficients are as follows: Finally, we can obtain the bottom-hole pressure: Following the methods used by Xu et al. [25], the zeroorder perturbation solution of the bottom-hole pressure in the Laplace space considering the wellbore storage C D and the skin factor S is obtained: Finally, by using Stehfest numerical inversion methods proposed by Stehfest [26], the zero-order perturbation solution of the bottom-hole pressure in real space is obtained, and the real bottom-hole pressure m wD can be obtained by Equation (9).

Model Verification
In this section, we will simplify our proposed model to make a comparison with other models.Wu et al. [27] have established a similar multilinear flow model of multifractured horizontal wells in stress-sensitive SGR.Further analysis of the proposed model in this work reveals that if we do not subdivide the inner SRV subsections, our model can be simplified to be a four-linear flow model as proposed by Wu et al. [27].To verify our model, comparison is made with the model proposed by Wu et al. [27].Some basic parameter settings are listed in Table 1.As shown in Figure 2, there is perfect agreement between the results of the two models.
In addition, our new model has the ability to subdivide the inner SRV zones into small subsections with different permeability.Our model can be equivalent to standard four-linear flow models if all subsections have the same permeability.In a numerical experiment, we subdivide the SRV Zone III into 1 (without division), 5, and 20 subsections.Figure 3 shows the evolution of dimensionless pressure and its derivative curves.The number of subsections has no influence on the pressure curves, which indicates that our pro-posed models can be generalized to be any multilinear flow models (MLFM).

Results and Discussions
4.1.Transient Pressure Behavior Analysis.Figure 4 illustrates the typical pressure and pressure derivative curves of a multifractured horizontal well in stress-sensitive SGR on the basis of proposed MLFM.The physical models of SGR consisting of SRV and USRV can be seen in Figure 2. We take the reference length L ref as 0.1 m.The values of the relevant dimensionless variables used in this numerical experiments are as follows: C f D = 10, k 2D = 100, k 3D = 10, k 4D = 1, ω 3D = 0:01, λ = 0:15, m DL = 0:12, V DL = 0:1, γ D = 0:05,C D = 0:001, and S = 0:1.We can approximately observe five flow stages from this log-log of dimensionless pseudopressure and its derivative curves.
Stage I: wellbore storage flow regime.This is a very common flow regime at the early stage of hydrogen extraction through wells.The wellbore storage coefficient and the skin can be obtained by fitting the plotted pressure curves to the real data.The detailed explanation of this stage can be referred to the literature [28][29][30].
Stage II: the formation bilinear flow in SRV Zone I.This flow regime is a transition between fracture linear flow and formation linear flow.The existence of highly permeable hydraulic fractures speeds up the liquid supply from formation to the hydraulic fractures at the direction of being perpendicular to the fracture faces.In general, a 1/4-slope straight line (Figure 4) will be observed on the dimensionless derivative pressure curve.We also should note that    10 Geofluids due to the high permeability of hydraulic fractures, the fracture linear flow is also blurred in this stage.This phenomenon is also consistent with the observations from the realistic applications.
Stage III: the formation linear flow in SRV Zone II and the transition flow between SRV Zone II and Zone IV.In this stage, a 1/2-slope straight line (Figure 4) can be clearly identified from the dimensionless pseudopressure derivative curve.In addition, the fluid in Zone IV also starts to supply SRV Zone II.The fluid supply of Zone IV delays the pressure depletion evidently.We can noticeably observe that a "recess" shape exists in the pressure derivative curve.
Stage IV: the transition flow between SRV Zone II and SRV Zone III.The stress sensitivity begins to take effect in this regime, and the front of the pressure wave reaches the boundary of the SRV.However, due to the lower permeability of the unstimulated Zone IV, there is not enough fluid supply from the outer region.This stage can be used to identify the boundary-dominated flow.
Stage V: the formation linear flow in outer SRV Zone III.We have idealized the outer SRV zone as dual-porosity media in SGR.The fluid will feed from the matrix to the natural fracture in this region.In this stage, the typical feature of this flow behavior is that a 1/2-slope straight line (Figure 4) occurs on the dimensionless derivative pressure curve.The effects of stress sensitivity on the transient pressure response become evident in this stage, which will be systematically analyzed in the following part.

Sensitivity Analysis.
In this section, on the basis of our proposed multilinear flow model, comprehensive sensitivity analysis of key parameters to well performance is implemented to systematically investigate the transient pressure behaviors of the multifractured horizontal well in stresssensitive SGR.We will separately investigate these two situations.

Effect of Stress Sensitivity
When we make an assumption that there is a constant permeability module for all zones, three cases were studied in which the dimensionless permeability modulus γ D is equal to 0, 0.01, and 0.02, as illustrated in Figure 5(a).It can be seen from Figure 5(a) that as the dimensionless permeability modulus increases, the dimensionless pressure and its derivative curves rise gradually; the stress sensitivity mainly affects the flow behaviors in intermediate and later flow stages.In stress-sensitive shale gas reservoirs, as the fluid is produced, the gradual reduction of formation pressure will result in a decrease of the permeability of the system and growth of pressure depletion.When the dimensionless permeability modulus increases to a certain value, the pressure derivative curve rises up significantly in later periods, showing the characteristic of a closed boundary.
When we assume that each zone has different permeability modules, three cases also were studied as illustrated in Figure 5(a).As we can see from Figure 5(b), the stress sensitivity will have significant influences on the entire flow stages, which is significantly different from the previous situation.This finding is very important for us to enhance the gas recovery from shale gas reservoir.We need to design a reasonable production scheme to maintain the formation pressure.Both Zone I and SRV Zone II are closed to the wellbore; as a result, the pressure drops are significant and hence causes severe stress sensitivity.Thus, larger stress sensitivity coefficients in the SRV zone generally induce larger pressure drops.4.2.2.Effect of the Mobility Capacity of Zone I and SRV Zone II.Transient pressure and pressure derivative curves for different mobility capacity of Zone I and SRV Zone II are illustrated in Figure 6.The values of relevant parameters are listed as follows: k 3D = 10, k 4D = 1, ω 3D = 0:01, λ = 0:15, m DL = 0:12, V DL = 0:1, γ D = 0:05,, C D = 0:001, and S = 0:1.By analysis, the mobility capacities of Zone I and SRV Zone II are represented by dimensionless fracture conductivity C f D and dimensionless formation permeability k 2D .In this numerical experiment, the SRV Zone III is not divided into subsections, e.g., N = 1.It can be seen from Figure 6(a) that dimensionless fracture conductivity C f D has significant influence on the flow behaviors of the earlier time, especially the wellbore storage and skin effect flow.Other flow stages are almost not affected.This result is consistent with our previous observations that fracture linear flow is easily blurred by the wellbore storage and skin effect flow stage due to its high mobility capacity.Specifically, as C f D decreases, the duration of wellbore storage and skin effect flow period will be shortened and the starting time of the formation bilinear flow period in the SRV Zone II will be advanced.It is mainly because the mobility capacity determines flow capacity contrast of the SRV and hydraulic fracture.Therefore, the pressure and derivative responses associated with the SRV zone will surely enlarge; the curves of the transient response will rise in the wellbore storage and skin effect flow.
The mobility capacity k 2D determines gas flow capacity in the SRV Zone II.We set the mobility capacity k 2D separately to be 100, 200, and 300 to investigate its effects on the pressure curves.It can be seen from Figure 6(b) that the dimensionless formation permeability k 2D almost influences the entire flow stages.Specifically, the flow regimes of the shape of type curves are not distorted, the pressure only rises up as the dimensionless formation permeability k 2D decreases.This means that a large value of k 2D increases the gas flow capacity in the SRV zone and therefore leads to a small pressure drop.An effective maintenance of formation pressure is very significant to yield long-term gas production from a shale gas reservoir.In addition, we also find a very interesting phenomenon that a "concave" will occur at the transition flow stage between SRV Zone II and Zone III as the k 2D increases.Our proposed multilinear flow model will show a     13 Geofluids special flow stage for the dual-porosity model.The occurrence of a "concave" also can be an indicator that the permeability of the induced SRV zone is much larger than that of the outer zone.

Effect of the Mobility Capacity of Zone III and Zone IV.
Transient pressure and pressure derivative curves for different mobility capacities of the SRV Zone III and the outer region Zone IV are illustrated in Figure 7.The values of relevant parameters are listed as follows: C f D = 10, k 2D = 100, ω 3D = 0:01, λ = 0:15, m DL = 0:12, V DL = 0:1, γ D = 0:05, C D = 0:001, and S = 0:1.In this numerical experiment, the SRV Zone III is not divided into subsections as wells, e.g., N = 1.We can observe from Figure 7(a) that the mobility capacity of Zone III almost has no influence on the flow behaviors of the early and intermediate flow stages; on the contrary, the transition flow and formation linear flow regimes are severely dominated.These results are also in agreement with our common understanding.Specifically, as the mobility capacity of Zone III k 3D decreases, the duration of transition flow period will be enlarged and the starting time of the formation linear flow period and boundary-dominated flow stage in the SRV Zone III will be delayed.
The mobility capacity k 4D determines gas flow capacity in the outer Zone IV.We set mobility capacity k 4D separately to be 10, 20, and 30 to investigate its effects on the pressure curves.Figure 7(b) apparently demonstrates that the dimensionless formation permeability k 4D has great influence on the formation bilinear flow stage in SRV Zone II.This is very different from the effects of mobility capacity k 3D in SRV Zone III.In the process of model development, we assume that the outer Zone IV simultaneously supplies gas to SRV Zone II and SRV Zone III.Because SRV Zone II is adjacent to the wellbore and thus yields larger pressure drop, the outer Zone IV is more prone to supply gas to SRV Zone II, as illustrated in Figure 7(b).As the mobility capacity k 4D decreases, the gas supply will be reduced as well; as a result, the dimensionless pressure curves will move upward.8, the width and length of the outer region have mainly affected different flow regimes in the SRV Zone II and SRV Zone III.Specifically speaking, the width of the outer region mainly has influence on the formation linear flow regime in SRV Zone III, while the width of the outer region mainly impacts the formation bilinear and transition flow regimes.As the size of the outer region decreases, the boundary conditions will interfere with the flow regimes in the SRV zone in advance; as a result, the dimensionless pressure and its derivatives will be higher.This phenomenon is very similar to the effects of the mobility capacity of the SRV Zone III and Zone IV.Enlarging the size of the outer region can decrease the negative effects of the boundary condition, and increasing the mobility capacity of Zone III and Zone IV can decrease the flow resistance in SRV zone, which is beneficial to obtain a high production rate.and the dimensionless size of SRV Zone II L D is 300, 500, and 700 (Figure 9).It can be seen from Figure 9 that the dimensionless size of SRV Zone II can affect all the flow stages after the formation bilinear flow regimes.The larger size of SRV Zone II has (1) a later end of formation bilinear in the SRV Zone II, (2) a postponed beginning of formation linear flow in SRV Zone III, and (3) the lower values of dimensionless pressure and its derivatives in the transition flow regime between SRV Zone II and SRV Zone III.All these phenomena might provide useful information to identify the sizes of SRV after massive hydraulic fracturing.Figure 9 also shows that the larger the SRV size, the smaller the dimensionless pressure and its derivatives.Smaller dimensionless pressure indicates lower pressure depletion in the formation, which is beneficial to obtain a high production rate.In this paper, Zone II is a highly permeable SVR zone which is adjacent to the main hydraulic fracture as illustrated in Figure 2(a).Zone II is the gas flow path from SGR to the wellbore.Increasing the size of this zone will delay the negative effects of the outer boundary, which is beneficial to obtain a high production rate as well.4.2.6.Effect of Coefficient of Permeability Variation.In this paper, to characterize spatially varying permeability due to the nonuniform distribution of proppant within SRV Zone III, we further divide Zone III into N small zones as illustrated in Figure 1 To systematically represent the spatial variation of permeability in Zone III, two mathematical formulas, including linear and logarithmic functions, are used to describe the reduction of permeability away from the main hydraulic fracture.Figures 10(a) and 10(b) show the dimensionless permeability for each subsection.Figures 10(c) and 10(d) depict the dimensionless permeability as a function of distance to the main hydraulic fracture.In contrast to the linear function, the logarithmic function yields a more rapid reduction of dimensionless permeability.What is more, increasing the number of subsections can generate a continuous reduction of permeability.
Figures 10(e) and 10(f) show the effect of coefficient of permeability variation in SRV Zone III on the transient behavior.
Because the linear function generates a continuous reduction 14 Geofluids of permeability, large parts of SRV Zone III still preserve relatively high permeability, and only a small part of regions which is adjacent to the outer boundary will get small permeability.
Under this condition, the formation linear flow regime will be mainly impacted as shown in Figure 10(e); as the number of subsections increases, the reduction rate of permeability  15 Geofluids will decrease; as a result, the negative of the boundary condition will be delayed.On the contrary, the logarithmic function easily leads to a rapid reduction of permeability as illustrated in Figure 10(d).Almost the entire flow regimes are severely influenced.This phenomenon can be considered to be a reliable indicator to judge whether the variation of  Geofluids permeability follows a logarithmic function or not.In the real-field applications, a linear variation of permeability induced by a massive hydraulic fracturing operation has the potential to obtain a high production rate.

Conclusion
An analytical multilinear flow model is proposed for multifractured horizontal wells with SRV in shale gas reservoir.
The transient pressure curves are also established to analyze the performance of this new model.Some key results could be summarized as follows: (1) On the basis of perturbation technique and Laplace transformation, an easy-to-simulate model with the consideration of stress sensitivity and SRV is obtained in Laplace space, and then the solution of transient pressure behaviors for a multifractured horizontal well with SRV is finally obtained by Stehfest numerical inversion algorithm (4) SRV Zone II provides the main flow path to supply gas from the shale gas reservoir to the wellbore.Increasing the size of the outer region could delay the negative influences of boundary conditions.The permeability of SRV Zone III mainly impacts the formation linear flow regimes, while the permeability of Zone IV has significant influence on the formation linear flow regime in SRV Zone II (5) To systematically represent the spatial variation of permeability in Zone III, both linear and logarithmic functions are used to describe the reduction of permeability.In contrast to the linear function, the logarithmic function yields a more rapid reduction of dimensionless permeability.As the number of subsections increases, the reduction rate of permeability will be decreased; as a result, the negative of the boundary condition will be delayed.The logarithmic function easily leads to a rapid reduction of permeability, and hence, almost the entire flow regimes are severely influenced

Appendix
Here, we take Zone I as an example to derive the governing equations as shown in the main content.Based on the mass conservation principle, Darcy law, the shale gas flow along the fracture direction can be presented as follows: ðA:4Þ The outer boundary condition: Peripheral zone of SRV (b) Illustration of a multilinear flow model in a shale gas reservoir

Figure 1 :
Figure 1: Physical model of multifractured horizontal well in a shale gas reservoir.
Wu et al, 2019) dm D (Wu et al, 2019) m D (new model) dm D (new model)

Figure 2 :
Figure 2: Comparison between the model in this paper and that in Wu et al. [27].

Figure 3 :
Figure 3: Investigation of number of subdivisions of SRV Zone III on the type curves.

Figure 4 :
Figure 4: Typical curves of multifractured horizontal well in stress-sensitive SGR formation.

Figure 5 :
Figure 5: The effect of permeability modulus on pressure and derivative curves: (a) constant permeability module for all zones; (b) each zone has different permeability module.

Figure 6 :
Figure 6: The effect of mobility capacity of Zone I and the SRV Zone II on pressure and derivative curves.

4. 2 . 5 .
Effect of the Size of SRV Zone II.

Figure 7 :
Figure 7: The effect of mobility capacity of Zone III and Zone IV on pressure and derivative curves.

Figure 8 :
Figure 8: The effect of size of the outer region on pressure and derivative curves.

( 2 )
This newly proposed multilinear flow model (MLFM) is validated by comparing with an existing model, and a perfect agreement has been obtained (3) Approximately five flow stages can be identified: wellbore storage and skin effect flow, the formation bilinear flow in the inner SRV, the formation linear flow in the SRV zone and the transition flow between SRV Zone II and Zone IV, the transition flow between SRV Zone II and Zone III, and the formation linear flow in Zone III

Figure 9 :
Figure 9: The effect of the size of the SRV Zone II on pressure and derivative curves.

Table 1 :
Basic data used for the model.
Introducing the dimensionless definition, the term −ð2πk ref hx f /q gsc μ gi B gi Þ is multiplied to the two sides: