**1. 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 noted that the isopotential line approximates to an elongated elliptical shape for a long horizontal wellbore, as shown in Figure 1b.

**Figure 1.** The flow regime characteristics of MHFW, where (**a**) sequence of flow-regimes for a MFHW, completed in a tight sand gas reservoir modified from Clarkson [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 di fference is the way in which these regions are coupled [9–11]. Therefore, it is di fficult 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–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 grea<sup>t</sup> di fficulties 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–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 e ffects 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 di fferent 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–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 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 e ffect on

fluid flow owing to fluid-particle interactions. According to the conclusions presented by Zhang and Tahmasebi [25,26], the e ffect 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 e ffective 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 e ffects have to be taken into account to investigate the e ffect of nonlinear behaviors on the performance of a vertical well and fractured well in the infinite-acting reservoir. Limited e fforts have been made to examine the e ffects 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 e ffect 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 e ffect, 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 e ffects 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 e ffect 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 e ffect 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 e ffect and non-Darcy flow condition, are taken into account, which results in the nonlinearity in partial di fferential equation. Then, the technique of dimension transformation is used to render the resulting nonlinear equations amenable to linear treatment, and an e ffective and e fficient algorithm is developed to solve the model. Finally, a study for sensitivity analysis is introduced to investigate the e ffect of fracture properties and nonlinear flow mechanisms on the productivity-index behavior.

## **2. Model Assumption**
