Study on Deliverability Evaluation of Staged Fractured Horizontal Wells in Tight Oil Reservoirs

: At present, the existing deliverability evaluation models mainly consider the impact of speciﬁc factors on production, and the description of the complex fracture network structure primarily remains at the stage of an ideal dual-pore medium with uniform distribution. However, this cannot reﬂect the actual fracture network structure and ﬂuid ﬂow law of fractured horizontal wells. Thus, in this paper, a non-uniform fracture network structure is proposed considering the inﬂuence of the threshold pressure gradient and stress sensitivity characteristics on the production performance of horizontal wells. The stress sensitivity and the fractal theory are combined to characterize the permeability of the complex fracture network, and a three-zone compound unsteady deliverability model for staged fractured horizontal wells in tight oil reservoirs is successfully developed. Laplace transformation, perturbation theory, and numerical inversion are applied to obtain the semi-analytical solution of the proposed deliverability model. The reliability and accuracy of the analytical solution are veriﬁed by the classical tri-linear ﬂow model and an oil ﬁeld example. The effects of related inﬂuential parameters on the production of horizontal wells are analyzed. The deliverability evaluation method proposed in this paper can provide a theoretical basis for formulating rational development technology policies for tight oil reservoirs.


Introduction
Tight oil reserves are abundant and the development of tight oil reservoirs has long been a popular study topic in petroleum research. At present, the technology of "long horizontal section horizontal well + large-scale volume fracturing" is often employed to develop tight oil reservoirs [1,2]. A complex fracture network is formed after large-scale volume fracturing of horizontal wells, which shortens the seepage distance of oil and gas from the matrix to the fracture [3,4], and changes the fluid seepage pattern of the area around the horizontal well, thus improving the development effect of tight oil reservoirs.
Numerous studies have been undertaken globally on the productivity evaluation of fractured horizontal wells in tight oil reservoirs, including steady-state and unsteady-state productivity models. Numerical simulation, analytical, and semi-analytical methods are extensively applied in relevant studies [5][6][7][8][9][10][11][12][13][14]. Li Longlong et al. [15] applied potential theory and the superposition principle to derive the steady-state productivity linear equations of staged multi-cluster fractured horizontal wells in the two cases of infinite and A cluster of primary fractures is first formed during the process of fracturing horizontal wells, and then the primary fractures connect with the natural fractures in the reservoir to generate secondary fractures, thus forming a fracture network around the primary fractures. It has been demonstrated in relevant studies that the density of the fracture network decreases away from the primary fracture. The Fractal Random Fracture Network Generation Algorithm (FRFNA) proposed by Sheng Guanglong et al. [24] based on microseismic data inverted the complex fracture network for horizontal wells in Eagle Ford shows that the further away from hydraulic fractures, the fewer the number of secondary fractures and the lower their density.
Based on this, our paper proposes a non-uniform fracture network structure model, as depicted in Figure 1 (the primary fracture is located at the left side of the model, and the secondary fracture extends to the right direction). As shown in Figure 1, the matrix rock block is unevenly divided by the fracture network, and the fracture network is nonuniformly distributed in the horizontal direction. The farther from the primary fracture, the smaller the fracture density. The gradual change in fracture density along the horizontal well leads to the uneven change in fluid flow capacity. Therefore, to consider the stress sensitivity of the fracture network, the fractal theory is introduced to characterize the complex fluid seepage law under the structure of a non-uniform fracture network.
Energies 2021, 14, x FOR PEER REVIEW 3 of 22 this method, the fractured horizontal well in the tight reservoir is divided into three seepage areas. The fracture network in the volume area between primary fractures is unevenly distributed along the horizontal wellbore, which is different from the classical tri-linear flow model. Laplace transformation, perturbation theory, and numerical inversion are applied to solve the established mathematical model for fractured horizontal wells to obtain the semi-analytical productivity formula. The accuracy of the analytical solution of the model was verified by comparing the analytical results of the simplified model and the classic model. Finally, the reliability of the model was further proved by matching the model calculation results with the actual production data of an oil well in the X Oilfield (in China).

Non-Uniform Fracture Network Structure Model
A cluster of primary fractures is first formed during the process of fracturing horizontal wells, and then the primary fractures connect with the natural fractures in the reservoir to generate secondary fractures, thus forming a fracture network around the primary fractures. It has been demonstrated in relevant studies that the density of the fracture network decreases away from the primary fracture. The Fractal Random Fracture Network Generation Algorithm (FRFNA) proposed by Sheng Guanglong et al. [24] based on micro-seismic data inverted the complex fracture network for horizontal wells in Eagle Ford shows that the further away from hydraulic fractures, the fewer the number of secondary fractures and the lower their density.
Based on this, our paper proposes a non-uniform fracture network structure model, as depicted in Figure 1 (the primary fracture is located at the left side of the model, and the secondary fracture extends to the right direction). As shown in Figure 1, the matrix rock block is unevenly divided by the fracture network, and the fracture network is nonuniformly distributed in the horizontal direction. The farther from the primary fracture, the smaller the fracture density. The gradual change in fracture density along the horizontal well leads to the uneven change in fluid flow capacity. Therefore, to consider the stress sensitivity of the fracture network, the fractal theory is introduced to characterize the complex fluid seepage law under the structure of a non-uniform fracture network.

Seepage Physical Model and Condition Assumptions for Fractured Horizontal Well
According to the fluid seepage characteristics of fractured horizontal wells in tight reservoirs, we established an unsteady productivity seepage physical model for a fractured horizontal well, as shown in Figure 2. The quarter reservoir controlled by a single fracture is viewed as the research object, and the seepage physical model is composed of three parts containing the main fracture area (region I), the volume area between primary fractures (region II), and the unstimulated volume area at the reservoir boundary (region III). Region I represents the primary fracture, which is a linear flow in this area, and the fluid flows directly into the horizontal wellbore along the primary fracture. Region II is the volume area between the primary fractures, which is a dual-porous medium area. There are two kinds of media in zone II: the matrix and the fracture network. The fluid in region II flows in the direction perpendicular to the primary fracture into region I. The complex fracture network in this area is composed of secondary fractures and natural fractures. Fractures are unevenly distributed in zone II. The density of secondary fractures near the primary fracture is the largest and gradually decreases in the y-direction. There are almost no secondary fractures and only a small number natural fractures at the position distant from the primary fracture. Therefore, the reservoir volume distant from the primary fracture can be regarded as the unstimulated reservoir volume. Region III is an unstimulated volume area at the reservoir boundary with a small number of natural fractures. The threshold pressure gradient is introduced to characterize the non-Darcy seepage of the fluid in this area. The fluid in region III flows in the direction parallel to the primary fracture into region II.
According to the fluid seepage characteristics of fractured horizontal wells in tight reservoirs, we established an unsteady productivity seepage physical model for a fractured horizontal well, as shown in Figure 2. The quarter reservoir controlled by a single fracture is viewed as the research object, and the seepage physical model is composed of three parts containing the main fracture area (region I), the volume area between primary fractures (region II), and the unstimulated volume area at the reservoir boundary (region III). Region I represents the primary fracture, which is a linear flow in this area, and the fluid flows directly into the horizontal wellbore along the primary fracture. Region II is the volume area between the primary fractures, which is a dual-porous medium area. There are two kinds of media in zone II: the matrix and the fracture network. The fluid in region II flows in the direction perpendicular to the primary fracture into region I. The complex fracture network in this area is composed of secondary fractures and natural fractures. Fractures are unevenly distributed in zone II. The density of secondary fractures near the primary fracture is the largest and gradually decreases in the y-direction. There are almost no secondary fractures and only a small number natural fractures at the position distant from the primary fracture. Therefore, the reservoir volume distant from the primary fracture can be regarded as the unstimulated reservoir volume. Region III is an unstimulated volume area at the reservoir boundary with a small number of natural fractures. The threshold pressure gradient is introduced to characterize the non-Darcy seepage of the fluid in this area. The fluid in region III flows in the direction parallel to the primary fracture into region II. The basic assumptions of the physical model are as follows: (1) The tight oil reservoir is a homogeneous box-shaped formation, of which the upper and lower boundaries are closed. The horizontal wellbore passes through the center of the reservoir. The half-width of the reservoir is and the thickness of the reservoir is h. (2) The half distance between primary fractures is . (3) The width of the hydraulic fracture is and the half-length of the fracture is . The hydraulic fractures are symmetrically distributed perpendicular to the horizontal wellbore, and the height of the fracture is equal to the thickness of the reservoir. (4) There is a single-phase isothermal unsteady seepage in the reservoir, without considering the influence of gravity and capillary force. (5) Both rock and fluid are slightly compressible and have a constant compression coefficient. (6) The fracture network in region II is a fractal reservoir with a fractal dimension of D, which is embedded in a twodimensional Euclidean rock block. (7) The pressure loss in the horizontal section is negligible. (8) The stress sensitivity of the fracture system and the non-Darcy seepage characteristics of the fluid in region III are taken into consideration. (9) The skin and wellbore storage effects are not considered in the model. (10) The horizontal well produces with the fixed bottom hole pressure.  Considering the stress-sensitive characteristics of the fracture network, the fractal theory is introduced to characterize the permeability and porosity of the complex nonuniform fracture network. It is assumed that the fluid storage space in the fractal dualmedium is at the node with the unit volume of V s , each node has the same volume, the radial distance between the node and the wellbore is x, and the density of the node is N(x). where α-A parameter related to fracture porosity; D-Fractal dimension. There are many calculation methods for the fractal dimension, among which the most commonly used is the box-counting dimension method [25]. This method obtains the box-counting dimension by processing the image and then using the least-squares fitting function. Therefore, the fractal dimension of the fracture network can be obtained by the box-counting dimension method after acquiring the pattern of the fracture network based on the inversion of micro-seismic data [24].
The fracture fractal porosity is defined as φ 2 f . G represents the corresponding symmetry relationship. The equation at (x + dx) is obtained as follows: Substituting Equation (1) into Equation (2), we can get: Rearranging the formula, we can obtain: where d-dimension of Euclidean space embedded by the fractal. When x = x w , the equation is given by: Combining Equations (4) and (5), the fracture fractal porosity is obtained as follows (d equals 2 in two-dimensional space): Considering the pressure sensitivity of the fracture system, its permeability changes with the effective confining pressure and has an exponential relationship with its effective overlying pressure. Referring to the definition of fractal porosity and the abnormal diffusion coefficient, the fracture fractal permeability considering stress sensitivity can be obtained as: The non-uniform fracture network structure model shows that the closer to the primary fracture, the denser the distribution of secondary fractures. The boundary of the hydraulic fracture is taken as the reference point in this paper (the location with the largest fracture density, i.e., x w = w/2), and the value of fractal dimension D is generally 0-2 [22].

Definition of Dimensionless Parameters
The dimensionless seepage mathematical model is obtained by dimensionless treatment of the relevant parameters.

Establishment of the Dimensionless Seepage Mathematical Model
(1) Mathematical Model for Region III Region III is the unstimulated volume area at the reservoir boundary. The fluid in this area flows into the volume area between primary fractures (region II) in the x-direction, which is a one-dimensional linear flow. Considering the influence of the threshold pressure gradient, the continuity equation is obtained as: Incorporating the boundary conditions and initial conditions, we can obtain the dimensionless seepage mathematical model for region III after dimensionless transformation: (2) Mathematical Model for Region II Region II is the volume area between primary fractures. The two kinds of media in this area are the matrix and the fracture network, and the fluid seepage law under the complex fracture network is very complicated. Fluid exchange occurs due to the pressure difference between the matrix and the fracture. A portion of the fluid in this area flows from the matrix to the fracture network, which is regarded as a stable process. The mass of fluid discharged from the matrix to the fracture per unit time is described as follows: where q -Mass of fluid flowing out in per unit rock volume and per unit time, kg/ m 3 ·s , ρ 0 -Density of crude oil, kg/m 3 . Therefore, the governing equation for the matrix system is given by: Secondary fractures are unevenly distributed in the y-direction in this area. The farther away from the primary fracture, the smaller the density of secondary fractures, thus leading to the change in fluid flow capacity in region II in the y-axes. Considering the stress sensitivity, fractal theory based on the self-similarity and the scale invariance of the fracture network are introduced to obtain the expression of the porosity and permeability of the complex fracture network. The change in fluid flow capacity is characterized by permeability. With the exception that part of the fluid flows from the matrix system to the fracture system, the fluid in region III flows into region II. Therefore, the governing equation of the fluid in the fracture system is given by: Combining the boundary conditions and initial conditions, the dimensionless seepage mathematical model for region II can be obtained after dimensionless transformation: In the above mathematical model, ω 2 is the elastic storage capacity ratio, and its expression is as follows: λ 2 is the dimensionless interporosity flow coefficient, and its expression is given by: (3) Mathematical Model for Region I Region I is the area of the primary fracture. The fluid in this area flows into the horizontal wellbore along the fracture, which is a one-dimensional linear flow. Considering the influence of stress sensitivity, the permeability of the main fracture changes with the effective confining pressure, and has an exponential function relationship with the overlying pressure. We can acquire the expression of the permeability of the main fracture as: Energies 2021, 14, 5857 The tip of the fracture is closed, and the fluid in region II flows into region I. Without considering the skin effect of the primary fracture and the wellbore storage effect, the governing equation of fluid in this area can be expressed as: Incorporating the boundary and initial conditions, the dimensionless seepage mathematical model for region I can be obtained after dimensionless treatment:

Solution of the Mathematical Model
The mathematical equations for region I and region II have strong nonlinear characteristics due to consideration of the influence of the stress sensitivity of the fracture system. The transformation method proposed by Pedrosa [26] can eliminate the nonlinearity of the equation, and its transformation formula is given by: The value of γ D is small in general. According to the perturbation theory, the following equations are defined as follows: Considering that γ D is very small, zero-order perturbation can provide the required accuracy of the solution. The above method is employed to solve the seepage mathematical model for area I and area II.

Solution of the Seepage Mathematical Model for Region III
Laplace transformation is applied to the dimensionless seepage mathematical model for region III to obtain the corresponding mathematical model in Laplace space: Solving the mathematical model, we can get: where F 32 = r 31 r 32 e r 31 −e (r 31 −r 32 )x eD e r 32 r 32 e r 31 −r 31 e (r 31 −r 32 )x eD e r 32 ,

Solution of the Seepage Mathematical Model for Region II
By performing zero-order perturbation transformation and Laplace transformation on the dimensionless seepage mathematical model for region II, the mathematical model in Laplace space is obtained as follows: Substituting Equation (33) into Equation (34) and simplifying the governing equations, we can acquire: Combined with the boundary conditions, the model is solved according to the Bessel function, and the result is obtained as follows: In above equation,

Solution of the Seepage Mathematical Model for Region I
Laplace transformation is applied to the dimensionless seepage mathematical model for region I to obtain the corresponding mathematical model in Laplace space: In above equation, The productivity at the inner boundary of the primary fracture is a quarter of the production of the bottom hole, and its relational expression is as follows: Substituting Equation (38) into the above equation, the production of a single fracture of a horizontal well under Laplace space is obtained as follows:

Solution of the Production of a Single Primary Fracture for Horizontal Wells
We obtained the production expression of a single fracture of the horizontal well in Laplace space in Section 3.1.3. The production expression in real space can be obtained by the Stehfest [27] numerical inversion according to the productivity formula in Laplace space. Combining with Equation (40), the corresponding dimensionless production expression of a single fracture at the bottom hole for any dimensionless time point t D is represented as follows: In the above expression, the Laplace variable s = ln 2 t D i.
Therefore, the dimensionless production expression at the bottom hole of a single fracture in real space is given by:

Solution of the Production of Staged Multi-Cluster Fractured Horizontal Wells with Non-Uniformly Distributed Fractures
Taking into account the current situation of tight oil exploitation with bridge plug staged multi-cluster perforation fracturing [28], a productivity prediction model for staged multi-cluster fractured horizontal wells was established in this study, and the fracture distribution diagram is represented in the following  Taking into account the current situation of tight oil exploitation with bridge plug staged multi-cluster perforation fracturing [28], a productivity prediction model for staged multi-cluster fractured horizontal wells was established in this study, and the fracture distribution diagram is represented in the following Figure 3. As shown in Figure 3, fractures of the horizontal well are unevenly arranged. The number of fractures in each fracturing section, the distance between fracturing sections, and the space between each fracture of fracturing sections are different. Taking the fracture at any position as the central axis, its control range is the sum of half the distance from the left and right adjacent fractures. (In particular, the range that can be controlled by fractures at both ends of the horizontal well is composed of two parts: the distance from the fractures at both ends of the horizontal well to the reservoir boundary and 1/2 of the distance between the fractures and their adjacent fractures.) We divide the fractures into three types based on their position at the horizontal well: the fractures at both ends of the horizontal well, the fractures at both ends of the fracturing segment, and the fractures inside the fracturing segment. The shape parameters of these are defined as follows.
The fractures at both ends of the horizontal well (∆ is the distance between the fractures at both ends of the horizontal well and the reservoir boundary): As shown in Figure 3, fractures of the horizontal well are unevenly arranged. The number of fractures in each fracturing section, the distance between fracturing sections, and the space between each fracture of fracturing sections are different. Taking the fracture at any position as the central axis, its control range is the sum of half the distance from the left and right adjacent fractures. (In particular, the range that can be controlled by fractures at both ends of the horizontal well is composed of two parts: the distance from the fractures at both ends of the horizontal well to the reservoir boundary and 1/2 of the distance between the fractures and their adjacent fractures.) We divide the fractures into three types based on their position at the horizontal well: the fractures at both ends of the horizontal well, the fractures at both ends of the fracturing segment, and the fractures inside the fracturing segment. The shape parameters of these are defined as follows.
The fractures at both ends of the horizontal well (∆y d is the distance between the fractures at both ends of the horizontal well and the reservoir boundary): The fractures at both ends of the fracturing segment (∆y ni is the distance between the last fracture of the ith fracturing segment and the first fracture of the (i + 1)th fracturing segment, i represents the sequence of the fracturing segment): The fractures inside the fracturing segment (∆y mij is the distance between the jth fracture inside the ith fracturing segment and the (j + 1)th fracture inside the same fracturing segment, j is the sequence of the fracture): According to the fracture production superposition calculation method proposed by Mayer [29], the single well production formula of staged multi-cluster fractured horizontal wells can be obtained as follows: In the above equation, ns expresses the number of fracturing segments and nc means the number of fractures in the fracturing segment. The total number of fracturing clusters in all fracturing sections is the total number of fractures.

Description of Symbols
p 0 is initial formation pressure, P a ; p w is bottom hole pressure, P a ; p is current fluid flow pressure, P a ; N F is the number of fracture; x f is the half-length of fracture, m; w f is the width of fracture, m; x e is the half-width of formation, m; y e is width of region II, m; B is oil volume factor; h is formation thickness, m; µ is oil viscosity, P a ·s; k is permeability, m 2 ; φ is porosity, %; G is threshold pressure gradient, P a /m; γ is stress sensitivity coefficient, P −1 a ; D is fractal dimension; θ is anomalous diffusion coefficient; α is shape factor of reservoir, m −2 ; C t is total compression coefficient, P −1 a ; C l is liquid phase compression coefficient, P −1 a ; q is the production of a single fracture of horizontal wells, m 3 /s; Q is the productivity of a single staged multi-cluster fractured horizontal well, m 3 /s; t is the production time of horizontal well, s; I, K represent the Bessel function; the subscripts 1, 2, and 3, respectively, represent the regions I, II, and III; the subscripts m, f respectively represent the matrix system and fracture system, respectively; the subscript 0 represents the initial value; subscript D means dimensionless; the superscript -means Laplace space; s is the Laplace variable.

Verification of the Analytical Solution
Without considering the non-uniform structure characteristics of the fracture network, the influence of threshold pressure gradient on fluid flow, and the stress sensitivity of the fracture system, our simplified model is similar to the classical trilinear flow model proposed by M. Brown et al. [17]. According to the expression of the dimensionless bottom hole pressure of the classical model, the dimensionless production at the bottom hole of a single fracture can be obtained as [30]: Setting G 3 = 0, γ = 0, D = 2, θ = 0 in the original mathematical model, the results of the simplified model were compared with the results of the classical model to verify the accuracy of the analytical solution of our original model. The simulation parameters are listed in Table 1. As shown in Figure 4, the results of the simplified model are highly consistent with the results of the classical model, which proves the accuracy of the analytical solution of our original model.  Setting = 0, = 0, = 2, = 0 in the original mathematical model, the results of the simplified model were compared with the results of the classical model to verify the accuracy of the analytical solution of our original model. The simulation parameters are listed in Table 1. As shown in Figure 4, the results of the simplified model are highly consistent with the results of the classical model, which proves the accuracy of the analytical solution of our original model.

An Example in an Oil Field
We can calculate the production of a single fractured horizontal well according to Equation (46). To verify the reliability of our model, the basic parameters of an oil well in the X Oilfield (in China) are substituted into the mathematical model for calculation, and the calculation results of the model fitted well with the production data of the well. The basic parameters and fracturing design parameters of the oil well, respectively, are shown in Tables 2 and 3.
The fitting results of the model and production data are shown in Figure 5. It can be seen from the figure that the calculation results of this model are basically consistent with the actual production data of the oil field within the allowable error range, which confirms the reliability of the model. The calculation results of the model deviate from the actual production data without considering the influence of stress sensitivity or the threshold

An Example in an Oil Field
We can calculate the production of a single fractured horizontal well according to Equation (46). To verify the reliability of our model, the basic parameters of an oil well in the X Oilfield (in China) are substituted into the mathematical model for calculation, and the calculation results of the model fitted well with the production data of the well. The basic parameters and fracturing design parameters of the oil well, respectively, are shown in Tables 2 and 3.
The fitting results of the model and production data are shown in Figure 5. It can be seen from the figure that the calculation results of this model are basically consistent with the actual production data of the oil field within the allowable error range, which confirms the reliability of the model. The calculation results of the model deviate from the actual production data without considering the influence of stress sensitivity or the threshold pressure gradient. In addition, the fitting effect is obviously inferior to the original model when simultaneously ignoring both stress sensitivity and the threshold pressure gradient, which further proves the accuracy of the original model.

Analysis of Factors Affecting Productivity of Fractured Horizontal Wells
To facilitate the analysis of the factors affecting productivity, it is assumed that the fractures are evenly distributed along the horizontal well, with a total of 36 fracturing segments, four fractures within each segment, 28 m between each fracturing segment, and 7 m between each fracture inside the fracturing segment. Equation (46) can be rewritten as follows: Based on Equation (48) and the reference data in Table 2, the factors influencing the productivity of horizontal wells were investigated according to the deliverability evaluation method proposed in this paper, mainly including fracturing parameters (fracture half-length, fracture number) and reservoir characteristic parameters (including matrix permeability, threshold pressure gradient, and stress sensitive coefficient).

The Effect of Fracture Number on Productivity
The number of fracturing sections was set as 30, 32, 34, and 36, with other parameters in Table 2 unchanged. The impact of the number of fractures on the productivity of the horizontal wells was analyzed, and the daily production and cumulative production of horizontal wells were obtained, as shown in Figure 6. As depicted in Figure 6, the production of horizontal wells rises with the increase in the number of fractures. Fractures can improve the fluid seepage environment. The number of fractures increases, leading to the increase in fluid seepage channels and a significant improvement in the seepage environment, which results in an increase in the production of the horizontal well. At the same time, the phenomenon is reflected in Figure 6 that the production growth of the horizontal wells slows with the increase in the number of fractures. There is mutual interference between hydraulic fractures, and the increase in the number of fractures leads to the stronger interference between the fractures, thus resulting in the decrease in the production growth of horizontal wells. In addition, the increase in the number of fractures leads to the growth in the fracturing construction costs. Therefore, the number of fractures should be controlled within a reasonable range during the process of fracturing design to ensure the positive growth in production benefits.

The Effect of Fracture Half-Length on Productivity
The value of was assumed to be 60, 90, 120, and 150 m, with other parameters in Table 2 unchanged. The impact of the half-length of fractures on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal wells were obtained, as shown in Figure 7. As depicted in Figure 6, the production of horizontal wells rises with the increase in the number of fractures. Fractures can improve the fluid seepage environment. The number of fractures increases, leading to the increase in fluid seepage channels and a significant improvement in the seepage environment, which results in an increase in the production of the horizontal well. At the same time, the phenomenon is reflected in Figure 6 that the production growth of the horizontal wells slows with the increase in the number of fractures. There is mutual interference between hydraulic fractures, and the increase in the number of fractures leads to the stronger interference between the fractures, thus resulting in the decrease in the production growth of horizontal wells. In addition, the increase in the number of fractures leads to the growth in the fracturing construction costs. Therefore, the number of fractures should be controlled within a reasonable range during the process of fracturing design to ensure the positive growth in production benefits.

The Effect of Fracture Half-Length on Productivity
The value of x f was assumed to be 60, 90, 120, and 150 m, with other parameters in Table 2 unchanged. The impact of the half-length of fractures on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal wells were obtained, as shown in Figure 7

The Effect of Fracture Half-Length on Productivity
The value of was assumed to be 60, 90, 120, and 150 m, with other parameters in Table 2 unchanged. The impact of the half-length of fractures on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal wells were obtained, as shown in Figure 7.  As depicted in Figure 7a, the productivity of horizontal wells rises with the increase in fracture length. At the initial stage during the exploitation of tight oil reservoirs, the half-length of fractures increases by 30 m, and the daily production of the horizontal well rises by about 10-50 m 3 . As presented in Figure 7b, the growth rate of final cumulative production decreases obviously with the increase in fracture length, and the fracturing cost expands with the increase in the hydraulic fracture length. Therefore, there is an optimal fracture half-length, which can not only improve the production of horizontal wells, but also save expenditure.

The Effect of Matrix Permeability on Productivity
The value of matrix permeability was set as 0.13, 0.15, 0.17, and 0.19 mD, with other parameters in Table 2 unchanged. The influence of matrix permeability on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal well were acquired, as depicted in Figure 8. As depicted in Figure 7a, the productivity of horizontal wells rises with the increase in fracture length. At the initial stage during the exploitation of tight oil reservoirs, the half-length of fractures increases by 30 m, and the daily production of the horizontal well rises by about 10-50 m³. As presented in Figure 7b, the growth rate of final cumulative production decreases obviously with the increase in fracture length, and the fracturing cost expands with the increase in the hydraulic fracture length. Therefore, there is an optimal fracture half-length, which can not only improve the production of horizontal wells, but also save expenditure.

The Effect of Matrix Permeability on Productivity
The value of matrix permeability was set as 0.13, 0.15, 0.17, and 0.19 mD, with other parameters in Table 2 unchanged. The influence of matrix permeability on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal well were acquired, as depicted in Figure 8. It can be seen from Figure 8a that, in the early and middle stages of exploitation, the daily productivity of the horizontal well rises with the increase in matrix permeability, and the higher the matrix permeability, the bigger the deceleration rate of the daily production of horizontal wells. As depicted in Figure 8b, the trend of the rising cumulative productivity slows with the increase in matrix permeability. After the horizontal well has produced for a certain time, the cumulative productivity remains largely unchanged with the increase in matrix permeability, which indicates that increasing matrix permeability cannot effectively improve the final cumulative productivity of the horizontal well, but can shorten the time to reach the same cumulative productivity. In general, matrix perme- It can be seen from Figure 8a that, in the early and middle stages of exploitation, the daily productivity of the horizontal well rises with the increase in matrix permeability, and the higher the matrix permeability, the bigger the deceleration rate of the daily production of horizontal wells. As depicted in Figure 8b, the trend of the rising cumulative productivity slows with the increase in matrix permeability. After the horizontal well has produced for a certain time, the cumulative productivity remains largely unchanged with the increase in matrix permeability, which indicates that increasing matrix permeability cannot effectively improve the final cumulative productivity of the horizontal well, but can shorten the time to reach the same cumulative productivity. In general, matrix permeability has a significant impact on the early stage and, in particular, the middle stage, of the whole producing process.

The Effect of the Stress Sensitivity Coefficient on Productivity
The value of the stress sensitivity coefficient was assumed to be 0.00, 0.01, 0.03, and 0.05 MP −1 a , with other parameters in Table 2 unchanged. The effect of the stress sensitivity coefficient on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of horizontal wells were obtained, as represented in Figure 9. As presented in Figure 9, the stress sensitivity coefficient has a significant impact on the productivity of horizontal wells. The production of the horizontal well decreases with the increase in the stress sensitivity coefficient, and the trend of the reduction in productivity slows. The fracture system is compressed because of the increase in the effective pressure in the process of development and exploration of the tight oil reservoir, thus resulting in the decrease in permeability, which further reduces the productivity of horizontal wells. As depicted in the figures, the productivity of the horizontal well when considering the stress sensitivity of the fracture system is significantly lower than that without stress sensitivity, which shows that stress sensitivity is an important factor that cannot be ignored in the deliverability evaluation of fractured horizontal wells in tight reservoirs.

The Effect of Threshold Pressure Gradient on Productivity
was assumed to be 0.00, 0.02, 0.04, and 0.06 MP m ⁄ , and other parameters in Table   2 were unchanged. The influence of the threshold pressure gradient on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal well were obtained, as presented in Figure 10. The variation in the daily production of horizontal wells with time is demonstrated as a double logarithm curve in Figure 10a. As presented in Figure 9, the stress sensitivity coefficient has a significant impact on the productivity of horizontal wells. The production of the horizontal well decreases with the increase in the stress sensitivity coefficient, and the trend of the reduction in productivity slows. The fracture system is compressed because of the increase in the effective pressure in the process of development and exploration of the tight oil reservoir, thus resulting in the decrease in permeability, which further reduces the productivity of horizontal wells. As depicted in the figures, the productivity of the horizontal well when considering the stress sensitivity of the fracture system is significantly lower than that without stress sensitivity, which shows that stress sensitivity is an important factor that cannot be ignored in the deliverability evaluation of fractured horizontal wells in tight reservoirs.

The Effect of Threshold Pressure Gradient on Productivity
G 3 was assumed to be 0.00, 0.02, 0.04, and 0.06 MP a /m, and other parameters in Table 2 were unchanged. The influence of the threshold pressure gradient on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal well were obtained, as presented in Figure 10. The variation in the daily production of horizontal wells with time is demonstrated as a double logarithm curve in Figure 10a. was assumed to be 0.00, 0.02, 0.04, and 0.06 MP m ⁄ , and other parameters in Table   2 were unchanged. The influence of the threshold pressure gradient on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of the horizontal well were obtained, as presented in Figure 10. The variation in the daily production of horizontal wells with time is demonstrated as a double logarithm curve in Figure 10a.
(a) (b) Figure 10. The effect of the threshold pressure gradient on productivity: (a) effect of the threshold pressure gradient on well production rate; (b) effect of the threshold pressure gradient on cumulative production rate. Figure 10. The effect of the threshold pressure gradient on productivity: (a) effect of the threshold pressure gradient on well production rate; (b) effect of the threshold pressure gradient on cumulative production rate.
It can be seen from Figure 10 that the threshold pressure gradient has an obvious influence on the productivity of horizontal wells. The greater the threshold pressure gradient, the bigger the reduction rate of the daily productivity of horizontal wells and the smaller the final cumulative productivity. This is because a larger threshold pressure gradient leads to stronger seepage resistance, and difficulty in obtaining the fluid flow results for the reduction in productivity. The threshold pressure gradient was used to characterize the non-Darcy seepage of fluid in this study. As shown in Figure 10, the threshold pressure gradient is a non-negligible factor leading to the reduction in productivity of horizontal wells. Therefore, the fluid seepage process cannot be simply viewed as a flow following Darcy's law in the process of the deliverability evaluation for staged multi-cluster fractured horizontal wells in tight reservoirs.  Table 2 were unchanged. The impact of the fractal dimension on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of horizontal wells were obtained, as shown in Figure 11. The variation in the daily production of horizontal wells with time is demonstrated as a double logarithm curve in Figure 11a. It can be seen from Figure 10 that the threshold pressure gradient has an obvious influence on the productivity of horizontal wells. The greater the threshold pressure gradient, the bigger the reduction rate of the daily productivity of horizontal wells and the smaller the final cumulative productivity. This is because a larger threshold pressure gradient leads to stronger seepage resistance, and difficulty in obtaining the fluid flow results for the reduction in productivity. The threshold pressure gradient was used to characterize the non-Darcy seepage of fluid in this study. As shown in Figure 10, the threshold pressure gradient is a non-negligible factor leading to the reduction in productivity of horizontal wells. Therefore, the fluid seepage process cannot be simply viewed as a flow following Darcy's law in the process of the deliverability evaluation for staged multi-cluster fractured horizontal wells in tight reservoirs.  Table 2 were unchanged. The impact of the fractal dimension on the productivity of horizontal wells was analyzed, and the daily production and cumulative production of horizontal wells were obtained, as shown in Figure 11. The variation in the daily production of horizontal wells with time is demonstrated as a double logarithm curve in Figure 11a.
(a) (b) Figure 11. The effect of the fractal dimension on productivity: (a) effect of the fractal dimension on well production rate; (b) effect of the fractal dimension on cumulative production rate.
In this study, the complex fracture network considering the fractal characteristics is called the non-uniform fracture network. As presented in Figure 11, the fractal dimension has a significant influence during the whole process of exploration of horizontal wells. The productivity of horizontal wells rises with the increase in the fractal dimension of the non-uniform fracture network, and the trend of rising productivity is also more obvious. In this study, the complex fracture network considering the fractal characteristics is called the non-uniform fracture network. As presented in Figure 11, the fractal dimension has a significant influence during the whole process of exploration of horizontal wells. The productivity of horizontal wells rises with the increase in the fractal dimension of the non-uniform fracture network, and the trend of rising productivity is also more obvious.

Conclusions
This paper proposes a new method of deliverability evaluation for staged fractured horizontal wells. The solution of the productivity model is obtained using the semianalytical method, and the reliability of the model was verified with a classic model and an oilfield example. According to the research into the deliverability evaluation method and factors affecting productivity, several conclusions of this paper are represented as follows.
(1) The non-uniform distribution fracture network structure along the horizontal wellbore is proposed in this paper. Combining fractal theory with stress sensitivity to characterize the permeability and porosity, a new three-zone composite unsteady deliverability evaluation model for staged multi-cluster fractured horizontal wells is established. (2) The number of fractures, fracture half-length, and matrix permeability can affect the productivity of horizontal wells to a certain extent. The oil production rises with the increase in fracture number, fracture half-length, and matrix permeability, but the growth rate of productivity decreases at the same time. (3) The stress sensitivity coefficient, threshold pressure gradient, and fractal dimension have obvious effects on the productivity of horizontal wells. The smaller the stress sensitivity coefficient and threshold pressure gradient, the greater the oil production. A larger fractal dimension leads to greater productivity of horizontal wells. (4) The deliverability evaluation method for staged fractured horizontal wells in tight reservoirs based on the non-uniform fracture network structure can effectively predict the actual production of oil wells. Additionally, it can provide a theoretical basis for the optimization of fracturing parameters of horizontal wells and the formulation of technical policies for rational development of tight reservoirs, which has great significance for realizing the effective development of tight reservoirs.
Although some innovations were achieved in this study, the productivity model does not consider the imbibition phenomenon in the development of tight oil reservoirs and the pressure loss along the horizontal wellbore. Therefore, a more comprehensive deliverability evaluation model will be established by considering both the imbibition effect and the pressure drop effect in our future research work.  Data Availability Statement: Restrictions apply to the availability of these data. Data was partly obtained from PetroChina Changqing Oilfield Company and partly generated during the process of calculations. Oilfield data are available from the corresponding author with the permission of PetroChina Changqing Oilfield Company and calculation data can be acquired from the lead author.