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

Authors: Guangdong Wang; Ailin Jia; Yunsheng Wei; Cong Xiao
Journal: Geofluids
 (2020)

---

## Abstract

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 first 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 flow model (MLFM) is obtained and validated as a comparison with the previous models. The flow regimes are identified, 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 significance for an effective recovery of SGR resources.


---

## Body

## 1. 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 high-conductivity 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 reservoir 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 complexity 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.


## 2. Mathematical Model and Analytical Solution

### 2.1. Model Descriptions

As illustrated in Figure1(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 xf and the width and length of the entire zone are xe and ye, 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 k1, k2, k3, and k4.Figure 1
Physical model of multifractured horizontal well in a shale gas reservoir.
(a)
The illustration of microseismic response
(b)
Illustration of a multilinear flow model in a shale gas reservoir
(c)
Illustration of a further subdivision of Zone IIAs 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 intoN 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. Single-phase and microcompressible gas is assumed to flow in SGR, which follows Darcy’s Law.


### 2.2. Mathematical Formula of Multilinear Model

To begin with, we define the pseudopressure as follows:(1)mp=μgZgipi∫prefppμgZgdp,and then,
(2)mnD=2πkrefhmpi−mpnqgscμgiBgi.According to the theory proposed by Pedrosa Jr. [22], stress-sensitive permeability could be described as follows:
(3)kim=k0e−γim0−mi,i=1,2,3,4,where the reference permeability k0 is the initial permeability in SGR, m2, 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
(4)VE=VLmmL+m,where VE is the adsorption equilibrium concentration, sm3/m3; VL is the Langmuir adsorption concentration, sm3/m3; PL is the Langmuir pressure, MPa; and mL is the Langmuir pressure, MPa2/(mPa·s).For the convenience and simplicity of deducing formulas, some dimensionless parameters are introduced first:(5)mi=krefhm0−mi1.842qgscμgiBgi,i=1,2,3,4,tD=3.6×10−3kreftμgiφCgrefxf2,xD=xxf,yD=yxf,lj,D=ljxf,wfD=wfxf,γiD=1.842qscμgBgkrefhγi,CfD=k1wfkrefxf,kiD=kikref,ωiD=φiCgiφCgref,ω3mD=φ3mCgiφCgref,λ=σk3mkrefLref2,mDL=krefhmL−mi1.842qgscμgiBgi,VDL=2pscTμiTSCkrefhmL−miφ3mVL1.842qgscμgiBgi.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.


#### 2.2.1. 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:(6)∂2m1∂x2+γ1∂m1∂x2+2wfk2,1k1∂m2,1∂yy=wf/2=eλ1m0−m1φ1μCgk1∂m1∂t,with the inner and outer boundary conditions:
(7)e−γ1m0−m1∂m1∂xx=0=qgscμgiBgi2wfhk1,∂m1∂xx=xf=0.Substituting predefined dimensionless variables to Equations (6) and (7), their dimensionless formula is as follows:
(8)∂2m1D∂xD2+γ1D∂m1D∂xD2+2k2,1DC1D∂m2,1D∂yDyD=wfD/2=eγ1Dm1Dω1Dk1D∂m1D∂tD,with e−γ1Dm1D∂m1D/∂xDxD=0=−π/C1D and ∂m1D/∂xDxD=1=0.Equation (8) is a strongly nonlinear partial differential equation, which is not convenient to be solved analytically. 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:
(9)miD=−1γiDln1−γiDηiD,i=1,2,3,4.According to the simplified method proposed by [6], due to the dimensionless permeability modulus γD which is usually a small value, ηjDcan be expanded as a power series in the parameter γD.
(10)ηjD=ηjD0+γDηjD1+γD2ηjD2+γD3ηjD3+⋯,11−γDηjD=1+γDηjD+γD2ηjD2+γD3ηjD3+⋯.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.
(11)∂2η1D∂xD2+2k2,1DC1D∂η2,1D∂yDyD=wfD/2=ω1Dk1D∂η1D∂tD,with ∂η1D/∂xDxD=0=−π/C1D and ∂η1D/∂xDxD=1=0.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.(12)d2η1D¯dxD2+2k2,1DC1Ddη2,1D¯dyDyD=wfD/2=sω1Dk1Dη1D¯,with dη1D¯/dxDxD=0=−1/sπ/C1D and dη1D¯/dxDxD=1=0.


#### 2.2.2. Zone II

Zone II is a highly permeable SVR zone which is adjacent to the main hydraulic fracture as illustrated in Figure1(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:
(13)∂2m2∂y2+γ2∂mm2∂y2+k4k2∂m4∂yx=xf/2=eγ2m0−m2φ2μCgk2∂m2∂t,wf2<x<l1.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:(14)m2y=wf/2=m1y=wf/2,m2y=l1=m3y=l1.Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows:(15)∂2η2D∂yD2+k4Dk2D∂η4D∂xDxD=1=ω2Dk2D∂η2D∂tD,wfD2<yD<l1,Dwith the following boundary conditions:
(16)γ1Dη1Dy=wfD/2=γ2Dη2y=wfD/2,γ2Dη2Dy=l1,D=γ3,1Dη3,1Dy=l1,D.Similarly, Laplace transformation can be used to transform these equations from the time domain to the Laplace domain.(17)d2η2D¯dyD2+k4Dk2Ddη4D¯dxDxD=1=sω2Dk2Dη2D¯,wfD2<yD<l1,D,with the following boundary conditions:
(18)γ1Dη1D¯y=wfD/2=γ2Dη2D¯y=wfD/2,γ2Dη2D¯y=l1,D=γ3,1Dη3,1D¯y=l1,D.


#### 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 intoN 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 considered, the governing function can be presented as follows:
(19)Fracturesystem:∂2m3,j∂y2+γ3,j∂mm3,j∂y2+k4k3,j∂m4∂yx=xf+σk3mk3,jm3m−m3,j=eγ3,jm0−m3,jφ3,jμCgk3,j∂m3,j∂t,lj<x<lj+1,(20)Matrixsystem:−σk3mm3m−m3,j=φ3mμCg∂m3m∂t+pscTTSC∂V∂t.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.(21)m3,jy=lj=m3,j+1y=lj,k3,je−γ3,jm0−m3,j∂m3,j∂yy=lj=k3,j+1e−γ3,j+1m0−m3,j+1∂m3,j+1∂yy=lj,j=2,⋯,N−1.In particular, the first small zone will directly connect to Zone II, and the right side of theN-th small zone is considered to be impermeable. The boundary conditions could be explicitly presented as follows:
(22)m2y=l1=m3,1y=l1,k2e−γ2m0−m2∂m2∂yy=l1=k3,1e−γ3,1m0−m3,1∂m3,1∂yy=l1,∂m3,N∂yy=ye=0.Obtained from Fick’s first law of diffusion, the diffusion rate of shale gas can be expressed as(23)∂V∂t=−ϕ3m∂VE∂t.Substituting the Langmuir isotherm Equation (4) into Equation (23), we can obtain
(24)∂V∂t=−φ3mVLmLmL+m3m2∂m3m∂t.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:
(25)Fracturesystem:∂2η3,jD∂yD2+k4Dk3,jD∂m4D∂yDxD=1+λk3,jDη3mD−η3,jD=ω3,jDk3,jD∂η3,jD∂tD,lj,D<yD<lj+1,D,Matrixsystem:−λη3mD−η3,jD=ω3mD∂η3mD∂tD+VDLmDLmDL+mDi−m3mD2∂η3mD∂tD,with the following boundary conditions:
(26)γ3,jDη3,jDxD=lj,D=γ3,j+1Dη3,j+1DxD=lj,D,∂η3,jD∂yDyD=lj,D=k3,j+1Dk3,jD∂η3,j+1D∂yDyD=lj,D,j=2,⋯,N−1,∂η3,1D∂yDyD=l1,D=k2Dk3,1D∂η2D∂yDyD=l2,D,∂η3,ND∂yDyD=lN,D=0.Similarly, Laplace transformation can be used to transform these equations as follows:(27)Fracturesystem:d2η3,jD¯dyD2+k4Dk3,jDdm4D¯dyDxD=1+λk3,jDη3mD¯−η3,jD¯=ω3,jDk3,jDsη3,jD¯,lj,D<yD<lj+1,D,(28)Matrixsystem:−λη3mD¯−η3,jD¯=ω3mDsη3mD¯,with the following boundary conditions:
(29)γ3,jDη3,jD¯xD=lj,D=γ3,j+1Dη3,j+1D¯xD=lj,D,dη3,jD¯dyDyD=lj,D=k3,j+1Dk3,jDdη3,j+1D¯dyDyD=lj,D,j=2,⋯,N−1,γ3,1Dη3,1D¯=γ2Dη2D¯,dη3,1D¯dyDyD=l1,D=k2Dk3,1Ddη2D¯dyDyD=l2,D,dη3,ND¯dyDyD=lN,D=0.


#### 2.2.4. Zone IV

Zone IV is assumed to be an USVR area without any fracturing stimulation as illustrated in Figure2(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:
(30)∂2m4∂x2+γ4∂mm4∂x2=eγ4m0−m4φ4μCgk4∂m4∂t,xf<x<xe.Figure 2
Comparison between the model in this paper and that in Wu et al. [27].The upper condition is assumed to be impermeable:(31)∂m4∂xx=xe=0.The lower boundary is adjacent to Zone II and Zone III along with they direction. Specifically, this lower boundary condition can be presented as follows:
(32)m4x=xf=m2x=xf,wf2<y<l1,m3,jx=xf,lj−1<y<lj.Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows:(33)∂2η4D∂xD2=ω4dk4D∂η4D∂tD,1<xD<xeD,with the following upper and lower boundary conditions:
(34)∂η4D∂xDxD=xeD=0,γ4Dη4DxD=1=γ2Dη2DxD=1,wfD2<yD<l1,D,γ3,jDη3,jDxD=1,lj−1,D<yD<lj,D.Laplace transformation is used to transform these equations to the Laplace domain:(35)d2η4D¯dxD2=sω4dk4Dη4D¯,1<xD<xeD,(36)dη4D¯dxDxD=xeD=0,γ4Dη4D¯xD=1=γ2Dη2D¯xD=1,wfD/2<yD<l1,D,γ3,jDη3,jD¯xD=1,lj−1,D<yD<lj,D.


### 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:
(37)η4D¯=A4esω4d/k4DxD+B4e−sω4d/k4DxD.The upper boundary condition can derive the relationship betweenA4 and B4:
(38)B4=A4e2sω4d/k4DxeD.Two lower boundary conditions as shown in Equation (31) could get different solutions:
(39)η4D¯=γ2Dcoshsω4d/k4DxD−xeDγ4Dcoshsω4d/k4D1−xeDη2D¯,wfD2<yD<l1,D,γ3jDcoshsω4d/k4DxD−xeDγ4Dcoshsω4d/k4D1−xeDη3jD¯,lj−1,D<yD<lj,D.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:
(40)ZoneII:d2η2D¯dyD2=γ2Dk4Dγ4Dk2Dsω4Dk4Dtanhsω4dk4D1−xeD+sω2Dk2Dη2D¯,wfD2<yD<l1,D,ZoneIII:d2η3,jD¯dyD2=γ3jDk4Dγ4Dk3,jDsω4Dk4Dtanhsω4dk4D1−xeD+λω3mDsλ+ω3mDsk3,jD+sω3,jDk3,jDη3,jD¯,lj,D<yD<lj+1,D.To simplify the notation,(41)β=γ2Dk4Dγ2Dk3Dsω4Dk4Dtanhsω4dk4D1−xeD+sω2Dk2D,αj=γ3jDk4Dγ4Dk2,jDsω4Dk4Dtanhsω4dk4D1−xeD+λω3mDsλ+ω3mDsk3,jD+sω3,jDk3,jD,j=1,2,⋯,N.Similarly, the general solutions for Equations (38) and (39) can be described as follows:
(42)ZoneII:η2D¯=A2eβxD+B2e−βxD,wfD2<yD<l1,D,ZoneIII:η3,jD¯=A3,jeαjxD+B3,je−αjxD,lj,D<yD<lj+1,D.There are in total2N+1 coefficients that need to be determined. That is to say, we should write 2N+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,(43)γ3,jDA3,jeαjlj,D+B3,je−αjlj,D=γ3,j+1DA3,j+1eαj+1lj,D+B3,j+1e−αj+1lj,D,k3,jDαjA3,jeαjlj,D−αjB3,je−αjlj,D=k3,j+1Dαj+1A3,j+1eαj+1lj,D−αj+1B3,j+1e−αj+1lj,D.At the interface of Zone II and Zone III,(44)γ3,1DA3,1eα1l1,D+B3,1e−α1l1,D=γ2DA2eβl1,D+B2e−βl1,D,k3,1Dα1A3,1eα1l1,D−α1B3,1e−α1l1,D=k2DβA2eβl1,D−βB2e−βl1,D.Two other equations can be further obtained at the inner boundary of Zone II and outer boundary of Zone III:(45)γ1Dη1D¯=γ2DA2eα1wfD/2+B2e−α1wfD/2,A3,NeβyeD−B3,Ne−βyeD=0.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:
(46)d2η1D¯dxD2=sω1Dk1D−2k2DC1DβA2eβwfD/2−B2e−βwfD/2η1D¯.The general solutions should be(47)η1D¯=A1eχxD+B1e−χxD,0<xD<1where χ=sω1D/k1D−2k2D/C1DβA2eβwfD/2−B2e−βwfD/2.Conditioning to the inner and outer boundary conditions, these two coefficients are as follows:(48)A1=πsC1Dχe2χ−1,B1=πe2χsC1Dχe2χ−1.Finally, we can obtain the bottom-hole pressure:(49)ηwD¯=η1D¯xD=0=πe2χ+1sC1Dχe2χ−1.Following the methods used by Xu et al. [25], the zero-order perturbation solution of the bottom-hole pressure in the Laplace space considering the wellbore storage CD and the skin factor S is obtained:
(50)η¯wDS,CD=sη¯wD+Ss+CDs2sη¯wD+S.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 mwD can be obtained by Equation (9).


## 2.1. Model Descriptions

As illustrated in Figure1(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 xf and the width and length of the entire zone are xe and ye, 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 k1, k2, k3, and k4.Figure 1
Physical model of multifractured horizontal well in a shale gas reservoir.
(a)
The illustration of microseismic response
(b)
Illustration of a multilinear flow model in a shale gas reservoir
(c)
Illustration of a further subdivision of Zone IIAs 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 intoN 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. Single-phase and microcompressible gas is assumed to flow in SGR, which follows Darcy’s Law.


## 2.2. Mathematical Formula of Multilinear Model

To begin with, we define the pseudopressure as follows:(1)mp=μgZgipi∫prefppμgZgdp,and then,
(2)mnD=2πkrefhmpi−mpnqgscμgiBgi.According to the theory proposed by Pedrosa Jr. [22], stress-sensitive permeability could be described as follows:
(3)kim=k0e−γim0−mi,i=1,2,3,4,where the reference permeability k0 is the initial permeability in SGR, m2, 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
(4)VE=VLmmL+m,where VE is the adsorption equilibrium concentration, sm3/m3; VL is the Langmuir adsorption concentration, sm3/m3; PL is the Langmuir pressure, MPa; and mL is the Langmuir pressure, MPa2/(mPa·s).For the convenience and simplicity of deducing formulas, some dimensionless parameters are introduced first:(5)mi=krefhm0−mi1.842qgscμgiBgi,i=1,2,3,4,tD=3.6×10−3kreftμgiφCgrefxf2,xD=xxf,yD=yxf,lj,D=ljxf,wfD=wfxf,γiD=1.842qscμgBgkrefhγi,CfD=k1wfkrefxf,kiD=kikref,ωiD=φiCgiφCgref,ω3mD=φ3mCgiφCgref,λ=σk3mkrefLref2,mDL=krefhmL−mi1.842qgscμgiBgi,VDL=2pscTμiTSCkrefhmL−miφ3mVL1.842qgscμgiBgi.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.


### 2.2.1. 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:(6)∂2m1∂x2+γ1∂m1∂x2+2wfk2,1k1∂m2,1∂yy=wf/2=eλ1m0−m1φ1μCgk1∂m1∂t,with the inner and outer boundary conditions:
(7)e−γ1m0−m1∂m1∂xx=0=qgscμgiBgi2wfhk1,∂m1∂xx=xf=0.Substituting predefined dimensionless variables to Equations (6) and (7), their dimensionless formula is as follows:
(8)∂2m1D∂xD2+γ1D∂m1D∂xD2+2k2,1DC1D∂m2,1D∂yDyD=wfD/2=eγ1Dm1Dω1Dk1D∂m1D∂tD,with e−γ1Dm1D∂m1D/∂xDxD=0=−π/C1D and ∂m1D/∂xDxD=1=0.Equation (8) is a strongly nonlinear partial differential equation, which is not convenient to be solved analytically. 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:
(9)miD=−1γiDln1−γiDηiD,i=1,2,3,4.According to the simplified method proposed by [6], due to the dimensionless permeability modulus γD which is usually a small value, ηjDcan be expanded as a power series in the parameter γD.
(10)ηjD=ηjD0+γDηjD1+γD2ηjD2+γD3ηjD3+⋯,11−γDηjD=1+γDηjD+γD2ηjD2+γD3ηjD3+⋯.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.
(11)∂2η1D∂xD2+2k2,1DC1D∂η2,1D∂yDyD=wfD/2=ω1Dk1D∂η1D∂tD,with ∂η1D/∂xDxD=0=−π/C1D and ∂η1D/∂xDxD=1=0.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.(12)d2η1D¯dxD2+2k2,1DC1Ddη2,1D¯dyDyD=wfD/2=sω1Dk1Dη1D¯,with dη1D¯/dxDxD=0=−1/sπ/C1D and dη1D¯/dxDxD=1=0.


### 2.2.2. Zone II

Zone II is a highly permeable SVR zone which is adjacent to the main hydraulic fracture as illustrated in Figure1(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:
(13)∂2m2∂y2+γ2∂mm2∂y2+k4k2∂m4∂yx=xf/2=eγ2m0−m2φ2μCgk2∂m2∂t,wf2<x<l1.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:(14)m2y=wf/2=m1y=wf/2,m2y=l1=m3y=l1.Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows:(15)∂2η2D∂yD2+k4Dk2D∂η4D∂xDxD=1=ω2Dk2D∂η2D∂tD,wfD2<yD<l1,Dwith the following boundary conditions:
(16)γ1Dη1Dy=wfD/2=γ2Dη2y=wfD/2,γ2Dη2Dy=l1,D=γ3,1Dη3,1Dy=l1,D.Similarly, Laplace transformation can be used to transform these equations from the time domain to the Laplace domain.(17)d2η2D¯dyD2+k4Dk2Ddη4D¯dxDxD=1=sω2Dk2Dη2D¯,wfD2<yD<l1,D,with the following boundary conditions:
(18)γ1Dη1D¯y=wfD/2=γ2Dη2D¯y=wfD/2,γ2Dη2D¯y=l1,D=γ3,1Dη3,1D¯y=l1,D.


### 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 intoN 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 considered, the governing function can be presented as follows:
(19)Fracturesystem:∂2m3,j∂y2+γ3,j∂mm3,j∂y2+k4k3,j∂m4∂yx=xf+σk3mk3,jm3m−m3,j=eγ3,jm0−m3,jφ3,jμCgk3,j∂m3,j∂t,lj<x<lj+1,(20)Matrixsystem:−σk3mm3m−m3,j=φ3mμCg∂m3m∂t+pscTTSC∂V∂t.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.(21)m3,jy=lj=m3,j+1y=lj,k3,je−γ3,jm0−m3,j∂m3,j∂yy=lj=k3,j+1e−γ3,j+1m0−m3,j+1∂m3,j+1∂yy=lj,j=2,⋯,N−1.In particular, the first small zone will directly connect to Zone II, and the right side of theN-th small zone is considered to be impermeable. The boundary conditions could be explicitly presented as follows:
(22)m2y=l1=m3,1y=l1,k2e−γ2m0−m2∂m2∂yy=l1=k3,1e−γ3,1m0−m3,1∂m3,1∂yy=l1,∂m3,N∂yy=ye=0.Obtained from Fick’s first law of diffusion, the diffusion rate of shale gas can be expressed as(23)∂V∂t=−ϕ3m∂VE∂t.Substituting the Langmuir isotherm Equation (4) into Equation (23), we can obtain
(24)∂V∂t=−φ3mVLmLmL+m3m2∂m3m∂t.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:
(25)Fracturesystem:∂2η3,jD∂yD2+k4Dk3,jD∂m4D∂yDxD=1+λk3,jDη3mD−η3,jD=ω3,jDk3,jD∂η3,jD∂tD,lj,D<yD<lj+1,D,Matrixsystem:−λη3mD−η3,jD=ω3mD∂η3mD∂tD+VDLmDLmDL+mDi−m3mD2∂η3mD∂tD,with the following boundary conditions:
(26)γ3,jDη3,jDxD=lj,D=γ3,j+1Dη3,j+1DxD=lj,D,∂η3,jD∂yDyD=lj,D=k3,j+1Dk3,jD∂η3,j+1D∂yDyD=lj,D,j=2,⋯,N−1,∂η3,1D∂yDyD=l1,D=k2Dk3,1D∂η2D∂yDyD=l2,D,∂η3,ND∂yDyD=lN,D=0.Similarly, Laplace transformation can be used to transform these equations as follows:(27)Fracturesystem:d2η3,jD¯dyD2+k4Dk3,jDdm4D¯dyDxD=1+λk3,jDη3mD¯−η3,jD¯=ω3,jDk3,jDsη3,jD¯,lj,D<yD<lj+1,D,(28)Matrixsystem:−λη3mD¯−η3,jD¯=ω3mDsη3mD¯,with the following boundary conditions:
(29)γ3,jDη3,jD¯xD=lj,D=γ3,j+1Dη3,j+1D¯xD=lj,D,dη3,jD¯dyDyD=lj,D=k3,j+1Dk3,jDdη3,j+1D¯dyDyD=lj,D,j=2,⋯,N−1,γ3,1Dη3,1D¯=γ2Dη2D¯,dη3,1D¯dyDyD=l1,D=k2Dk3,1Ddη2D¯dyDyD=l2,D,dη3,ND¯dyDyD=lN,D=0.


### 2.2.4. Zone IV

Zone IV is assumed to be an USVR area without any fracturing stimulation as illustrated in Figure2(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:
(30)∂2m4∂x2+γ4∂mm4∂x2=eγ4m0−m4φ4μCgk4∂m4∂t,xf<x<xe.Figure 2
Comparison between the model in this paper and that in Wu et al. [27].The upper condition is assumed to be impermeable:(31)∂m4∂xx=xe=0.The lower boundary is adjacent to Zone II and Zone III along with they direction. Specifically, this lower boundary condition can be presented as follows:
(32)m4x=xf=m2x=xf,wf2<y<l1,m3,jx=xf,lj−1<y<lj.Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows:(33)∂2η4D∂xD2=ω4dk4D∂η4D∂tD,1<xD<xeD,with the following upper and lower boundary conditions:
(34)∂η4D∂xDxD=xeD=0,γ4Dη4DxD=1=γ2Dη2DxD=1,wfD2<yD<l1,D,γ3,jDη3,jDxD=1,lj−1,D<yD<lj,D.Laplace transformation is used to transform these equations to the Laplace domain:(35)d2η4D¯dxD2=sω4dk4Dη4D¯,1<xD<xeD,(36)dη4D¯dxDxD=xeD=0,γ4Dη4D¯xD=1=γ2Dη2D¯xD=1,wfD/2<yD<l1,D,γ3,jDη3,jD¯xD=1,lj−1,D<yD<lj,D.


## 2.2.1. 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:(6)∂2m1∂x2+γ1∂m1∂x2+2wfk2,1k1∂m2,1∂yy=wf/2=eλ1m0−m1φ1μCgk1∂m1∂t,with the inner and outer boundary conditions:
(7)e−γ1m0−m1∂m1∂xx=0=qgscμgiBgi2wfhk1,∂m1∂xx=xf=0.Substituting predefined dimensionless variables to Equations (6) and (7), their dimensionless formula is as follows:
(8)∂2m1D∂xD2+γ1D∂m1D∂xD2+2k2,1DC1D∂m2,1D∂yDyD=wfD/2=eγ1Dm1Dω1Dk1D∂m1D∂tD,with e−γ1Dm1D∂m1D/∂xDxD=0=−π/C1D and ∂m1D/∂xDxD=1=0.Equation (8) is a strongly nonlinear partial differential equation, which is not convenient to be solved analytically. 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:
(9)miD=−1γiDln1−γiDηiD,i=1,2,3,4.According to the simplified method proposed by [6], due to the dimensionless permeability modulus γD which is usually a small value, ηjDcan be expanded as a power series in the parameter γD.
(10)ηjD=ηjD0+γDηjD1+γD2ηjD2+γD3ηjD3+⋯,11−γDηjD=1+γDηjD+γD2ηjD2+γD3ηjD3+⋯.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.
(11)∂2η1D∂xD2+2k2,1DC1D∂η2,1D∂yDyD=wfD/2=ω1Dk1D∂η1D∂tD,with ∂η1D/∂xDxD=0=−π/C1D and ∂η1D/∂xDxD=1=0.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.(12)d2η1D¯dxD2+2k2,1DC1Ddη2,1D¯dyDyD=wfD/2=sω1Dk1Dη1D¯,with dη1D¯/dxDxD=0=−1/sπ/C1D and dη1D¯/dxDxD=1=0.


## 2.2.2. Zone II

Zone II is a highly permeable SVR zone which is adjacent to the main hydraulic fracture as illustrated in Figure1(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:
(13)∂2m2∂y2+γ2∂mm2∂y2+k4k2∂m4∂yx=xf/2=eγ2m0−m2φ2μCgk2∂m2∂t,wf2<x<l1.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:(14)m2y=wf/2=m1y=wf/2,m2y=l1=m3y=l1.Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows:(15)∂2η2D∂yD2+k4Dk2D∂η4D∂xDxD=1=ω2Dk2D∂η2D∂tD,wfD2<yD<l1,Dwith the following boundary conditions:
(16)γ1Dη1Dy=wfD/2=γ2Dη2y=wfD/2,γ2Dη2Dy=l1,D=γ3,1Dη3,1Dy=l1,D.Similarly, Laplace transformation can be used to transform these equations from the time domain to the Laplace domain.(17)d2η2D¯dyD2+k4Dk2Ddη4D¯dxDxD=1=sω2Dk2Dη2D¯,wfD2<yD<l1,D,with the following boundary conditions:
(18)γ1Dη1D¯y=wfD/2=γ2Dη2D¯y=wfD/2,γ2Dη2D¯y=l1,D=γ3,1Dη3,1D¯y=l1,D.


## 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 intoN 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 considered, the governing function can be presented as follows:
(19)Fracturesystem:∂2m3,j∂y2+γ3,j∂mm3,j∂y2+k4k3,j∂m4∂yx=xf+σk3mk3,jm3m−m3,j=eγ3,jm0−m3,jφ3,jμCgk3,j∂m3,j∂t,lj<x<lj+1,(20)Matrixsystem:−σk3mm3m−m3,j=φ3mμCg∂m3m∂t+pscTTSC∂V∂t.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.(21)m3,jy=lj=m3,j+1y=lj,k3,je−γ3,jm0−m3,j∂m3,j∂yy=lj=k3,j+1e−γ3,j+1m0−m3,j+1∂m3,j+1∂yy=lj,j=2,⋯,N−1.In particular, the first small zone will directly connect to Zone II, and the right side of theN-th small zone is considered to be impermeable. The boundary conditions could be explicitly presented as follows:
(22)m2y=l1=m3,1y=l1,k2e−γ2m0−m2∂m2∂yy=l1=k3,1e−γ3,1m0−m3,1∂m3,1∂yy=l1,∂m3,N∂yy=ye=0.Obtained from Fick’s first law of diffusion, the diffusion rate of shale gas can be expressed as(23)∂V∂t=−ϕ3m∂VE∂t.Substituting the Langmuir isotherm Equation (4) into Equation (23), we can obtain
(24)∂V∂t=−φ3mVLmLmL+m3m2∂m3m∂t.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:
(25)Fracturesystem:∂2η3,jD∂yD2+k4Dk3,jD∂m4D∂yDxD=1+λk3,jDη3mD−η3,jD=ω3,jDk3,jD∂η3,jD∂tD,lj,D<yD<lj+1,D,Matrixsystem:−λη3mD−η3,jD=ω3mD∂η3mD∂tD+VDLmDLmDL+mDi−m3mD2∂η3mD∂tD,with the following boundary conditions:
(26)γ3,jDη3,jDxD=lj,D=γ3,j+1Dη3,j+1DxD=lj,D,∂η3,jD∂yDyD=lj,D=k3,j+1Dk3,jD∂η3,j+1D∂yDyD=lj,D,j=2,⋯,N−1,∂η3,1D∂yDyD=l1,D=k2Dk3,1D∂η2D∂yDyD=l2,D,∂η3,ND∂yDyD=lN,D=0.Similarly, Laplace transformation can be used to transform these equations as follows:(27)Fracturesystem:d2η3,jD¯dyD2+k4Dk3,jDdm4D¯dyDxD=1+λk3,jDη3mD¯−η3,jD¯=ω3,jDk3,jDsη3,jD¯,lj,D<yD<lj+1,D,(28)Matrixsystem:−λη3mD¯−η3,jD¯=ω3mDsη3mD¯,with the following boundary conditions:
(29)γ3,jDη3,jD¯xD=lj,D=γ3,j+1Dη3,j+1D¯xD=lj,D,dη3,jD¯dyDyD=lj,D=k3,j+1Dk3,jDdη3,j+1D¯dyDyD=lj,D,j=2,⋯,N−1,γ3,1Dη3,1D¯=γ2Dη2D¯,dη3,1D¯dyDyD=l1,D=k2Dk3,1Ddη2D¯dyDyD=l2,D,dη3,ND¯dyDyD=lN,D=0.


## 2.2.4. Zone IV

Zone IV is assumed to be an USVR area without any fracturing stimulation as illustrated in Figure2(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:
(30)∂2m4∂x2+γ4∂mm4∂x2=eγ4m0−m4φ4μCgk4∂m4∂t,xf<x<xe.Figure 2
Comparison between the model in this paper and that in Wu et al. [27].The upper condition is assumed to be impermeable:(31)∂m4∂xx=xe=0.The lower boundary is adjacent to Zone II and Zone III along with they direction. Specifically, this lower boundary condition can be presented as follows:
(32)m4x=xf=m2x=xf,wf2<y<l1,m3,jx=xf,lj−1<y<lj.Their dimensionless formulas with the consideration of stress sensitivity effects can be as follows:(33)∂2η4D∂xD2=ω4dk4D∂η4D∂tD,1<xD<xeD,with the following upper and lower boundary conditions:
(34)∂η4D∂xDxD=xeD=0,γ4Dη4DxD=1=γ2Dη2DxD=1,wfD2<yD<l1,D,γ3,jDη3,jDxD=1,lj−1,D<yD<lj,D.Laplace transformation is used to transform these equations to the Laplace domain:(35)d2η4D¯dxD2=sω4dk4Dη4D¯,1<xD<xeD,(36)dη4D¯dxDxD=xeD=0,γ4Dη4D¯xD=1=γ2Dη2D¯xD=1,wfD/2<yD<l1,D,γ3,jDη3,jD¯xD=1,lj−1,D<yD<lj,D.


## 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:
(37)η4D¯=A4esω4d/k4DxD+B4e−sω4d/k4DxD.The upper boundary condition can derive the relationship betweenA4 and B4:
(38)B4=A4e2sω4d/k4DxeD.Two lower boundary conditions as shown in Equation (31) could get different solutions:
(39)η4D¯=γ2Dcoshsω4d/k4DxD−xeDγ4Dcoshsω4d/k4D1−xeDη2D¯,wfD2<yD<l1,D,γ3jDcoshsω4d/k4DxD−xeDγ4Dcoshsω4d/k4D1−xeDη3jD¯,lj−1,D<yD<lj,D.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:
(40)ZoneII:d2η2D¯dyD2=γ2Dk4Dγ4Dk2Dsω4Dk4Dtanhsω4dk4D1−xeD+sω2Dk2Dη2D¯,wfD2<yD<l1,D,ZoneIII:d2η3,jD¯dyD2=γ3jDk4Dγ4Dk3,jDsω4Dk4Dtanhsω4dk4D1−xeD+λω3mDsλ+ω3mDsk3,jD+sω3,jDk3,jDη3,jD¯,lj,D<yD<lj+1,D.To simplify the notation,(41)β=γ2Dk4Dγ2Dk3Dsω4Dk4Dtanhsω4dk4D1−xeD+sω2Dk2D,αj=γ3jDk4Dγ4Dk2,jDsω4Dk4Dtanhsω4dk4D1−xeD+λω3mDsλ+ω3mDsk3,jD+sω3,jDk3,jD,j=1,2,⋯,N.Similarly, the general solutions for Equations (38) and (39) can be described as follows:
(42)ZoneII:η2D¯=A2eβxD+B2e−βxD,wfD2<yD<l1,D,ZoneIII:η3,jD¯=A3,jeαjxD+B3,je−αjxD,lj,D<yD<lj+1,D.There are in total2N+1 coefficients that need to be determined. That is to say, we should write 2N+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,(43)γ3,jDA3,jeαjlj,D+B3,je−αjlj,D=γ3,j+1DA3,j+1eαj+1lj,D+B3,j+1e−αj+1lj,D,k3,jDαjA3,jeαjlj,D−αjB3,je−αjlj,D=k3,j+1Dαj+1A3,j+1eαj+1lj,D−αj+1B3,j+1e−αj+1lj,D.At the interface of Zone II and Zone III,(44)γ3,1DA3,1eα1l1,D+B3,1e−α1l1,D=γ2DA2eβl1,D+B2e−βl1,D,k3,1Dα1A3,1eα1l1,D−α1B3,1e−α1l1,D=k2DβA2eβl1,D−βB2e−βl1,D.Two other equations can be further obtained at the inner boundary of Zone II and outer boundary of Zone III:(45)γ1Dη1D¯=γ2DA2eα1wfD/2+B2e−α1wfD/2,A3,NeβyeD−B3,Ne−βyeD=0.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:
(46)d2η1D¯dxD2=sω1Dk1D−2k2DC1DβA2eβwfD/2−B2e−βwfD/2η1D¯.The general solutions should be(47)η1D¯=A1eχxD+B1e−χxD,0<xD<1where χ=sω1D/k1D−2k2D/C1DβA2eβwfD/2−B2e−βwfD/2.Conditioning to the inner and outer boundary conditions, these two coefficients are as follows:(48)A1=πsC1Dχe2χ−1,B1=πe2χsC1Dχe2χ−1.Finally, we can obtain the bottom-hole pressure:(49)ηwD¯=η1D¯xD=0=πe2χ+1sC1Dχe2χ−1.Following the methods used by Xu et al. [25], the zero-order perturbation solution of the bottom-hole pressure in the Laplace space considering the wellbore storage CD and the skin factor S is obtained:
(50)η¯wDS,CD=sη¯wD+Ss+CDs2sη¯wD+S.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 mwD can be obtained by Equation (9).


## 3. 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.Table 1
Basic data used for the model.
Parameters
Value
Units
Reference lengthLref
0.1
m
Hydraulic fracture permeability in Zone I,k1
1000
mD
Permeability in SRV Zone II,k2
0.1
mD
Permeability in SRV Zone III,k3
0.01
mD
Permeability in Zone IV,k4
0.001
mD
Initial formation pressure,pi
30
MPa
Fluid viscosity
2.7×10‐3
mPa·s
Compressibility in Zone I,Ct1
5×10‐4
MPa-1
Compressibility in SRV Zone II,Ct2
5×10‐4
MPa-1
Compressibility in SRV Zone III,Ct3
5×10‐4
MPa-1
Compressibility in Zone IV,Ct4
5×10‐4
MPa-1
Porosity in Zone I
0.15
Fraction
Porosity in SRV Zone II
0.15
Fraction
Porosity in SRV Zone III
0.15
Fraction
Porosity in Zone IV
0.15
Fraction
Half-length of the fractures,xf
100
m
Width of the SRV zone,ye
100
m
Length of the SRV zone,xe
150
mIn 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. Figure3 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 proposed models can be generalized to be any multilinear flow models (MLFM).Figure 3
Investigation of number of subdivisions of SRV Zone III on the type curves.


## 4. Results and Discussions

### 4.1. Transient Pressure Behavior Analysis

Figure4 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 Lref as 0.1 m. The values of the relevant dimensionless variables used in this numerical experiments are as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05,CD=0.001, and S=0.1. We can approximately observe five flow stages from this log-log of dimensionless pseudopressure and its derivative curves.Figure 4
Typical curves of multifractured horizontal well in stress-sensitive SGR formation.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–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 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.


### 4.2. 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 stress-sensitive SGR.


#### 4.2.1. Effect of Stress Sensitivity

Figure5 shows the effect of stress sensitivity on the pressure transient curves. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, CD=0.001, and S=0.1. This newly proposed multilinear flow model can be used to analyze two situations: (a) there is a constant permeability module for all zones; (b) each zone has different permeability modules. We will separately investigate these two situations.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.
(a)
(b)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 Figure5(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 Figure6. The values of relevant parameters are listed as follows: k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05,, CD=0.001, and S=0.1. By analysis, the mobility capacities of Zone I and SRV Zone II are represented by dimensionless fracture conductivity CfD and dimensionless formation permeability k2D. 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 CfD 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 CfD 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.Figure 6
The effect of mobility capacity of Zone I and the SRV Zone II on pressure and derivative curves.
(a)
Zone I
(b)
SRV Zone IIThe mobility capacityk2D determines gas flow capacity in the SRV Zone II. We set the mobility capacity k2D 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 k2D 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 k2D decreases. This means that a large value of k2D 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 k2D increases. Our proposed multilinear flow model will show a 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.


#### 4.2.3. 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 Figure7. The values of relevant parameters are listed as follows: CfD=10, k2D=100, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=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 k3D 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.Figure 7
The effect of mobility capacity of Zone III and Zone IV on pressure and derivative curves.
(a)
Zone III
(b)
Zone IVThe mobility capacityk4D determines gas flow capacity in the outer Zone IV. We set mobility capacity k4D 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 k4D has great influence on the formation bilinear flow stage in SRV Zone II. This is very different from the effects of mobility capacity k3D 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 k4D decreases, the gas supply will be reduced as well; as a result, the dimensionless pressure curves will move upward.


#### 4.2.4. Effect of the Size of the Outer Region

Three cases with different sizes of outer regions shown in Figure8 are studied in this section to identify the effect of size of outer regions on the transient behaviors. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1. As shown in Figure 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.Figure 8
The effect of size of the outer region on pressure and derivative curves.
(a)
Length of the outer region
(b)
Width of the outer region


#### 4.2.5. Effect of the Size of SRV Zone II

Figure9 illustrates the effect of the size of SRV Zone II on the transient behaviors. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1, and the dimensionless size of SRV Zone II LD 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.Figure 9
The effect of the size of the SRV Zone II on pressure and derivative curves.


#### 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 intoN small zones as illustrated in Figure 1(b). Figure 10 shows the effect of the coefficient of permeability variation in SRV Zone III on the transient behavior; the values of the parameters related to the Figure 10 are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1.Figure 10
The effect of coefficient of permeability variation on pressure and derivative curves.
(a)
Linear function
(b)
Logarithmic functions
(c)
Linear function
(d)
Logarithmic functions
(e)
Linear function
(f)
Logarithmic functionsTo 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. Figures10(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.Figures10(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 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 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 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.


## 4.1. Transient Pressure Behavior Analysis

Figure4 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 Lref as 0.1 m. The values of the relevant dimensionless variables used in this numerical experiments are as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05,CD=0.001, and S=0.1. We can approximately observe five flow stages from this log-log of dimensionless pseudopressure and its derivative curves.Figure 4
Typical curves of multifractured horizontal well in stress-sensitive SGR formation.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–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 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.


## 4.2. 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 stress-sensitive SGR.


### 4.2.1. Effect of Stress Sensitivity

Figure5 shows the effect of stress sensitivity on the pressure transient curves. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, CD=0.001, and S=0.1. This newly proposed multilinear flow model can be used to analyze two situations: (a) there is a constant permeability module for all zones; (b) each zone has different permeability modules. We will separately investigate these two situations.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.
(a)
(b)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 Figure5(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 Figure6. The values of relevant parameters are listed as follows: k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05,, CD=0.001, and S=0.1. By analysis, the mobility capacities of Zone I and SRV Zone II are represented by dimensionless fracture conductivity CfD and dimensionless formation permeability k2D. 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 CfD 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 CfD 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.Figure 6
The effect of mobility capacity of Zone I and the SRV Zone II on pressure and derivative curves.
(a)
Zone I
(b)
SRV Zone IIThe mobility capacityk2D determines gas flow capacity in the SRV Zone II. We set the mobility capacity k2D 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 k2D 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 k2D decreases. This means that a large value of k2D 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 k2D increases. Our proposed multilinear flow model will show a 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.


### 4.2.3. 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 Figure7. The values of relevant parameters are listed as follows: CfD=10, k2D=100, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=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 k3D 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.Figure 7
The effect of mobility capacity of Zone III and Zone IV on pressure and derivative curves.
(a)
Zone III
(b)
Zone IVThe mobility capacityk4D determines gas flow capacity in the outer Zone IV. We set mobility capacity k4D 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 k4D has great influence on the formation bilinear flow stage in SRV Zone II. This is very different from the effects of mobility capacity k3D 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 k4D decreases, the gas supply will be reduced as well; as a result, the dimensionless pressure curves will move upward.


### 4.2.4. Effect of the Size of the Outer Region

Three cases with different sizes of outer regions shown in Figure8 are studied in this section to identify the effect of size of outer regions on the transient behaviors. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1. As shown in Figure 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.Figure 8
The effect of size of the outer region on pressure and derivative curves.
(a)
Length of the outer region
(b)
Width of the outer region


### 4.2.5. Effect of the Size of SRV Zone II

Figure9 illustrates the effect of the size of SRV Zone II on the transient behaviors. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1, and the dimensionless size of SRV Zone II LD 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.Figure 9
The effect of the size of the SRV Zone II on pressure and derivative curves.


### 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 intoN small zones as illustrated in Figure 1(b). Figure 10 shows the effect of the coefficient of permeability variation in SRV Zone III on the transient behavior; the values of the parameters related to the Figure 10 are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1.Figure 10
The effect of coefficient of permeability variation on pressure and derivative curves.
(a)
Linear function
(b)
Logarithmic functions
(c)
Linear function
(d)
Logarithmic functions
(e)
Linear function
(f)
Logarithmic functionsTo 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. Figures10(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.Figures10(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 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 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 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.


## 4.2.1. Effect of Stress Sensitivity

Figure5 shows the effect of stress sensitivity on the pressure transient curves. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, CD=0.001, and S=0.1. This newly proposed multilinear flow model can be used to analyze two situations: (a) there is a constant permeability module for all zones; (b) each zone has different permeability modules. We will separately investigate these two situations.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.
(a)
(b)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 Figure5(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 Figure6. The values of relevant parameters are listed as follows: k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05,, CD=0.001, and S=0.1. By analysis, the mobility capacities of Zone I and SRV Zone II are represented by dimensionless fracture conductivity CfD and dimensionless formation permeability k2D. 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 CfD 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 CfD 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.Figure 6
The effect of mobility capacity of Zone I and the SRV Zone II on pressure and derivative curves.
(a)
Zone I
(b)
SRV Zone IIThe mobility capacityk2D determines gas flow capacity in the SRV Zone II. We set the mobility capacity k2D 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 k2D 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 k2D decreases. This means that a large value of k2D 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 k2D increases. Our proposed multilinear flow model will show a 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.


## 4.2.3. 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 Figure7. The values of relevant parameters are listed as follows: CfD=10, k2D=100, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=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 k3D 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.Figure 7
The effect of mobility capacity of Zone III and Zone IV on pressure and derivative curves.
(a)
Zone III
(b)
Zone IVThe mobility capacityk4D determines gas flow capacity in the outer Zone IV. We set mobility capacity k4D 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 k4D has great influence on the formation bilinear flow stage in SRV Zone II. This is very different from the effects of mobility capacity k3D 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 k4D decreases, the gas supply will be reduced as well; as a result, the dimensionless pressure curves will move upward.


## 4.2.4. Effect of the Size of the Outer Region

Three cases with different sizes of outer regions shown in Figure8 are studied in this section to identify the effect of size of outer regions on the transient behaviors. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1. As shown in Figure 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.Figure 8
The effect of size of the outer region on pressure and derivative curves.
(a)
Length of the outer region
(b)
Width of the outer region


## 4.2.5. Effect of the Size of SRV Zone II

Figure9 illustrates the effect of the size of SRV Zone II on the transient behaviors. The values of relevant parameters are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1, and the dimensionless size of SRV Zone II LD 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.Figure 9
The effect of the size of the SRV Zone II on pressure and derivative curves.


## 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 intoN small zones as illustrated in Figure 1(b). Figure 10 shows the effect of the coefficient of permeability variation in SRV Zone III on the transient behavior; the values of the parameters related to the Figure 10 are listed as follows: CfD=10, k2D=100, k3D=10, k4D=1, ω3D=0.01, λ=0.15, mDL=0.12, VDL=0.1, γD=0.05, CD=0.001, and S=0.1.Figure 10
The effect of coefficient of permeability variation on pressure and derivative curves.
(a)
Linear function
(b)
Logarithmic functions
(c)
Linear function
(d)
Logarithmic functions
(e)
Linear function
(f)
Logarithmic functionsTo 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. Figures10(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.Figures10(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 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 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 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.


## 5. 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
(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
(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



---

*Source: 8348205-2020-02-07.xml*