Productivity-Index Behavior for a Horizontal Well Intercepted by Multiple Finite-Conductivity Fractures Considering Nonlinear Flow Mechanisms under Steady-State Condition

In this paper, a mathematical model is proposed to investigate the effect of nonlinear flow mechanisms on productivity-index (PI) behavior in hydraulically fractured reservoirs during steady-state condition. This approach focuses on the fact that PI approaches a constant value at a certain time, indicating the beginning of steady state. In this model, the reservoirs are considered as an elliptical-shaped drainage with constant-pressure boundary, which is depleted by a multiple-fractured horizontal well (MFHW), and various nonlinear flow mechanisms, such as the non-Darcy flow effect and pressure-dependency effect, control flow patterns in the hydraulic fractures. Then, an exact algorithm of solving the resulting nonlinear equations is developed to obtain the PI of MFHW using a semi-analytical approach. Next, type curves are generated to investigate the influences of flow mechanisms and fracture properties on PI. The most interesting points in this study are the following: (1) PI is determined by the properties of MHFW (i.e., dimensions and configuration), the reservoir geometry, and flow mechanism; (2) PI is deteriorated by non-Darcy flow caused by inertial forces; and (3) PI is reduced under the influence of pressure sensitivity caused by the degradation of dynamic conductivity. Generally, this study provides a significant insight into understanding the factors affecting the productivity of a MFHW with nonlinear flow mechanisms.


Introduction
With the success of Bakken tight oil resource in the United States, tight oil and gas reservoirs have become a significant source of hydrocarbon supply globally [1]. Subsequently, the technology of horizontal drilling combined with hydraulically fracturing has effectively enabled commercial production from the reservoir with extremely low permeability [2,3]. In practice, the characteristics for a stimulated well are identified by flow regimes; for example, the approach of pressure-derivative is used to identify flow-regimes for pressure transient analysis. For a multi-fractured horizontal well, the sequence of flow-regimes may be more complex (Figure 1) [4]. At some point after compound linear flow caused by the fracture interference, a compound elliptical flow regime (regime 5 in Figure 1a) would appear as the transient pressure waves move past the ends of the fractures [5,6]. Put another way, the compound elliptical response is the pseudo-steady state in an infinite drainage area. It is  [4]. Lines represent isopotential lines; arrows represent streamlines. Flow-regime 1 corresponds to linear flow; flow-regime 2 corresponds to elliptical flow; flow-regime 3 corresponds to fracture interference; flow-regime 4 corresponds to compound linear, and flow-regime 5 corresponds to compound elliptical flow. (b) simulation of flow regimes for a stimulated MFHW in a tight sand gas reservoir. (simulated through numerical simulator).
In physics, the performance of a fractured reservoir is determined by numerous factors such as formation petrophysical properties, fluid properties, and reservoir and wellbore configuration [7]. Numerous methods are provided to capture the physics of fluid flowing in the fractured reservoirs. Originating from the pioneer work suggested by Warren and Root [8], the fractures are represented by a 3D uniformly spaced fracture network, and the matrix blocks are contributing sources. It is a fictitious homogeneous porous medium. As a result, multiple-continuum models are used to compute the transient response of uniform naturally (continuously) fractured reservoirs. In analytical modeling, these models consist of a set of one-dimensional linear flow regions, and the main difference is the way in which these regions are coupled [9][10][11]. Therefore, it is difficult for the multiple-continuum model to capture the intrinsic behavior of complexity of nonuniform naturally fractured reservoirs, such as hydraulically fractured reservoirs. To capture the strong heterogeneity between fracture and matrix, various improved techniques are presented to subdivise the grids, such as the embedded discrete fracture model (EDFM), Galerkin finite-element method (FEM), and discrete-element method (DEM) [12][13][14]. However, it is still a huge challenge to obtain accurate transmissibility by computing the irregular connections between fracture segments. Considering the high computational burden and the great difficulties in gridding, the mesh-free semi-analytical method is necessarily developed as an alternative to the numerical simulation method. On the basis of the works suggested by Cinco-Ley et al. [15], various semi-analytical models, which have the flexibility of accounting for the individual fracture properties, are proposed to investigate the influence factors on the performance of the fractured reservoirs [16][17][18][19]. It is proven that the semianalytical methods have huge potential in analyzing transient performance response and evaluating productivity of wells.
From the viewpoint of petroleum engineering, the effects of these parameters on reservoir performance can be represented by the productivity index (PI). It is defined as the ratio of production rate to a certain pressure drop. PI is determined by the flow regime, with different behavior characteristics with time. It starts out with a high value at an early time when the transient response is dominant. PI declines with time to a small value, and approximates to a constant value, known as stabilized PI, at a late time period when the compound elliptical response is reached [20]. Currently, most studies have put the focus on the transient productivity index to identify the flow regimes [21][22][23]. However, they found that the designed fractured wells do not usually result in the expected performance, which is attributed to the existence of nonlinear flows. From the viewpoint of multi-   [4]. Lines represent isopotential lines; arrows represent streamlines. Flow-regime 1 corresponds to linear flow; flow-regime 2 corresponds to elliptical flow; flow-regime 3 corresponds to fracture interference; flow-regime 4 corresponds to compound linear, and flow-regime 5 corresponds to compound elliptical flow. (b) simulation of flow regimes for a stimulated MFHW in a tight sand gas reservoir. (simulated through numerical simulator).
In physics, the performance of a fractured reservoir is determined by numerous factors such as formation petrophysical properties, fluid properties, and reservoir and wellbore configuration [7]. Numerous methods are provided to capture the physics of fluid flowing in the fractured reservoirs. Originating from the pioneer work suggested by Warren and Root [8], the fractures are represented by a 3D uniformly spaced fracture network, and the matrix blocks are contributing sources. It is a fictitious homogeneous porous medium. As a result, multiple-continuum models are used to compute the transient response of uniform naturally (continuously) fractured reservoirs. In analytical modeling, these models consist of a set of one-dimensional linear flow regions, and the main difference is the way in which these regions are coupled [9][10][11]. Therefore, it is difficult for the multiple-continuum model to capture the intrinsic behavior of complexity of nonuniform naturally fractured reservoirs, such as hydraulically fractured reservoirs. To capture the strong heterogeneity between fracture and matrix, various improved techniques are presented to subdivise the grids, such as the embedded discrete fracture model (EDFM), Galerkin finite-element method (FEM), and discrete-element method (DEM) [12][13][14]. However, it is still a huge challenge to obtain accurate transmissibility by computing the irregular connections between fracture segments. Considering the high computational burden and the great difficulties in gridding, the mesh-free semi-analytical method is necessarily developed as an alternative to the numerical simulation method. On the basis of the works suggested by Cinco-Ley et al. [15], various semi-analytical models, which have the flexibility of accounting for the individual fracture properties, are proposed to investigate the influence factors on the performance of the fractured reservoirs [16][17][18][19]. It is proven that the semi-analytical methods have huge potential in analyzing transient performance response and evaluating productivity of wells.
From the viewpoint of petroleum engineering, the effects of these parameters on reservoir performance can be represented by the productivity index (PI). It is defined as the ratio of production rate to a certain pressure drop. PI is determined by the flow regime, with different behavior characteristics with time. It starts out with a high value at an early time when the transient response is dominant. PI declines with time to a small value, and approximates to a constant value, known as stabilized PI, at a late time period when the compound elliptical response is reached [20]. Currently, most studies have put the focus on the transient productivity index to identify the flow regimes [21][22][23]. However, they found that the designed fractured wells do not usually result in the expected performance, which is attributed to the existence of nonlinear flows. From the viewpoint of multi-physics, the flowing Energies 2020, 13, 2015 3 of 22 in the underground porous media is the result of fluid-thermal-solid coupling [24]. The reservoir production is generally regarded as the isothermal process; therefore, the thermal is not considered in the production simulation. Strictly speaking, the rock deformation may have a significant effect on fluid flow owing to fluid-particle interactions. According to the conclusions presented by Zhang and Tahmasebi [25,26], the effect of solid deformation on fluid flow, especially in hydraulic fractures, could be incorporated by regarding the porosity and permeability as a function of stress changes. In other words, the conductivity of hydraulic fracture changes because of a variety of mechanisms, including fracture closure under high effective stress and permeability loss as a result of proppant embedment and crushing. Besides, non-Darcy flow may occur in the fracture when the well is produced at a high flow rate [27]. Therefore, the nonlinear flow effects have to be taken into account to investigate the effect of nonlinear behaviors on the performance of a vertical well and fractured well in the infinite-acting reservoir. Limited efforts have been made to examine the effects of non-Darcy flow and pressure-sensitive conductivity on the performance of horizontal wells with multiple fractures. Wu et al. [28] used a numerical approach to investigate the non-Darcy effect in an unconventional gas reservoir. Lin and Zhu [29] developed a slab source method to evaluate the performance of horizontal wells with or without fractures using an analytical approach. Besides the non-Darcy effect, pressure sensitivity in the fracture may also exert an important influence on the transient response.
A transient pressure analysis for a vertical well intercepting a dynamic-conductivity fracture has been conducted using analytical and numerical methods. Pedroso and Correa [30] developed a new model to investigate the effects of pressure-sensitive permeability on the flow behavior of a fractured well. Zhang et al. [31] presented the results of a study that analyzed the build-up of pressure in a vertically fractured well with a stress-sensitive conductivity in the fracture. Yao et al. [32] studies the effect of dynamic conductivity on the transient performance of multiple fractured horizontal well.
To the best of our knowledge, few papers are published to study the productivity index of a multiple-fractured horizontal well (MFHW) in a finite drainage area with an elliptical boundary under nonlinear flow conditions. The objective of this work is to quantify the effect of non-Darcy flow and pressure sensitivity on the inflow performance of a MFHW. First, we establish a hybrid model that describes Darcy flow in the reservoirs toward hydraulic fractures, as well as non-Darcy flow along the fracture with pressure-dependent conductivity under steady-state conditions. First, a mathematical model is established to describe the fluid flow depleted by a MFHW in a reservoir with elliptical-shaped boundary, with the flexibility of accounting for fracture properties (i.e., fracture dimensions and configuration). Next, nonlinear flow mechanisms in the fracture, including stress-sensitivity effect and non-Darcy flow condition, are taken into account, which results in the nonlinearity in partial differential equation. Then, the technique of dimension transformation is used to render the resulting nonlinear equations amenable to linear treatment, and an effective and efficient algorithm is developed to solve the model. Finally, a study for sensitivity analysis is introduced to investigate the effect of fracture properties and nonlinear flow mechanisms on the productivity-index behavior.

Physical Model
In this physical model, a horizontal well is located in the center of a homogeneous reservoir with an elongated elliptical boundary. The whole drainage area consists of several different flow patterns in four contiguous flow regions, as sketched in Figure 2, including (1) the hydraulic fractures, (2) the inner stimulated rock volume (SRV) region impacted by hydraulic fracturing, (3) the outer linear region beyond the extent of the SRV, and (4) the outer elliptical region beyond the tips of the horizontal wellbore. The wellbore is produced under constant-rate condition.  Flow in the reservoir is assumed to be single-phase fluid (i.e., pure gas) with slight compressibility and constant viscosity.  Fluid flowing in the fracture is assumed to be nonlinear.  Fluids from the reservoir enter the wellbore only through hydraulic fractures, and the pressure loss in the wellbore is neglected.

Variables Definitions
To introduce the mathematical model conveniently and concisely, the dimensionless variables are defined, respectively, by When the flow in the fracture satisfies non-Darcy effect owing to high velocity, the corresponding dimensionless variables are defined as When the fracture conductivity behaves pressure sensitivity, the corresponding dimensionless variables are defined as where pm is the pressure in the reservoir, Pa; pf is the pressure in the fracture, Pa; pi is the pressure on the elliptical boundary, Pa; pex is the pressure on the width of inner SRV region, Pa; pey is the pressure on the length of inner SRV region, Pa; (xofm, yofm) is the location of m-th fracture tip, m; xe is the width of inner SRV region, m; ye is the length of inner SRV region, m; xR is the width of the whole drainage region, m 3 ; h is the reservoir thickness, m; wf is the fracture width, m; Lf is the fracture length, m; rw is the radius of wellbore, m; Lref is the reference length, m; qfm is the function of inflow distribution along  Different from the anisotropic characteristics in a naturally fractured reservoir, where a large number of nature fractures are uniformly and globally distributed [33][34][35], this paper assumes that the reservoir is stimulated by several hydraulic fractures, so the resulting hydraulically fractured reservoir is regarded as isotropic. In addition, there are some necessary assumptions, as follows: • The reservoir is horizontal and of uniform thickness, with impermeable lower and upper boundaries. The pressure on the outer boundary in the x-y plane keeps constant.

•
Along the horizontal wellbore, multiple hydraulic fractures are evenly distributed. Hydraulic fractures are considered to be a finite conductivity.

•
The wellbore is produced under constant-rate condition.

•
Flow in the reservoir is assumed to be single-phase fluid (i.e., pure gas) with slight compressibility and constant viscosity.

•
Fluid flowing in the fracture is assumed to be nonlinear.

•
Fluids from the reservoir enter the wellbore only through hydraulic fractures, and the pressure loss in the wellbore is neglected.

Variables Definitions
To introduce the mathematical model conveniently and concisely, the dimensionless variables are defined, respectively, by When the flow in the fracture satisfies non-Darcy effect owing to high velocity, the corresponding dimensionless variables are defined as When the fracture conductivity behaves pressure sensitivity, the corresponding dimensionless variables are defined as where p m is the pressure in the reservoir, Pa; p f is the pressure in the fracture, Pa; p i is the pressure on the elliptical boundary, Pa; p ex is the pressure on the width of inner SRV region, Pa; p ey is the pressure on the length of inner SRV region, Pa; (x ofm , y ofm ) is the location of m-th fracture tip, m; x e is the width of inner SRV region, m; y e is the length of inner SRV region, m; x R is the width of the whole drainage region, m 3 ; h is the reservoir thickness, m; w f is the fracture width, m; L f is the fracture length, m; r w is the radius of wellbore, m; L ref is the reference length, m; q fm is the function of inflow distribution along the m-th fracture, m 2 /s; q wm is the production rate, m 3 /s; q c is the influx rate, m 3 /s; µ is the viscosity, Pa s; B is the volume factor; γ f is the permeability modulus, Pa −1 ; and β is the Beta factor, m −1 . Here, the subscript D represents dimensionless variables.

Fluid Flow in the Inner SRV Region
The SRV region is regarded as a rectangular reservoir with constant-pressure boundary containing a set of hydraulic fractures, as shown in Figure 3. The diffusion equation governing fluid flow is given as the following dimensionless form: subscript D represents dimensionless variables.

Fluid Flow in the Inner SRV Region
The SRV region is regarded as a rectangular reservoir with constant-pressure boundary containing a set of hydraulic fractures, as shown in Figure 3. The diffusion equation governing fluid flow is given as the following dimensionless form: The outer boundary conditions are satisfied by where there are Nf fractures in the SRV region, the dimensionless length of m-th fracture is LfDm, and the dimensionless coordinate of m-th fracture tip is (xofDm, yofDm). The pressure on the outer boundary parallel to the fracture is peDx, and the pressure on the outer boundary perpendicular to the fracture is peDy.
Using Fourier transformation and inverse transformation (Appendix A), the dimensionless pressure distribution caused by Nf fracture is given by

Fluid Flow in the Outer Region
The fluids stored in the outer region are depleted and then flow into the inner SRV region through the interface (denoted by the yellow block in Figure 4). The streamlines satisfy the hyperbolic shape (orange solid line in Figure 4), which are the counterpart of elliptical-shaped isopotential line (black dashed line). Once the fluids enter the inner region, the streamlines are in a linear shape around the interface. Therefore, the coordinate systems are divided into two subsets: 1D elliptical coordinate in the outer region (blue region), and 1D linear coordinate (yellow region).
The outer boundary satisfies the constant pressure of peD.  The outer boundary conditions are satisfied by where there are N f fractures in the SRV region, the dimensionless length of m-th fracture is L fDm , and the dimensionless coordinate of m-th fracture tip is (x ofDm , y ofDm ). The pressure on the outer boundary parallel to the fracture is p eDx , and the pressure on the outer boundary perpendicular to the fracture is p eDy . Using Fourier transformation and inverse transformation (Appendix A), the dimensionless pressure distribution caused by N f fracture is given by

Fluid Flow in the Outer Region
The fluids stored in the outer region are depleted and then flow into the inner SRV region through the interface (denoted by the yellow block in Figure 4). The streamlines satisfy the hyperbolic shape (orange solid line in Figure 4), which are the counterpart of elliptical-shaped isopotential line (black dashed line). Once the fluids enter the inner region, the streamlines are in a linear shape around the interface. Therefore, the coordinate systems are divided into two subsets: 1D elliptical coordinate in the outer region (blue region), and 1D linear coordinate (yellow region). In the outer elliptic region, as shown in Figure 4, there are two sets of coordinates: the Cartesian coordinates (x',y') and the elliptic coordinates (ξ,η). The relationship is given by We used the elliptic coordinates to simulate the flowing process. The outer and inner boundaries of the elliptic region, represented by ξ in the elliptic coordinates, could be expressed using spatial variables in the Cartesian coordinates: The corresponding pressure gradient is written as the following dimensionless expression, Note that the flux inflow in Equation (9) is continuous on the interface between the outer elliptical region and inner SRV region. Meanwhile, the flowing is modeled using linear coordinates. As a result, the continuous condition is written in the form of the average integral, which is given by It could be expressed as the function of peDx and peDy (Appendix B), which is written as follows: Likewise, the continuity condition of flux inflow between outer linear region and inner SRV region contributes to We can obtain the dimensionless pressure on the interface by combing Equations (11) and (12), y=y'+y e x=x'+x e /2 The outer boundary satisfies the constant pressure of p eD . In the outer elliptic region, as shown in Figure 4, there are two sets of coordinates: the Cartesian coordinates (x',y') and the elliptic coordinates (ξ,η). The relationship is given by We used the elliptic coordinates to simulate the flowing process. The outer and inner boundaries of the elliptic region, represented by ξ in the elliptic coordinates, could be expressed using spatial variables in the Cartesian coordinates: The corresponding pressure gradient is written as the following dimensionless expression, Note that the flux inflow in Equation (9) is continuous on the interface between the outer elliptical region and inner SRV region. Meanwhile, the flowing is modeled using linear coordinates. As a result, the continuous condition is written in the form of the average integral, which is given by Energies 2020, 13, 2015 It could be expressed as the function of p eDx and p eDy (Appendix B), which is written as follows: Likewise, the continuity condition of flux inflow between outer linear region and inner SRV region contributes to We can obtain the dimensionless pressure on the interface by combing Equations (11) and (12), where the relationship between influx distribution and production rate of m-th fracture is given as. Substituting Equation (13) into Equation (6) could yield the pressure distribution caused by N f fractures, which is the function with regard to influx distribution along the fracture.

Fluid Flow in the Fracture
For hydraulic fracture, the flow pattern is widely assumed to be 1D linear and incompressible [36]. As shown in Figure 5, the governing equation of fluid flow on the m-th fracture is described as the following dimensionless form: Energies 2017, 10, x FOR PEER REVIEW 7 of 22 where the relationship between influx distribution and production rate of m-th fracture is given as Substituting Equation (13) into Equation (6) could yield the pressure distribution caused by Nf fractures, which is the function with regard to influx distribution along the fracture.

Fluid Flow in the Fracture
For hydraulic fracture, the flow pattern is widely assumed to be 1D linear and incompressible [36]. As shown in Figure 5, the governing equation of fluid flow on the m-th fracture is described as the following dimensionless form: The corresponding boundary conditions are described as Integrating Equation (14) from 0 to xDm with respect to xDm based on Equation (15) twice, the result is given by where is the integral of Heaviside function.
As shown in Figure 5, the fluid flow is regarded as a 1D linear pattern in the region far from wellbore, and becomes a 1D radial pattern around the wellbore. It is different from the vertical fracture, where the vertical fracture is in lateral contact with the wellbore and only linear flow is considered. The additional pressure loss owing to the convergence of fluids into the horizontal  The corresponding boundary conditions are described as Integrating Equation (14) from 0 to x Dm with respect to x Dm based on Equation (15) twice, the result is given by Energies 2020, 13, 2015 As shown in Figure 5, the fluid flow is regarded as a 1D linear pattern in the region far from wellbore, and becomes a 1D radial pattern around the wellbore. It is different from the vertical fracture, where the vertical fracture is in lateral contact with the wellbore and only linear flow is considered. The additional pressure loss owing to the convergence of fluids into the horizontal wellbore should be considered. The effect of flow convergence is regarded as a skin factor, namely convergence flow skin: Equation (16) is rewritten as the final expression by

Model of a Conductivity Fracture with Non-Darcy Flow
Non-Darcy flow may occur in the fracture when the inertial forces may no longer be neglected compared with viscous forces. It is very common near production wells where local velocities can be very high. Extra pressure drop is required to overcome the inertial forces.
To account for the effect of extra pressure loss caused by high-velocity flow, which is proportional to the square of the velocity, the Forchheimer flow equation is presented [27]: The apparent fracture permeability in Equation (19) is defined as The corresponding apparent permeability is rewritten by Substituting Equation (21) into Equation (18) yields the following form: where the apparent conductivity is written as , and the flux on the cross section is given as . The corresponding integral form is given by Note that Equation (22) presents a general model for a uniform-conductivity under the non-Darcy-flow condition.

Model of a Conductivity Fracture with Pressure Sensitivity Effect
The variation of permeability caused by a stress change can be expressed as a function of pore pressure. For hydraulic fracture, a general model presented by Zhang et al. [31] is introduced to describe the relationship between fracture permeability and pore pressure next: According to the definitions of the dimensionless parameters, the pressure-sensitive conductivity of a fracture is Substituting Equation (25) into Equation (14) yields the following equation governing flow in the m-th fracture flowing: Note that Equation (26) presents a general model for a uniform-conductivity with the pressure sensitivity effect.

Dimension Transformation
Compared with the equation of uniform conductivity fracture described by Equations (14), (22) and (26) are subject to nonlinear equations. In the case of considering non-Darcy flow, the fracture conductivity is a function of flowing rate; in the case of considering pressure-sensitivity effect, the fracture conductivity is a function of pressure. In essence, both flow rate and pressure chang with the spatial variable. Here, we attempt to solve the problem in the same way as the linear treatment in Equation (14) by a dimension transformation, which is defined as Then, the fracture-flow equations (Equations (22) and (26)) can be written as Equation (28) has the same form as the fracture-flow equations for a uniform-conductivity fracture with constant conductivity given by Equation (18). The difference is that the variable ξ D is a function of the variable x D and depends on the distribution of conductivity C fD (x D ).

Discretization
In this study, the productivity index is achieved by discretizing the fracture panel into M segments with equal length ∆x Di , as shown in Figure 6a. Note that the equal-length fracture segment would be transformed into an unequal-length segment after using dimension transformation, as shown in Figure 6b. In the reservoirs, the dimensionless pressure drop on the j-th segment of the n-th fracture caused by the i-th segment of the m-th fracture is given by According to the continuity condition that the pressure must be continuous along the fracture surface, the following conditions must hold along the fracture plane: Moreover, the sum of all the segments' flow rate equals to the total flow rate at the m-th fracture: Combining Equations (29)-(32), we can obtain the fundamental coupled solutions by computing the apparent linear equations [36,37].

Computation Consideration
According to the previous statement, the dimensionless conductivity of CfDn is the function of pressure pfDn in the fracture with pressure-dependent conductivity, and the dimensionless conductivity of CfDn is the function of flowing rate qcDn in the fracture under non-Darcy flow. In a mathematical context, CfDn is a function with regard to the spatial variable. At a given k-th step, if the distribution of pressure pfDn or flowing rate qcDn was known, CfDn along the fracture would be obtained. In the reservoirs, the dimensionless pressure drop on the j-th segment of the n-th fracture caused by the i-th segment of the m-th fracture is given by where BD n,j m = IB n,j , and ∆p n,j In the fracture, the dimensionless pressure on the j-th segment of the n-th fracture is expressed as where I(ξ D ) = ξ D 0 dζ ζ 0 q f Dn (ς)dς. According to the continuity condition that the pressure must be continuous along the fracture surface, the following conditions must hold along the fracture plane: Moreover, the sum of all the segments' flow rate equals to the total flow rate at the m-th fracture: Combining Equations (29)-(32), we can obtain the fundamental coupled solutions by computing the apparent linear equations [36,37].

Computation Consideration
According to the previous statement, the dimensionless conductivity of C fDn is the function of pressure p fDn in the fracture with pressure-dependent conductivity, and the dimensionless conductivity of C fDn is the function of flowing rate q cDn in the fracture under non-Darcy flow. In a mathematical context, C fDn is a function with regard to the spatial variable. At a given k-th step, if the distribution of pressure p fDn or flowing rate q cDn was known, C fDn along the fracture would be obtained. Thus, the nonlinear governing equation for the (k + 1)-th step could be linearized on the assumption of known conductivity distribution on the k-th step.
In the case of considering the non-Darcy flow effect (Case 1), the formation for the k-th step yields In the case of considering the pressure-dependent conductivity effect (Case 2), the formation for the k-th step yields The flow distribution q fDn and pressure distribution p fDn along the fracture at the k-th time step would be achieved based on coupled solution. Next, the calculated p fDn was used to update the fracture conductivity (Case 1); the calculated q cDn was used to update the fracture conductivity (Case 2). The iterative procedure is repeated until the wellbore-pressure was converged. The calculation procedure is given as follows:

•
Step 1: Initial calculation, with k = 0, the fracture is assumed to be uniform (Case 1 and Case 2). By combing through Equation (29) to Equation (32), we can obtain (q fDn ) k . The fracture pressure (p fDn ) k (Case 1) and flowing rate (q fDn ) k (Case 2) would be then achieved from Equation (28).

Results and Sensitivity Analysis
The productivity index (PI) is defined as the ratio of the production rate to pressure drawdown, as follows: Put another way, the dimensionless PI is the reciprocal of flow resistivity. Here, A is the drainage area, γ is the Eular constant, and C A is the shape factor. Specially, r we is the effective wellbore radius determined by the geometry of drainage area, well configuration, and nonlinear flow mechanism.

Influence of Fracture Properties on PI
In this subsection, the nonlinear flow mechanisms are firstly not taken into account, and the focus is put on the impacts of the dimensions and configuration of MFHW on PI behavior. For the purpose of discussion, homogeneous fracture properties and MFHW configuration are assumed; the fractures are evenly spaced along horizontal wellbore, and all fractures are of identical properties. As a result, the PI is mainly determined by five factors, including dimensionless conductivity (C fD ), fracture number (N f ), penetration ratio of fracture length to inner SRV width (I x = L f /x e ), penetration ratio of fractured horizontal wellbore to inner SRV length (I y = D f /y e ) (D f is defined as the distance between two outmost fractures), and the drainage ratio of inner drainage to the whole drainage (I e = x e /x R ). Figure 7 shows the effect of dimensionless conductivity on PI in the case of different penetration ratios. We further assume that the fracture is symmetrical (i.e., the fracture length on the right hand side equals to the left hand side with regard to wellbore), and all fractures are located in the center of the drainage. The fracture conductivity changes from 0.1 to 1000 (i.e., C fD = 0.1, 0.5, 1, 5, 10, 50, 100, 500, and 1000). As shown in Figure 7a, at a given penetration ratio (I x ), PI is increased with the increase in dimensionless conductivity. However, after a certain dimensionless conductivity, the increase could be miniscule. In other words, beyond a certain dimensionless conductivity, the PI increase with conductivity essentially equals the PI decrease caused by the inner interference within the fracture. Generally, the threshold value is regarded as C fD,threshold = 300. When the C fD is larger than 300, there is no significant pressure drop within the fracture; the PI reaches the maximum, and does not increase with the conductivity. Figure 7b shows the PI derivative with regard to ln(C fD ) in the semi-log plot. For each given I x , the conductivity, corresponding to the maximum PI derivative, is defined as certain conductivity. Taking I x = 1, for example (J Dmax = 7.25), when the derivative is at a maximum, and that the certain value is C fD,certain = 7.5, the corresponding value of PI is J D = 4.54. This indicates that the ratio of J D /J Dmax equals 62.6. This means that PI of MFHW at C fD = 7.5 could reach 62.6% of the maximum, but PI only increases by 37.4% when the C fD is further increased from 7.5 to 300. Therefore, the optimum dimensionless conductivity is defined to as the certain value. The optimum conductivity indicates that the inflow from the reservoir could match the outflow of the fracture. As analyzed from Figure 7b, the larger the I x , the larger the optimum C fD .
Energies 2017, 10, x FOR PEER REVIEW 12 of 22 side equals to the left hand side with regard to wellbore), and all fractures are located in the center of the drainage. The fracture conductivity changes from 0.1 to 1000 (i.e., CfD = 0.1, 0.5, 1, 5, 10, 50, 100, 500, and 1000). As shown in Figure 7a, at a given penetration ratio (Ix), PI is increased with the increase in dimensionless conductivity. However, after a certain dimensionless conductivity, the increase could be miniscule. In other words, beyond a certain dimensionless conductivity, the PI increase with conductivity essentially equals the PI decrease caused by the inner interference within the fracture. Generally, the threshold value is regarded as CfD,threshold = 300. When the CfD is larger than 300, there is no significant pressure drop within the fracture; the PI reaches the maximum, and does not increase with the conductivity. Figure 7b shows the PI derivative with regard to ln(CfD) in the semi-log plot. For each given Ix, the conductivity, corresponding to the maximum PI derivative, is defined as certain conductivity. Taking Ix = 1, for example (JDmax = 7.25), when the derivative is at a maximum, and that the certain value is CfD,certain = 7.5, the corresponding value of PI is JD = 4.54. This indicates that the ratio of JD/JDmax equals 62.6. This means that PI of MFHW at CfD = 7.5 could reach 62.6% of the maximum, but PI only increases by 37.4% when the CfD is further increased from 7.5 to 300. Therefore, the optimum dimensionless conductivity is defined to as the certain value. The optimum conductivity indicates that the inflow from the reservoir could match the outflow of the fracture. As analyzed from Figure 7b, the larger the Ix, the larger the optimum CfD.  Figure 8 shows the effect of the number of fractures on PI behavior. As expected, the cases with a larger number exhibit a higher value of PI. That is, increasing the number contributes to an increased connected area of MHFW with the reservoirs, and an increase in the fracture-fracture interference. However, as the number of fractures increases beyond a certain number, the increase of PI could be miniscule. For example, PI increases by 8.97% when the number is increased from Nf = 2 to Nf = 3. In comparison, incremental PI would decline to 0.39% when the fracture number is further increased from Nf = 9 to Nf = 10. Overall, the effect of the increased connected area could offset the effect of the increase of fracture interference, which results from the increased fracture number. Besides, the optimum conductivity is decreased with the increase in fracture number, as shown in Figure 8b.  Figure 8 shows the effect of the number of fractures on PI behavior. As expected, the cases with a larger number exhibit a higher value of PI. That is, increasing the number contributes to an increased connected area of MHFW with the reservoirs, and an increase in the fracture-fracture interference. However, as the number of fractures increases beyond a certain number, the increase of PI could be miniscule. For example, PI increases by 8.97% when the number is increased from N f = 2 to N f = 3. In comparison, incremental PI would decline to 0.39% when the fracture number is further increased from N f = 9 to N f = 10. Overall, the effect of the increased connected area could offset the effect of the increase of fracture interference, which results from the increased fracture number. Besides, the optimum conductivity is decreased with the increase in fracture number, as shown in Figure 8b Figure 9 presents the effect of the ratio of inner SRV area to the whole area, i.e., Ie. Here the value that Ie=1 means that the area of inner SRV equals to the whole drainage area. Compared with the PI that Ie=1, the more approaching unity the value is, the larger the PI is, which is caused by the decrease of the proportion of inner SRV. Moreover, the tendency is more noticeable in the condition of high conductivity. When CfD = 0.1, the PI at Ie = 0.5 is that JD@Ie = 0.5 = 0.45, while the PI at Ie = 1 is that JD@Ie = 1 = 0.65. Hence, the ratio of Ie = 1 to Ie = 0.5 in PI value is 1.44. As a comparison, when CfD = 1000, the PI at Ie=0.5 is that JD@Ie = 0.5 = 0.72, while the PI at Ie = 1 is that JD@Ie=1 = 1.95. The ratio in PI is 2.7.

Influence of Non-Darcy Flow on PI
In this subsection, the non-Darcy flow mechanism is taken into account, and we investigate the effect of non-Darcy flow on PI behavior. Figure 10 shows the flow distribution along the fracture, where CfD = 1. In this figure, the conductivity is set to be relatively low. It is shown that the influx density, which is defined as the influx rate per length along fracture face, is concentrated around the wellbore and fracture tips. It is noted that the concentration region of the influx density is controlled by a larger pressure drop/depletion. The variable of (qDND)f, defined in Equation (2), is the Forchheimer number for a given inertial factor β. When the (qDND)f is increased, the influx density is increasingly more concentrated towards wellbore. Put another way, an extra pressure drop is  Figure 9 presents the effect of the ratio of inner SRV area to the whole area, i.e., I e . Here the value that I e = 1 means that the area of inner SRV equals to the whole drainage area. Compared with the PI that I e = 1, the more approaching unity the value is, the larger the PI is, which is caused by the decrease of the proportion of inner SRV. Moreover, the tendency is more noticeable in the condition of high conductivity. When C fD = 0.1, the PI at I e = 0.5 is that J D@Ie=0.5 = 0.45, while the PI at I e = 1 is that J D@Ie=1 = 0.65. Hence, the ratio of I e = 1 to I e = 0.5 in PI value is 1.44. As a comparison, when C fD = 1000, the PI at I e = 0.5 is that J D@Ie=0.5 = 0.72, while the PI at I e = 1 is that J D@Ie=1 = 1.95. The ratio in PI is 2.7.  Figure 9 presents the effect of the ratio of inner SRV area to the whole area, i.e., Ie. Here the value that Ie=1 means that the area of inner SRV equals to the whole drainage area. Compared with the PI that Ie=1, the more approaching unity the value is, the larger the PI is, which is caused by the decrease of the proportion of inner SRV. Moreover, the tendency is more noticeable in the condition of high conductivity. When CfD = 0.1, the PI at Ie = 0.5 is that JD@Ie = 0.5 = 0.45, while the PI at Ie = 1 is that JD@Ie = 1 = 0.65. Hence, the ratio of Ie = 1 to Ie = 0.5 in PI value is 1.44. As a comparison, when CfD = 1000, the PI at Ie=0.5 is that JD@Ie = 0.5 = 0.72, while the PI at Ie = 1 is that JD@Ie=1 = 1.95. The ratio in PI is 2.7.

Influence of Non-Darcy Flow on PI
In this subsection, the non-Darcy flow mechanism is taken into account, and we investigate the effect of non-Darcy flow on PI behavior. Figure 10 shows the flow distribution along the fracture, where CfD = 1. In this figure, the conductivity is set to be relatively low. It is shown that the influx density, which is defined as the influx rate per length along fracture face, is concentrated around the wellbore and fracture tips. It is noted that the concentration region of the influx density is controlled by a larger pressure drop/depletion. The variable of (qDND)f, defined in Equation (2), is the Forchheimer number for a given inertial factor β. When the (qDND)f is increased, the influx density is increasingly more concentrated towards wellbore. Put another way, an extra pressure drop is

Influence of Non-Darcy Flow on PI
In this subsection, the non-Darcy flow mechanism is taken into account, and we investigate the effect of non-Darcy flow on PI behavior. Figure 10 shows the flow distribution along the fracture, where C fD = 1. In this figure, the conductivity is set to be relatively low. It is shown that the influx density, which is defined as the influx rate per length along fracture face, is concentrated around the wellbore and fracture tips. It is noted that the concentration region of the influx density is controlled by a larger pressure drop/depletion. The variable of (q DND ) f , defined in Equation (2), is the Forchheimer number for a given inertial factor β. When the (q DND ) f is increased, the influx density is increasingly more concentrated towards wellbore. Put another way, an extra pressure drop is required to offset the effect of non-Darcy flow caused by the inertial force. Figure 11 presents the effect of fracture conductivity on influx-density distribution along the fracture. Without the non-Darcy effect (Forchheimer number  Figure 11a, the increase in conductivity makes the distribution of fluid influx more gentle. When C fD > 300 (i.e., infinite conductivity), which indicates that the pressure on any point in fracture panel is same as the pressure on the wellbore, the inflow flux is concentrated on the fracture tips. With the non-Darcy effect (Figure 11b), when C fD > 1, more volume of influx fluid is concentrated on the nearby region of wellbore. The phenomenon is consistent with the characteristics given in Figure 11a. As the conductivity increases, the effect of non-Darcy flow on the distribution of inflow flux becomes weak, until the distribution of inflow flux is independent of the non-Darcy flow effect.
Energies 2017, 10, x FOR PEER REVIEW 14 of 22 required to offset the effect of non-Darcy flow caused by the inertial force. Figure 11 presents the effect of fracture conductivity on influx-density distribution along the fracture. Without the non-Darcy effect (Forchheimer number equals zero) shown in Figure 11a, the increase in conductivity makes the distribution of fluid influx more gentle. When CfD > 300 (i.e., infinite conductivity), which indicates that the pressure on any point in fracture panel is same as the pressure on the wellbore, the inflow flux is concentrated on the fracture tips. With the non-Darcy effect (Figure 11b), when CfD > 1, more volume of influx fluid is concentrated on the nearby region of wellbore. The phenomenon is consistent with the characteristics given in Figure 11a. As the conductivity increases, the effect of non-Darcy flow on the distribution of inflow flux becomes weak, until the distribution of inflow flux is independent of the non-Darcy flow effect.   Figure 12 shows the effect of non-Darcy flow on PI behavior. As shown in Figure 12a, at any given conductivity, PI is the highest in the condition of (qDND)f = 0. PI increases with the increase of CfD, until it reaches the maximum value (JDmax) at CfD = 300. However, the increasing rate of PI is gradually slowed down with CfD; likewise, the maximum PI is independent of (qDND)f. In other words, non-Darcy flow plays a negative role in determining PI for the range of low-and intermediateconductivity (CfD < 300), but has no effect on the maximum PI, which corresponds to the infinite conductivity. The relationship between PI and the Forchheimer number is further shown in Figure  12b. It is shown that, although PI declined with the Forchheimer number, the decreasing rate slows down.  required to offset the effect of non-Darcy flow caused by the inertial force. Figure 11 presents the effect of fracture conductivity on influx-density distribution along the fracture. Without the non-Darcy effect (Forchheimer number equals zero) shown in Figure 11a, the increase in conductivity makes the distribution of fluid influx more gentle. When CfD > 300 (i.e., infinite conductivity), which indicates that the pressure on any point in fracture panel is same as the pressure on the wellbore, the inflow flux is concentrated on the fracture tips. With the non-Darcy effect (Figure 11b), when CfD > 1, more volume of influx fluid is concentrated on the nearby region of wellbore. The phenomenon is consistent with the characteristics given in Figure 11a. As the conductivity increases, the effect of non-Darcy flow on the distribution of inflow flux becomes weak, until the distribution of inflow flux is independent of the non-Darcy flow effect.   Figure 12 shows the effect of non-Darcy flow on PI behavior. As shown in Figure 12a, at any given conductivity, PI is the highest in the condition of (qDND)f = 0. PI increases with the increase of CfD, until it reaches the maximum value (JDmax) at CfD = 300. However, the increasing rate of PI is gradually slowed down with CfD; likewise, the maximum PI is independent of (qDND)f. In other words, non-Darcy flow plays a negative role in determining PI for the range of low-and intermediateconductivity (CfD < 300), but has no effect on the maximum PI, which corresponds to the infinite conductivity. The relationship between PI and the Forchheimer number is further shown in Figure  12b. It is shown that, although PI declined with the Forchheimer number, the decreasing rate slows down.   Figure 12 shows the effect of non-Darcy flow on PI behavior. As shown in Figure 12a, at any given conductivity, PI is the highest in the condition of (q DND ) f = 0. PI increases with the increase of C fD , until it reaches the maximum value (J Dmax ) at C fD = 300. However, the increasing rate of PI is gradually slowed down with C fD ; likewise, the maximum PI is independent of (q DND ) f . In other words, non-Darcy flow plays a negative role in determining PI for the range of low-and intermediate-conductivity (C fD < 300), but has no effect on the maximum PI, which corresponds to the infinite conductivity. The relationship between PI and the Forchheimer number is further shown in Figure 12b. It is shown that, although PI declined with the Forchheimer number, the decreasing rate slows down. Figure 13 shows the effect of configuration and dimension of fractures on PI behavior with the consideration of the non-Darcy flow effect. As shown in Figure 13a, PI loss in the case of a lower fracture number is relatively more noticeable than that in the case of a larger number. For example, when N f = 2 and (q DND ) f = 0, PI is that J D = 0.96; when N f = 2 and (q DND ) f = 10, PI is that J D = 0.75. The relative loss in PI is up to 22%, which is defined as 1 − J D@(qDND)f=10 /J D@(qDND)f=0 . In comparison, the relative loss in PI is only 9%. Figure 13b shows the penetration ratio of fracture length with respect to the inner SRV width (I x ). The effect of Forchheimer number is weaker when the I x is relatively small, but its effect would be amplified when the I x is relatively small. Figure 13c shows the penetration ratio of the inner SRV region with regard to whole drainage (I e ). In the small range (I e < 0.9), the relationship between I e and PI exhibits an approximate linear behavior. When I e > 0.9, PI is increased rapidly with the I e .  Figure 13 shows the effect of configuration and dimension of fractures on PI behavior with the consideration of the non-Darcy flow effect. As shown in Figure 13a, PI loss in the case of a lower fracture number is relatively more noticeable than that in the case of a larger number. For example, when Nf = 2 and (qDND)f = 0, PI is that JD = 0.96; when Nf = 2 and (qDND)f = 10, PI is that JD = 0.75. The relative loss in PI is up to 22%, which is defined as 1 -JD@(qDND)f=10/JD@(qDND)f=0. In comparison, the relative loss in PI is only 9%. Figure 13b shows the penetration ratio of fracture length with respect to the inner SRV width (Ix). The effect of Forchheimer number is weaker when the Ix is relatively small, but its effect would be amplified when the Ix is relatively small. Figure 13c shows the penetration ratio of the inner SRV region with regard to whole drainage (Ie). In the small range (Ie < 0.9), the relationship between Ie and PI exhibits an approximate linear behavior. When Ie > 0.9, PI is increased rapidly with the Ie.   Figure 13 shows the effect of configuration and dimension of fractures on PI behavior with the consideration of the non-Darcy flow effect. As shown in Figure 13a, PI loss in the case of a lower fracture number is relatively more noticeable than that in the case of a larger number. For example, when Nf = 2 and (qDND)f = 0, PI is that JD = 0.96; when Nf = 2 and (qDND)f = 10, PI is that JD = 0.75. The relative loss in PI is up to 22%, which is defined as 1 -JD@(qDND)f=10/JD@(qDND)f=0. In comparison, the relative loss in PI is only 9%. Figure 13b shows the penetration ratio of fracture length with respect to the inner SRV width (Ix). The effect of Forchheimer number is weaker when the Ix is relatively small, but its effect would be amplified when the Ix is relatively small. Figure 13c shows the penetration ratio of the inner SRV region with regard to whole drainage (Ie). In the small range (Ie < 0.9), the relationship between Ie and PI exhibits an approximate linear behavior. When Ie > 0.9, PI is increased rapidly with the Ie.

Influence of the Pressure-Sensitivity Effect on PI
The characteristic of the PI behavior is similar to the behavior of non-Darcy flow. First, Figure 14 shows the influx-distribution along the fracture with the consideration of the pressure-sensitivity effect. Compared with the nonsensitive case (Figure 14a), more influx density is distributed around the wellbore in the sensitive case (Figure 14b), and the tendency would be more remarkable with the decrease of the conductivity. The reason is explained in that the apparent conductivity is degraded and the flow resistance is consequently increased, which is caused by the closing of the fracture under the influence of rock compaction.
including (a) the influence of the number of fractures, (b) the influence of penetration ratio of fracture, and (c) the influence of penetration ratio of horizontal well.

Influence of the Pressure-Sensitivity Effect on PI
The characteristic of the PI behavior is similar to the behavior of non-Darcy flow. First, Figure  14 shows the influx-distribution along the fracture with the consideration of the pressure-sensitivity effect. Compared with the nonsensitive case (Figure 14a), more influx density is distributed around the wellbore in the sensitive case (Figure 14b), and the tendency would be more remarkable with the decrease of the conductivity. The reason is explained in that the apparent conductivity is degraded and the flow resistance is consequently increased, which is caused by the closing of the fracture under the influence of rock compaction.  Figure 15 shows the effect of pressure sensitivity on PI behavior. The intensity is measured by the dimensionless variable γfD. As shown in Figure 15a, the effect of pressure sensitivity on PI becomes weak with the increase of the conductivity. When the magnitude of the conductivity reaches the level of infinite conductivity, the effect of pressure sensitivity is almost neglected. As expected, as presented in Figure 15b, PI declines until the minimum with the increasing γfD, and corresponding γfD is defined as the threshold. Meanwhile, the decreasing rate slows down. For example, if CfD = 1, PI reaches the minimum JDmin when γfDthreshold = 1.2; if CfD = 10, PI reaches the minimum JDmin when γfDthreshold = 2.7. Thus, a larger conductivity contributes to a larger γfDthreshold.  Figure 15 shows the effect of pressure sensitivity on PI behavior. The intensity is measured by the dimensionless variable γ fD . As shown in Figure 15a, the effect of pressure sensitivity on PI becomes weak with the increase of the conductivity. When the magnitude of the conductivity reaches the level of infinite conductivity, the effect of pressure sensitivity is almost neglected. As expected, as presented in Figure 15b, PI declines until the minimum with the increasing γ fD , and corresponding γ fD is defined as the threshold. Meanwhile, the decreasing rate slows down. For example, if C fD = 1, PI reaches the minimum J Dmin when γ fDthreshold = 1.2; if C fD = 10, PI reaches the minimum J Dmin when γ fDthreshold = 2.7. Thus, a larger conductivity contributes to a larger γ fDthreshold .
including (a) the influence of the number of fractures, (b) the influence of penetration ratio of fracture, and (c) the influence of penetration ratio of horizontal well.

Influence of the Pressure-Sensitivity Effect on PI
The characteristic of the PI behavior is similar to the behavior of non-Darcy flow. First, Figure  14 shows the influx-distribution along the fracture with the consideration of the pressure-sensitivity effect. Compared with the nonsensitive case (Figure 14a), more influx density is distributed around the wellbore in the sensitive case (Figure 14b), and the tendency would be more remarkable with the decrease of the conductivity. The reason is explained in that the apparent conductivity is degraded and the flow resistance is consequently increased, which is caused by the closing of the fracture under the influence of rock compaction.  Figure 15 shows the effect of pressure sensitivity on PI behavior. The intensity is measured by the dimensionless variable γfD. As shown in Figure 15a, the effect of pressure sensitivity on PI becomes weak with the increase of the conductivity. When the magnitude of the conductivity reaches the level of infinite conductivity, the effect of pressure sensitivity is almost neglected. As expected, as presented in Figure 15b, PI declines until the minimum with the increasing γfD, and corresponding γfD is defined as the threshold. Meanwhile, the decreasing rate slows down. For example, if CfD = 1, PI reaches the minimum JDmin when γfDthreshold = 1.2; if CfD = 10, PI reaches the minimum JDmin when γfDthreshold = 2.7. Thus, a larger conductivity contributes to a larger γfDthreshold.  Figure 16 shows the effect of configuration and dimension of fractures on PI in the consideration of pressure sensitivity in the fracture. The overall features are similar to the case of considering non-Darcy flow, so we will not rewrite again.
Energies 2017, 10, x FOR PEER REVIEW 17 of 22 Figure 16 shows the effect of configuration and dimension of fractures on PI in the consideration of pressure sensitivity in the fracture. The overall features are similar to the case of considering non-Darcy flow, so we will not rewrite again.

Conclusions
We derived a new steady-state productivity model for a multi-fractured horizontal well in a reservoir with pseudo-elliptical boundary. We summarize some important findings as follows: 1. PI increases with the increase of CfD, until it reaches the maximum (JDmax) at CfD = 300. 2. PI is deteriorated under the influence of nonlinear flow mechanisms. With the consideration of non-Darcy flow, for the small range of the penetration ratio of the inner SRV region with regard to whole drainage (Ie < 0.9), the relationship between Ie and PI exhibits an approximately linear behavior. When Ie > 0.9, PI is increased rapidly with the Ie. 3. With the consideration of pressure sensitivity, the apparent fracture is degraded from the initial conductivity to the minimal conductivity, which is caused by fracture closure. As a result, an extra pressure drop is acquired to offset the conductivity degradation, and PI would be declined to the minimal PI. 4. The disadvantage dimensions of fracture (such as small conductivity, penetration ratio, and less fractures) contribute to severe pressure depletion, while, in turn, the severe pressure depletion will strengthen the effect of the nonlinear flow mechanism on PI behavior. 5. If the conductivity in the fracture reaches the level of infinite conductivity, the influence of nonlinear flow mechanism on PI could be neglected.

Conclusions
We derived a new steady-state productivity model for a multi-fractured horizontal well in a reservoir with pseudo-elliptical boundary. We summarize some important findings as follows: 1.
PI increases with the increase of C fD , until it reaches the maximum (J Dmax ) at C fD = 300.

2.
PI is deteriorated under the influence of nonlinear flow mechanisms. With the consideration of non-Darcy flow, for the small range of the penetration ratio of the inner SRV region with regard to whole drainage (I e < 0.9), the relationship between I e and PI exhibits an approximately linear behavior. When I e > 0.9, PI is increased rapidly with the I e .

3.
With the consideration of pressure sensitivity, the apparent fracture is degraded from the initial conductivity to the minimal conductivity, which is caused by fracture closure. As a result, an extra pressure drop is acquired to offset the conductivity degradation, and PI would be declined to the minimal PI. 4.
The disadvantage dimensions of fracture (such as small conductivity, penetration ratio, and less fractures) contribute to severe pressure depletion, while, in turn, the severe pressure depletion will strengthen the effect of the nonlinear flow mechanism on PI behavior.

5.
If the conductivity in the fracture reaches the level of infinite conductivity, the influence of nonlinear flow mechanism on PI could be neglected.

Funding:
This work was supported by the National Science and Technology Major Project (No. 2017ZX05019005-006) and Heilongjiang Natural Science Foundation (LH2019F004). Special thanks to the anonymous reviewers and editors for their valuable comments.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Analytical Solution for the Inner SRV Region
The appendix is an optional section that can contain details and data supplemental to the main text. For example, explanations of experimental details that would disrupt the flow of the main text, but nonetheless remain crucial to understanding and reproducing the research shown; figures of replicates for experiments of which representative data is shown in the main text can be added here if brief, or as Supplementary data. Mathematical proofs of results not central to the paper can be added as an appendix.

Appendix B. Analytical Solution for the Outer Region
On the interface between the inner SRV region and outer elliptical region, the continuity condition is given as In detail, the right hand side (RHS) of Equation (A18) is expressed as sinh(υπx eD /y eD ) cosh In detail, the right hand side (RHS) of Equation (A22) is expressed as