Pressure Transient Analysis and Flux Distribution for Multistage Fractured Horizontal Wells in Triple-Porosity Reservoir Media with Consideration of Stress-Sensitivity Effect

Triple-porosity model is usually adopted to describe reservoirs with multiscaled pore spaces, including matrix pores, natural fractures, and vugs. Multiple fractures created by hydraulic fracturing can effectively improve the connectivity between existing natural fractures and thus increase well deliverability. However, little work has been done on pressure transient behavior of multistage fractured horizontal wells in triple-porosity reservoirs. Based on source/sink function method, this paper presents a triple-porosity model to investigate the transient pressure dynamics and flux distribution for multistage fractured horizontal wells in fractured-vuggy reservoirs with consideration of stress-dependent natural fracture permeability. The model is semianalytically solved by discretizing hydraulic fractures and Pedrosa’s transformation, perturbation theory, and integration transformation method. Type curves of transient pressure dynamics are generated, and flux distribution among hydraulic fractures for a fractured horizontal well with constant production rate is also discussed. Parametric study shows that major influential parameters on transient pressure responses are parameters pertinent to reservoir properties, interporosity mass transfer, and hydraulic fractures. Analysis of flux distribution indicates that flux density gradually increases from the horizontal wellbore to fracture tips, and the flux contribution of outermost fractures is higher than that of inner fractures. The model can also be extended to optimize hydraulic fracture parameters.


Introduction
The concept of dual-porosity media was first proposed by Barenblatt et al. [1] to describe naturally fractured media.Warren and Root [2] then extended this model to analyze pressure transient dynamics for vertical wells in naturally fractured oil reservoirs.Recently, Cai et al. [3,4] analyzed the imbibition mechanism in fractured-porous media based on fractal geometry.Although the dual-porosity model can describe most naturally fractured reservoirs, there are still some drawbacks when applying the dual-porosity model to describe reservoirs with multiple pore media, such as reservoirs with natural fractures and vugs (different from matrix pores).
In order to better describe reservoir heterogeneity in this type of naturally fractured-vuggy reservoirs, the concept of triple-porosity model was proposed.The first triple-porosity model in the literature was proposed by Abadassah and Ershaghi [5], in which the reservoir was represented by a combination of two matrix systems with different properties and uniformly distributed natural fractures.The model was then further extended by Al-Ghamdi and Ershaghi [6] to represent reservoirs with two fracture systems (microfractures and macrofractures) and one matrix system.After that, over the past several decades, numerous efforts have been made to investigate transient pressure behavior in so-called triple-porosity reservoirs.However, almost all the research focuses on vertical wells [7][8][9][10][11][12][13][14][15][16][17][18][19] or horizontal wells [20][21][22] in triple-porosity reservoirs, and few studies have been done on complex well-reservoir configuration, such as multistage fractured horizontal well in triple-porosity reservoirs.
Field practices have proven that horizontal wellbore combining multistage hydraulic fracturing can effectively enhance the permeability in the vicinity of horizontal wellbore.Almost all the triple-porosity models [23,24] for fractured horizontal wells assume linear flow, which cannot completely characterize the pressure responses during all flowing periods (such as pseudoradial flow period).On the other hand, in most existing models for fractured horizontal wells in triple-porosity reservoirs three porous media denote hydraulic fractures, natural fractures, and matrix pores, which means the reservoir itself is actually considered as a dual-porosity medium.
This study develops a semianalytical model based on source/sink function theory for analyzing transient pressure responses and flux distribution of naturally fractured-vuggy reservoirs, in which the reservoir itself is conceptualized as triple-porosity media and the interaction between hydraulic fractures and reservoir is considered.The model presented here can completely reflect transient pressure characteristics during all the possible flowing period as well as flux distribution among multiple fractures.In addition, this model also takes into account the stress-sensitivity of permeability caused by closure of natural fractures during production.

Physical Model of a Multistage Fractured Horizontal Well in a Triple-Porosity Reservoir
As shown in Figure 1, a horizontal well with multiple hydraulic fractures, producing at a constant rate () in a reservoir with multiscale storage spaces, is considered here.Both the cap rock and underlying rock formation of the reservoir are usually shales with extremely low permeability, which can be conceived as no-flow boundaries.The reservoir thickness and oil viscosity are assumed to be ℎ and , separately.Other basic assumptions involved in the derivation of the triple-porosity model are as follows: (1) The reservoir is assumed to be composed of three contiguous porous media, which are matrix, natural fractures, and vugs (representing larger voids in the reservoir which are not part of natural fractures).
(2) Natural fractures are assumed to be directly connected to hydraulic fractures or the horizontal wellbore, while the fluid stored in matrix pores or vugs does not have direct access to hydraulic fractures or the horizontal wellbore.Figure 2 presents a simple illustration of the fluid flowing path in triple-porosity reservoirs.
(3) The permeability of natural fractures is considered stress-sensitive to incorporate the effect of partial or complete closure of natural fractures during production.
(4) Both oil and rock are considered slightly compressible, and the reservoir has uniform pressure distribution at time zero.
(5) Consider single-phase isothermal flow and negligible gravity and capillary effects.
To obtain the pressure responses caused by the wellreservoir configuration shown in Figure 1, a mathematical model together with corresponding boundary and initial conditions is required.In the problem addressed here, the multiple fractures inner boundary condition is so complex that it is rather difficult to directly describe it with a mathematical equation.In order to solve this problem, source function theory, a powerful tool to solve complex transient flow problems in reservoirs, is adopted in this work.We first studied the transient pressure responses caused by a continuous line-sink in triple-porosity reservoirs and then adopted the line-sink solution with superposition principle to obtain the pressure transient dynamics and flux distribution of multistage fractured horizontal wells in triple-porosity reservoirs.

Line-Sink Model in Stress-Sensitive
Triple-Porosity Reservoirs  3, a continuous line-sink in a triple-porosity reservoir with stress-dependent fracture permeability is considered in this section.The strength of the line-sink is represented by q().

Mathematical Model for a Line-Sink in Triple-Porosity Reservoirs
(1) Governing Equations.Imposing mass-balance equation over a control volume in the triple-porosity reservoir and substituting equation of state and equation of motion, we can get the governing equations for matrix, natural fracture, and vug systems.For natural fracture system with stress-dependent permeability, the governing equation is as follows (i.e., (A.26b) in Appendix A): For matrix system, the governing equation is as follows (i.e., (A.30) in Appendix A): For vug system, the governing equation is as follows (i.e., (A.34) in Appendix A): (2) Boundary Conditions.According to the derivation in Appendix B, the inner boundary condition of the line-sink model is (i.e., (B.5) in Appendix B) lim Assuming infinitely large reservoirs, the outer boundary condition can be expressed as (3) Initial Condition.The initial pressure distribution in the reservoir is assumed to be uniform, which is With the definitions in Table 1, the dimensionless forms of (1)∼( 6) can be obtained as follows: lim 3.4.Solution of the Dimensionless Line-Sink Model 3.4.1.Pedrosa's Linearization.The stress-sensitivity of natural fracture permeability makes the above seepage model strongly nonlinear which cannot be solved analytically.Here, following Pedrosa Jr. [25], we introduced the following equation to alleviate the nonlinearity of ( 8) and (11): where  D is an intermediate dimensionless parameter which is also a function of radial distance and time.With ( 14), ( 8)∼(13) become lim  15) through ( 20) can be expanded as power series in dimensionless permeability modulus, and we can get the following equations: Given that the dimensionless permeability modulus  D is usually small (much smaller than 1), the zero-order approximate solution can satisfy the requirements of engineering precision, so (15)∼(20) become Equation ( 22) is the linearized line-sink model in tripleporosity reservoirs.

Laplace Transformation. Take the following Laplace transformation:
where  is the Laplace transformation variable.
Then (22) becomes lim 3.4.4.Line-Sink Solution.Equations ( 25) and ( 26) can be rewritten as Substituting ( 29) into (24) yields If we define then (30) can be expressed as follows: The general solution of (32) can be easily obtained as follows: where  0 and  0 are Bessel functions and  1 and  2 are unknown coefficients which can be determined by corresponding boundary conditions.If the line-sink is located at the origin of coordinates (center of the triple-porosity reservoir), then the dimensionless radial distance  D in (31) can be calculated by the following equation: If the line-sink is located at ( w ,  w ) instead, then  D is calculated by the following equation: where Substitution of (33) into the inner boundary condition ( 27) yields lim According to the properties of Bessel's function, which are lim  → 0  1 () → 1 and lim  → 0  1 () → 0, (39) becomes Similarly, with the outer boundary condition given in ( 28) and the properties of Bessel's function, lim  → ∞  0 () → 0 and lim  → ∞  0 () → ∞, the coefficient in (33),  2 , should satisfy the following equation: With ( 40) and (41), we can get the final form of (33) as follows: Equation ( 42) is the basic zero-order perturbation solution for a line-sink in a stress-sensitive triple-porosity reservoir.As mentioned above, zero-order perturbation solution is enough to approximate the exact solution of ( 15) to (20); that is,  D ≈  D0 .With the definition of qD in Table 1, (42) can be expressed as Equation ( 43) can be rewritten as where Equation ( 44) is the basic line-sink solution for stresssensitive triple-porosity reservoirs.With the basic line-sink solution and the superposition principle, we can obtain the pressure responses caused by multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs.
x y

Pressure Responses for Multistage Fractured Horizontal Wells in Stress-Sensitive Triple-Porosity Reservoirs
Figure 4 shows the schematic of a horizontal well with multiple (a total of ) hydraulic fractures in a stress-sensitive triple-porosity reservoir.The -axis is set to be aligned with the horizontal wellbore.The intersection of the th ( = 1, 2, . . ., ) hydraulic fracture and the -axis is represented by (0,   ), and the distance between every two adjacent fractures is Δ  .
To successfully obtain the pressure responses caused by the multiple fractures with the line-sink solution derived in the above section, it is necessary to discretize the multiple hydraulic fractures.As shown in Figure 4, each hydraulic fracture is divided into 2 segments along its length.The coordinate of the midpoint (i.e., discrete node) of the th segment on the th fracture is denoted by ( , ,  , ), and the coordinates of the two corresponding endpoints are denoted by ( , ,  , ) and ( ,+1 ,  ,+1 ).
With the discretization scheme shown in Figure 4, the coordinates of endpoints of discrete segments can be calculated by the following equations: The coordinates of discrete nodes (midpoint of each segment) can be calculated by the following equations: where  fL is the length of the left wing of th fracture, m;  fR is the length of the right wing of th fracture, m.In the model presented in this paper,  fL and  fR can be different and can vary from fracture to fracture.The distance between   and  −1 is It is obvious that the flux strength is different along the fracture length.However, in the same discrete segment (assuming the discrete segment is small enough), the flux strength can be considered as constant.If we use  , to represent the flux density per unit length, then, with the superposition principle, the pressure response at ( D ,  D ) in the reservoir caused by the discrete segment (, ) can be obtained by integrating the basic line-sink along the discrete segment, which is With the definitions of  D and  wD , that is, (35) and ( 37), (49) can be changed into The pressure response at ( D ,  D ) caused by 2 ×  segments can be obtained by applying the principle of superposition over all discrete fracture segments: Thus, the pressure response at the discrete segment (, ) (1 ≤  ≤ , 1 ≤  ≤ 2) can be obtained as follows: Since the permeability of hydraulic fractures is always much higher than the original reservoir permeability, the pressure drop in hydraulic fractures is much smaller than the pressure drop caused by oil flow in the reservoir.Thus, the hydraulic fractures can be considered as infinitely conductive, meaning the pressure in hydraulic fractures is equal to the bottom-hole pressure: Combining ( 55) and (56), we can get The subscripts  and  in (57) denote that this equation is written for discrete segment (, ).By writing (57) for all discrete segments, that is, letting ( = 1, 2, . . ., ,  = 1, 2, . . ., 2), we can get 2× equations.It should be noted that there are 2×+1 unknowns in these 2× equations, which are  wD and  D, ( = 1, 2, . . ., ,  = 1, 2, . . ., 2).To compose a closed system of equations, the assumption of constant total flow rate can give the following equation: where Δ , is the length of the discrete segment (, ) and can be calculated by Taking Laplace transformation of (58), one can get The dimensionless form of (60) is where Equations ( 57) and (61) represent a system of (2×+1) linear algebraic equations relating (2 ×  + 1) unknowns, which can be solved by direct methods, such as Gaussian elimination method, Gauss-Jordan reduction method.
wD calculated by (57) and (61) does not take into account the effects of wellbore storage and skin.According to van Everdingen [26] and Kucuk and Ayestaran [27], the effects of wellbore storage and skin can be incorporated in the calculated bottom-hole pressure response by where  D is the dimensionless wellbore storage coefficient and  is the skin factor defined by the following, respectively: where  is the wellbore storage coefficient, m 3 /Pa;  sf is the sandface flow rate of the multistage fractured horizontal well in the triple-porosity reservoir, m 3 /s; Δ  is an extra pressure drop near the wellbore, Pa. wDS in (63) is the bottom-hole pressure response obtained in the Laplace domain, and by Stehfest [28] numerical inversion algorithm we can calculate the pressure responses  wDS ( D ) in real time domain.Then, with (66), we can obtain the bottom-hole pressure response for multistage fractured horizontal wells in triple-porosity reservoirs incorporating the stress-sensitivity of natural fracture system:

Results and Discussion
In this section, a horizontal well with three hydraulic fractures (i.e.,  = 3) in a triple-porosity reservoir is considered, and the model introduced above is applied to investigate the dimensionless pressure responses of this well-reservoir configuration.Period 1.It is the early-time wellbore storage period.In this period, both pressure curves and their derivative curves exhibit unit slope on the log-log plots.

Flow Periods and Effect of Stress-Sensitivity
Period 2. It is the transition flow period after wellbore storage period.The pressure derivative curves exhibit a "hump, " representing the effect of skin.
Period 3. It is the first linear flow period.This period corresponds to the linear flow in the direction perpendicular  to hydraulic fracture length as illustrated in Figure 6(a).The pressure derivative curves exhibit a straight line with slope equal to "1/2." Period 4. It is the first pseudoradial flow period.This period can occur when the distances between adjacent fractures are relatively large.Oil flows in the reservoir in a way as shown in Figure 6(b), and pseudoradial flow occurs around each fracture.This period corresponds to a nearly horizontal line with a value of "1/(2)" in the derivative curve, where  is the number of hydraulic fractures.
Period 5.It is the second linear flow period.This period corresponds to linear flow perpendicular to the horizontal wellbore as shown in Figure 6(c).Another straight line with a slope of about "1/2" can be observed in the derivative curve.
When the permeability of natural fracture system remains constant during the whole production, corresponding to the no stress-sensitivity case, the slope of the derivative curves during this period is "1/2." When taking into account the stress-sensitivity of natural fracture system, the corresponding pressure derivative curve slightly deviates from the onehalf-slope line and exhibits a slight upward tendency.
Period 6.It is the second pseudoradial flow period, corresponding to radial flow in natural fracture system.During this period, multiple hydraulic fractures and the horizontal wellbore behave like an enlarged wellbore.The derivative curves exhibit a horizontal line with an upward tendency depending on the value of dimensionless permeability modulus.In the case that  D = 0, the derivative curves are horizontal with a value of 0.5 on the -axis.With the increase of the value of  D , the derivative curves turn upward gradually.
Period 7. It is the first interporosity flow period.This period reflects mass transfer between vugs and natural fracture system.Due to the oil flowing from vugs to natural fractures, the bottom-hole pressure drops slower and an obvious "dip" can be observed in the derivative curves.
Period 8.It is the third pseudoradial flow period, corresponding to compound radial flow in natural fracture system and vugs.Due to the existence of stress-sensitivity of natural fracture system, the pressure derivative curves are no longer horizontal.The position of the derivative curves becomes higher with larger value of  D .
Period 9.It is the second interporosity flow period.This period reflects oil flowing from matrix to natural fracture system.Another characteristic "dip" can be observed in the derivative curves.
Period 10.It is the late-time compound pseudoradial flow period.This period reflects compound radial flow in natural fracture system, matrix, and vugs as shown in Figure 6(d).
Similarly, the derivative curves exhibit upward tendencies due to the existence of stress-sensitivity of natural fracture system.

Effect of Skin
Factor. Figure 7 shows the effect of skin factor on dimensionless pressure curves and corresponding derivative curves.As stated above, the skin factor mainly affects the shape of type curves during the transition period (Period 2) after wellbore storage period.With the increase of the value of skin factor, the position of the "hump" in the derivative curves becomes higher, representing larger pressure drop in the formation.

Effect of Storativity Ratio of Natural Fracture System.
Figure 8 shows the effect of storativity ratio of natural fracture system on dimensionless pressure and pressure derivative curves.It is obvious in Figure 8 that when the storativity ratio of matrix (i.e.,  m ) and other parameters are kept constant, the storativity ratio of natural fracture system  f mainly affects the transient pressure dynamics during the first linear flow period, the first pseudoradial flow period, the second linear flow period, second pseudoradial flow period, and the first interporosity flow period.On the one hand, the position of the dimensionless pressure curves becomes higher during the abovementioned flowing periods with the decrease of  f .On the other hand, with the decrease of  f , the position of the pressure derivative curves during the first linear flow period, the first pseudoradial flow, and the second linear flow period becomes higher, and the first "dip" in pressure derivative curves becomes wider and deeper.

Effect of Storativity Ratio of Matrix.
Figure 9 shows the effect of storativity ratio of matrix on dimensionless pressure and derivative curves.It can be observed that, with fixed  f , the storativity ratio of matrix  m has primary effect on pressure dynamics during two interporosity flowing periods.
With the decrease of the value of  m , the fist "dip" in the pressure derivative curves becomes wider and deeper, while the second "dip" in the pressure derivative curves exhibits the opposite trend.interporosity flow period, that is, the appearance time of the first "dip" in pressure derivative curves.The smaller the value of  vf is, the later the first "dip" can be observed in the pressure derivative curves.When  vf decreases to a certain extent, the third pseudoradial flow period (Period 8 in Figure 5) may not be observed in the type curves.In addition, when  vf increases to a certain extent, the second pseudoradial flow

Effect of Interporosity Flow Coefficient between Vugs and Natural Fractures.
Figure 9: Type curves with different storativity ratios of matrix.period (Period 6 in Figure 5) and even the second linear flow period (Period 5 in Figure 5) may be masked by the following first interporosity flow period.

Effect of Interporosity Flow Coefficient between Matrix and
Natural Fractures.Figure 11 shows the effect of interporosity flow coefficient between matrix and natural fractures ( mf ) on dimensionless pressure and pressure derivative curves.It can be seen that  vf mainly affects the occurrence time of the second interporosity flow period, that is, the appearance time of the second "dip" in pressure derivative curves.The smaller the value of  mf is, the later the second "dip" can be observed in the pressure derivative curves.With the increase of  mf , the reflection of the third pseudoradial flow period on the derivative curves (Period 8 in Figure 5) may not be observed.When  mf increases to a certain degree, the second "dip" may merge with the first "dip, " meaning that simultaneous interporosity mass transfers happen between natural fractures and matrix as well as natural fractures and vugs.

Effect of Length of Hydraulic
Fractures.Figure 12 presents the effect of half-length of hydraulic fractures on type curves.
For the convenience of discussion,  fL =  fR =  f is assumed here.It can be observed that the length of hydraulic fracture has primary effect on the first linear flow period and the subsequent first pseudoradial flow period.With the increase of fracture length, the first linear flow period lasts longer with lower position of pressure derivative curves, and the first pseudoradial flow period occurs later with shorter duration.When the half-length of hydraulic fractures continues to increase, the horizontal line in the pressure derivative curves reflecting the first pseudoradial flow period around each fracture will gradually disappear.
5.8.Effect of Hydraulic Fracture Spacing. Figure 13 shows the effect of the distance between adjacent hydraulic fractures.The case discussed in Figure 13 assumes equal fracture spacing for the convenience of discussion.It can be observed that the distance between adjacent hydraulic fractures Δ  mainly affects the first pseudoradial flow period and the subsequent second linear flow period.The larger the distance between adjacent fractures, the longer the duration of the first pseudoradial flow period, and the later the occurrence of the second linear flow period.Another point should be addressed here is that when Δ  increases to a certain degree, the second pseudoradial flow period may not be observed in the type curves; that is, the first interporosity flow period appears after the second linear flow period.5.9.Flux Distribution.Figure 14 presents the dimensionless flux distribution along hydraulic fractures at a given time based on the calculation results obtained with the model proposed in this paper.In the case presented in Figure 14, three hydraulic fractures with equal half-length and fracture spacing are considered, and each hydraulic fracture is discretized into ten segments.It can be observed that, for the same hydraulic fracture, the flux distribution is symmetrical with respect to the horizontal wellbore, and the flux density at fracture tips is larger than that in the middle of the fracture.In addition, for different hydraulic fractures, it can be found that the flux contribution from hydraulic fracture in the middle (Fracture 3 in Figure 14) is always smaller than that from outer fractures (Fractures 1 and 2 in Figure 14).

Conclusion
Based on the mathematical model and discussion in this paper, the following conclusions can be warranted: (1) A mathematical model is presented for investigating transient pressure dynamics as well as flux distribution of multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs.The model can better describe fluid flow in naturally fractured reservoirs with multiple types of storage spaces and complex well type.(2) Semianalytical solution is developed in the Laplace domain by the integral transformation method and discretization of hydraulic fractures.With the semianalytical solution and the numerical inversion algorithm proposed by Stehfest, dimensionless pressure responses and corresponding pressure derivatives as well as flux distribution among all the fractures can be obtained in real time domain.Computation results prove that the model presented in this paper is stable and accurate.
(3) Analysis of transient pressure and derivative responses indicates that as many as ten flow periods can be identified for a multistage fractured horizontal well in triple-porosity reservoirs.Major factors affecting the transient pressure responses include skin factor, storativity ratios of different porous media type, interporosity flow coefficients, and parameters related to the geometry and placement of hydraulic fractures.With different combinations of these parameters, the reflection of some flow periods may be masked in the pressure derivative curves.
(4) Two characteristic "dips" can be observed in the pressure derivative curves for triple-porosity reservoirs, reflecting the mass transfer between natural fracturevugs and natural fracture-matrix, separately.
(5) The stress-sensitivity of natural fracture system results in larger pressure drop during intermediate and late flowing periods which is reflected by upward tendencies in both dimensionless pressure curves and corresponding derivative curves.
(6) Analysis of flux distribution among multiple hydraulic fractures indicates that the flux contribution from outer fractures is higher, and for a given fracture the flux contribution from fracture tips is higher.(7) The presented model can be used to interpret pressure signal for fractured horizontal wells in triple-porosity reservoirs and provide some important dynamic parameters for reservoir development.

A. Derivation of Governing Equations for Natural Fracture System, Matrix System, and Vug System in Triple-Porosity Reservoirs
Governing equations for natural fracture system, matrix system, and vug system in triple-porosity reservoirs can be obtained by combining equation of motion, equation of state, and mass conservation equation.
(1) Equation of Motion.Assuming radial flow in triple-porosity reservoirs, the oil flow velocity in natural fracture system is given by where V f is the radial oil velocity in natural fracture system, m/s;  f is the permeability of natural fracture system under pressure  f , m 2 ;  is the oil viscosity, Pa⋅s;  f is the pressure of natural fracture system, Pa;  is the radial distance, m.During the production process, the reservoir pressure gradually decreases and the effective stress increases.Consequently, the aperture of natural fractures gradually decreases, resulting in the reduction in the permeability of natural fracture system, which is so-called stress-sensitivity of permeability of natural fracture system.To account for this, a stress-dependent permeability is adopted in the natural fracture flow model.
Following Kikani and Pedrosa [29], the permeability modulus , which is used to describe the stress-sensitivity of natural fracture system, is defined by where  is the permeability modulus, Pa −1 .Equation (A.2) can be further written as where  fi is the permeability of natural fracture system under initial condition, m 2 ;  i is the reservoir pressure under initial condition, Pa.Solving (A.3), one can get  f =  fi e −( i − f ) .(A.4) Equation (A.4) is one of the most common relationships which are used to describe stress-dependent permeability.This exponential relationship corresponds to Type I rocks discussed by Raghavan and Chin [30].
Substituting (A.4) into (A.1)yields the equation of motion in natural fracture system with consideration of stresssensitivity effect as follows: fi e −( i − f )   f  . (A.5)

B. Derivation of Inner Boundary Condition for the Line-Sink
A vertical line-sink can be treated as the limiting case of a vertical cylinder-sink with radius  → 0. The flow rate across the cylinder wall with radius  ( → 0) can be expressed as lim  → 0 V f       = = q () , (B.1 where  is the lateral area of the cylinder-sink, m 2 , which can be calculated by (B.2); V f is oil velocity at the cylinder wall, m/s, which can be calculated by (B.3); q() is flow rate of the continuous line-sink under standard condition, m 3 /s;  is the oil formation volume factor, dimensionless:  Reference length, m  h : Length of horizontal wellbore, m : Total number of hydraulic fractures : Number of discretized segments for each fracture  f : Pressure of natural fracture system, Pa  m : Pressure of matrix system, Pa  v : Pressure of vug system, Pa  fD : Dimensionless pressure of natural fracture system, dimensionless  mD : Dimensionless pressure of matrix system, dimensionless  vD : Dimensionless pressure of vug system, dimensionless  i : Initial reservoir pressure, Pa  0 : Reference pressure, Pa Δ s : Extra pressure drop caused by skin effect near the wellbore, Pa q(): Production rate of the line-sink under standard condition, m 3 /s qD : Dimensionless production rate of the line-sink under standard condition, dimensionless : Constant surface production rate of the multistage fractured horizontal well, m 3 /s   : Flux density of the th segment in the th hydraulic fracture, m 3 /(s⋅m)   (): Laplace transformation of

Figure 1 :
Figure 1: Schematic of a multistage fractured horizontal well in a triple-porosity reservoir.

Figure 2 :
Figure 2: Illustration of fluid flowing paths in a triple-porosity reservoir.

Figure 3 :
Figure 3: A line-sink in a stress-sensitive triple-porosity reservoir.
Figure 5   presents the dimensionless pressure responses and corresponding pressure derivatives of a multistage fractured horizontal well in an infinitely large triple-porosity reservoir with consideration of stress-sensitivity of natural fractures.It can be observed from the figure that as many as ten flow periods can be identified at different time scales.Different flow periods and their corresponding characteristics are explained as follows.

Figure 7 :
Figure 7: Type curves with different skin factors.

Figure 8 :
Figure 8: Type curves with different storativity ratios of natural fracture system.

Figure 10 :
Figure 10: Type curves with different  vf .

Figure 11 :
Figure 11: Type curves with different  mf .
Oil formation volume factor, m 3 /sm 3 : Wellbore storage coefficient, m 3 /Pa  D : Dimensionless wellbore storage coefficient, dimensionless  o : Oil compressibility, Pa −1  f : Compressibility of natural fracture system, Pa −1  m : Compressibility of matrix system, Pa −1  v : Compressibility of vug system, Pa −1  ft : Total compressibility of natural fracture system, Pa −1  mt : Totalcompressibilityof matrixsystem, Pa −1  vt : Total compressibility of vug system, Pa −1 ℎ: Reservoir thickness, m  f : Permeability of natural fracture system, m 2  fi : Permeability of natural fracture system under initial condition, m 2  m : Permeabilityofmatrixsystem,m 2  v : Permeability of vug system, m 2  0 (): Modified Bessel function of the second kind, zero order  1 (): Modified Bessel function of the second kind, first order  0 (): Modified Bessel function of the first kind, zero order : mathematical model are Δ f =  i − f , Δ m =  i −  m , and Δ v =  i −  v ,and relevant dimensionless variables used in the mathematical model are listed in Table 1.Obviously,  f ,  m , and  v in Table 1 satisfy the following equation: