**4. Discussion**

The results of this study can be summarized as four types: (1) derivation of the approximate analytical solution; (2) validation of the solution against different numerical models; (3) introducing a step-by-step procedure to predict the values of physical parameters; (4) application of the analytical model in the field case.

Result 1 is one of the novelty of the article. The model contains three regions and effectively accounts for non-Darcy flow in the hydraulic fractures, and gas adsorption and slippage in the matrix. During the derivation process, the governing PDEs are transformed into ordinary differential equations (ODEs) by integration, instead of the Laplace transform. There is no doubt that Laplace transform works extremely well. However, there are still some problems: (1) the first step for Laplace transform is dimensionless transformation. For some dimensionless variables such as dimensionless time, dimensionless pressure, dimensionless production, and so on, more inputs are needed and many are even estimated, which will lead to calculation error. (2) The numerical inversion is an essential step for converting Laplace space into real-time space. Among them, Stefest numerical inversion algorithm is most commonly used. In this algorithm, the number of inversion items *N* is uncertain. The improper value will result in deviations in real time. Certainly, Result 1 also has some imperfections. For example, a dual-porosity model with Knudsen diffusion [6] would be more representative and the heterogeneity deserves further study.

Result 2 is the key section in the article. Seven numerical cases were set for verification. According to the fitting results, it shows that the analytical solution is feasible for the irregular and asymmetric stimulated regions in a multifractured horizontal well. Considering that the shapes of the enhanced fracture regions are unknown in real cases, we creatively set three types of stimulated regions: regular region, irregular region, and very irregular region. At least, the validation results show that our model is robust.

Result 3 is another novelty of our work. It introduced a step-by-step procedure to calculate inversion parameters. Comparing with the given volumes and calculated volumes from case 1 to case 7, the model is also verified to meet the engineering requirements. Considering the simplifying assumptions, the model may need to be improved and the accurate microseismic data would be required to make further verification.

Result 4 is the most important section. Considering the decline characteristics of shale gas, the constant rate case is not as important as the constant pressure case for long term performance of tight/shale formations. Therefore, our model is derived based on the constant bottom-hole pressure

condition, and namely one of limitations for our model is that it cannot be applied in a constant rate case. The second limitation is the data continuity, because the prerequisite for the model to make accurate prediction is grea<sup>t</sup> fitting. As for the discontinuous data like data missing, shut-in, or pressure/production jump, our model is not applicable. Combined with the assumption of single-phase flow in this work, the multiphase flow problem cannot be solved by this analytical solution. Gas–water two-phase flow is the most common in shale gas reservoirs, so our future work is to derive the new analytical solution for gas–water two-phase flow.
