**1. Introduction**

Unconventional hydrocarbon resources such as tight and shale oil/gas store in tight formations with ultra-low permeability. With the development of hydraulic fracturing technologies, multifractured horizontal wells (MFHWs) have rapidly emerged as the primary means for exploiting this type of resource. Meanwhile, some technologies are utilized, such as foam injection and carbon dioxide injection [1,2] on recovery enhancement, and photo-Fenton and floatation on sustainable managemen<sup>t</sup> of flow-back water after hydraulic fracturing [3]. In unconventional reservoirs, due to the propagation of fractures in different directions, branch fractures will be created around the main hydraulic fractures, which have a significant impact on the pressure and rate transient analysis for the fluid flow in the reservoirs.

In order to analyze production data and make long-term forecasts, analytical and numerical tools have been developed. Among them, a large number of numerical approaches [4–8], such as finite element method and boundary element method, are adopted to study the multiple flow mechanisms and fracture characteristics in shale gas reservoirs. Compared with the numerical models, the analytical models are simpler to set up and require fewer original data. Therefore, they are widely used in practice for field applications. Extensive studies have focused on analyzing production data of transient flow period for such complex systems [9–14]. Lee et al. [15] derived a tri-linear model to represent the transient flow behavior within a single fracture in an infinite homogeneous reservoir. Wattenbarger et al. [16] proposed an analytical solution accounting for transient linear flow contributions from each region for a fracture with infinite conductivity in a closed reservoir. El-Banbi [17] and Bello et al. [18] extended previous analytical models to develop a series of linear dual-porosity solutions for fluid flow in tight reservoirs, where constant-rate and constant-pressure inner boundaries were considered. Ozkan et al. [19] and Brown et al. [20] extended the concept of tri-linear flow and presented a solution to simulate the pressure-transient behavior in three flow regions, including the hydraulic fracture, the inner region between hydraulic fractures, and the outer region beyond the tip of hydraulic fractures. Stalgorova et al. [21] and Heidari et al. [22] introduced an enhanced-fracture-region (EFR) model to consider the cases where branch fractures cover the interfracture region. Mahdi et al. [23] firstly provided the transient shape factor among di fferent regions in the triple-porosity model.

The models represent an excellent attempt to capture reservoir parameters and flow characteristics for the complex system associated with the MFHWs. However, most analytical solutions were derived by Laplace transformation [17–23] and Boltzmann's transformation [24]. Numerical inversion algorithm [25] is needed to transform the Laplace space into real-time space. In order to avoid the inconvenience of Laplace transformation, several approximated analytical solutions have been presented. Shahamat et al. [26] introduced the concept of continuous succession of pseudo-steady states based on the capacitance-resistance methodology to analyze the production data in tight and shale reservoirs. In their work, the gas adsorption was ignored, which would underestimate the shale gas production. Ogunyomi [27,28] introduced a model with two compartments in which the first compartment represents the volume of the enhanced fracture region and the second one is the volume of the matrix. The approximate solution was obtained by converting governing partial-di fferential equations into a system of ordinary di fferential equations through integration. Then the system of ordinary di fferential equations was solved by calculating the eigenvalues and eigenvectors. Qiu et al. [29,30] extended this method to a triple-porosity model for production analysis. Nevertheless, they both did not take the non-Darcy flow into account and just derived the approximate solutions based on the oil phase rather than gas phase with the pressure dependent parameters. In this work, we extended the method to consider the gas flow and even account for non-Darcy flow in the hydraulic fractures, and gas adsorption and slippage in the matrix.

In this work, our goal is to derive a new integrative analytical solution for shale gas reservoirs bypassing the Laplace transformation, which can e ffectively account for non-Darcy flow in the hydraulic fractures, and gas adsorption and slippage in the matrix. The analytical model is created on the basis of a three-region model, which consists of the hydraulic fracture region, the inner region also called stimulated region connected to the hydraulic fracture, and the outer region representing the low permeability matrix. Then we introduce modified pseudo-variables [31,32], i.e., pseudo-pressure and pseudo-time, to eliminate the nonlinearity of governing equations. The whole derivation is based on the integration in real-time space bypassing Laplace transform and numerical inversion. The derived analytical solution is validated with numerical simulations. Further, several numerical cases are designed to verify the applicability in an irregular stimulated region using the new model. The proposed method is finally applied for shale gas production analysis and forecasting in the field.
