*4.2. Discretization*

In this study, the productivity index is achieved by discretizing the fracture panel into *M* segments with equal length Δ*xDi*, as shown in Figure 6a. Note that the equal-length fracture segmen<sup>t</sup> would be transformed into an unequal-length segmen<sup>t</sup> after using dimension transformation, as shown in Figure 6b.

**Figure 6.** Illustration of dimension transformation and fracture discretization.

In the reservoirs, the dimensionless pressure drop on the *j*-th segmen<sup>t</sup> of the *n*-th fracture caused by the *i*-th segmen<sup>t</sup> of the *m*-th fracture is given by

$$p\_{mDn,j} = \sum\_{m=1}^{N\_f} q\_{wDm} B D\_m^{n,j} + \sum\_{m=1}^{N\_f} \sum\_{i=1}^{M\_m} q\_{fDm,i} \Delta p\_{uDm,i}^{n,j} \frac{\Delta \xi\_{Dm,i}}{\Delta x\_{Dm,i}} \tag{29}$$

where *BD<sup>n</sup>*,*<sup>j</sup> m* = *IBn*,*<sup>j</sup> <sup>B</sup>*2*C*1,*m*<sup>−</sup>*B*1*C*2,*m A*1*B*2−*A*2*B*<sup>1</sup> + *IDn*,*<sup>j</sup> <sup>A</sup>*1*C*2,*m*<sup>−</sup>*A*2*C*1,*m A*1*B*2−*A*2*B*<sup>1</sup> , and <sup>Δ</sup>*p<sup>n</sup>*,*<sup>j</sup> uDm*,*i* = - *xo f Dm*,*i xo f Dm*,*i*−1 <sup>δ</sup>*puD*(*u*)*du*.

In the fracture, the dimensionless pressure on the *j*-th segmen<sup>t</sup> of the *n*-th fracture is expressed as

$$p\_{\rm wDn} - p\_{f\rm Dn,j} = \frac{2\pi}{\hat{\mathbb{C}}\_{f\rm Dn}} q\_{\rm wDn} G(\boldsymbol{\xi}\_{\rm Dn,j}, \boldsymbol{\xi}\_{\rm wDn}) - \left[ I(\boldsymbol{\xi}\_{\rm Dn,j}) - I(\boldsymbol{\xi}\_{\rm wD}) \right] - q\_{\rm wDn} \mathbf{S}\_{\rm c,n} \tag{30}$$

where *<sup>I</sup>*(ξ*D*) = - ξ*D* 0 *d*ζ- ζ0 .*qf 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:

$$p\_{mD}(\mathbf{x}\_{ofDn} + \mathbf{x}\_{Dn,j}, y\_{ofDn}) = p\_{fD}(\mathbf{x}\_{Dn,j}) \tag{31}$$

Moreover, the sum of all the segments' flow rate equals to the total flow rate at the *m*-th fracture:

$$q\_{wfDm} = \sum\_{i=1}^{M\_m} \left( \overline{q}\_{fDm,i} \frac{\Delta \xi\_{Dm,i}}{\Delta x\_{Dm,i}} \right) \tag{32}$$

Combining Equations (29)–(32), we can obtain the fundamental coupled solutions by computing the apparent linear equations [36,37].

## *4.3. 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. 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

$$\frac{\mathbb{C}\_{fDn}}{1 + (q\_{\rm NDD})\_{fn} q\_{\rm cDu}^{\langle k \rangle}} \frac{\partial p\_{fDn}^{\langle k+1 \rangle}}{\partial \mathbf{x}\_{\rm Dn}} + 2\pi q\_{\rm cDu}^{\langle k+1 \rangle} (\mathbf{x}\_{\rm Dn}) + 2\pi q\_{\rm wDn}^{\langle k+1 \rangle} H(\mathbf{x}\_{\rm Dn}, \mathbf{x}\_{\rm wDn}) = 0 \tag{33}$$

In the case of considering the pressure-dependent conductivity effect (Case 2), the formation for the *k*-th step yields

$$2\operatorname{\mathcal{L}}\_{f\operatorname{Dn}}(p\_{f\operatorname{Dn}'}^{\langle k\rangle}\mathbf{x}\_{\operatorname{Dn}})\frac{\partial p\_{f\operatorname{Dn}}^{\langle k+1\rangle}}{\partial \mathbf{x}\_{\operatorname{Dn}}} + 2\pi \Big[q\_{\operatorname{uDn}}^{\langle k+1\rangle}\mathbf{H}(\mathbf{x}\_{\operatorname{Dn},\mathbf{x}}\mathbf{x}\_{\operatorname{Dn}}) - \int\_{0}^{\operatorname{\mathcal{U}\mathcal{D}n}} q\_{f\operatorname{Dn}}^{\langle k+1\rangle}(\mathbf{x}')d\mathbf{x}'\Big] + 2\pi q\_{\operatorname{uDn}}^{\langle k+1\rangle}\mathbf{H}(\mathbf{x}\_{\operatorname{Dn}}\mathbf{x}\_{\operatorname{uDn}}) = 0 \tag{34}$$

The flow distribution *qfDn* and pressure distribution *pfDn* along the fracture at the *k*-th time step would be achieved based on coupled solution. Next, the calculated *pfDn* was used to update the fracture conductivity (Case 1); the calculated *qcDn* 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:


#### **5. Results and Sensitivity Analysis**

The productivity index (PI) is defined as the ratio of the production rate to pressure drawdown, as follows:

$$J\_D = \frac{\mu B}{2\pi k\_m h} \frac{q}{p\_i - p\_w} = \frac{1}{p\_{wD}} = \left(\frac{1}{2} \ln \frac{4A}{e^{\circ}C\_A r\_{w\epsilon}^2}\right)^{-1} \tag{35}$$

Put another way, the dimensionless PI is the reciprocal of flow resistivity. Here, *A* is the drainage area, γ is the Eular constant, and *CA* is the shape factor. Specially, *rwe* is the effective wellbore radius determined by the geometry of drainage area, well configuration, and nonlinear flow mechanism.

#### *5.1. 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 (*CfD*), fracture number (*Nf*), penetration ratio of fracture length to inner SRV width (*Ix* = *Lf*/*xe*), penetration ratio of fractured horizontal wellbore to inner SRV length (*Iy* = *Df*/*ye*) (*D*f is defined as the distance between two outmost fractures), and the drainage ratio of inner drainage to the whole drainage (*Ie* = *xe*/*xR*).

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., *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 *I*x, 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 7.** Productivity-index behavior of multiple fractures, including (**a**) the productivity index vs. dimensionless conductivity, and (**b**) the derivative of productivity index vs. dimensionless conductivity.

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.** Influence of number of fractures on PI behavior, including (**a**) the productivity index vs. dimensionless conductivity, and (**b**) the derivative of productivity index vs. dimensionless conductivity.

Figure 9 presents the effect of the ratio of inner SRV area to the whole area, i.e., *I*e. 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*=<sup>1</sup> = 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*=<sup>1</sup> = 1.95. The ratio in PI is 2.7.

**Figure 9.** Influence of penetration ratio of inner SRV region with regard to the whole region on productivity-index (PI) behavior.

#### *5.2. 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 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 e ffect (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 e ffect 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 e ffect.

**Figure 10.** Influx-flow distribution along the fracture under the influence of non-Darcy flow.

**Figure 11.** Influx-flow distribution along the fracture under the influence of conductivity (**a**) without and (**b**) with non-Darcy flow e ffect.

Figure 12 shows the e ffect 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 intermediate-conductivity (*CfD* < 300), but has no e ffect 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 e ffect of configuration and dimension of fractures on PI behavior with the consideration of the non-Darcy flow e ffect. 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 *J*D = 0.75. The relative loss in PI is up to 22%, which is defined as 1 − *JD@(qDND)f*=<sup>10</sup>/*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 e ffect of Forchheimer number is weaker when the *Ix* is relatively small, but its e ffect 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 *I*e and PI exhibits an approximate linear behavior. When *Ie* > 0.9, PI is increased rapidly with the *Ie*.

**Figure 12.** (**a**) Influence of fracture conductivity on PI behavior under the effect of non-Darcy flow, (**b**) influence of Forchheimer number on PI behavior with the consideration of finite conductivity.

**Figure 13.** Influence of dimensions of fracture on PI behavior in the consideration of non-Darcy flow, 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.

#### *5.3. Influence of the Pressure-Sensitivity E*ff*ect 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 14.** Influence of dimensionless conductivity on influx-flow distribution (**a**) with and (**b**) without the pressure-sensitivity effect.

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.** (**a**) Influence of fracture conductivity on PI behavior with the pressure-sensitivity effect, and (**b**) influence of permeability modulus on PI behavior with consideration of finite conductivity.

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.

**Figure 16.** Influence of dimensions of the fracture on PI behavior in the consideration of pressure sensitivity, including (**a**) number of fracture, (**b**) penetration index of fracture, and (**c**) penetration of horizontal wellbore.
