An Analytical Model for Production Analysis of Hydraulically Fractured Shale Gas Reservoirs Considering Irregular Stimulated Regions

: Shale gas reservoirs are typically developed by multistage, propped hydraulic fractures. The induced fractures have a complex geometry and can be represented by a high permeability region near each fracture, also called stimulated region. In this paper, a new integrative analytical solution coupled with gas adsorption, non-Darcy ﬂow e ﬀ ect is derived for shale gas reservoirs. The modiﬁed pseudo-pressure and pseudo-time are deﬁned to linearize the nonlinear partial di ﬀ erential equations (PDEs) and thus the governing PDEs are transformed into ordinary di ﬀ erential equations (ODEs) by integration, instead of the Laplace transform. The rate vs. pseudo-time solution in real-time space can be obtained, instead of using the numerical inversion for Laplace transform. The analytical model is validated by comparison with the numerical model. According to the ﬁtting results, the calculation accuracy of analytic solution is almost 99%. Besides the computational convenience, another advantage of the model is that it has been validated to be feasible to estimate the pore volume of hydraulic region, stimulated region, and matrix region, and even the shape of regions is irregular and asymmetrical for multifractured horizontal wells. The relative error between calculated volume and given volume is less than 10%, which meets the engineering requirements. The model is ﬁnally applied to ﬁeld production data for history matching and forecasting.


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 management 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][5][6][7][8], such as finite element method and boundary element method, are adopted to study the multiple flow to model such a complex system, the reservoir is simplified as a three-region system, where the first region is the hydraulic fracture region which is regarded as the sole connection to the well. The second one represents the stimulated region with aggregated volume of all the microfractures and the third one is the adjacent ultra-low-permeability matrix directly connected to the stimulated region. Figure 1a shows the schematic 3D model. Three regions are contained in the model: region 1 is the hydraulic fracture, region 2 (darker color) with higher permeability around each hydraulic fracture, and region 3 (lighter color) with lower permeability connected to the region 2. The arrows represent the flow directions. For this model, the flow directions are parallel, and the system is symmetrical with respect to the hydraulic fracture and horizontal-well. Thus, it is feasible to use one quarter of the reservoir shown in Figure 1b to replace the whole reservoir in order to simplify the derivation process. According to Anderson et al. [33], they verified that when the permeability of region 2 is less than or equal to 500 nD, the contribution of the region beyond fractures can be neglected for the 20-years production. For the case that the distance from the fracture face to permeability boundary (x 1 ) is less than the half distance between fractures (x 2 ), the contribution would be smaller. Meanwhile, Stalgorova et al. [21] also set numerical models to illustrate the negligible contribution of region beyond fractures, and they obtained that the difference of 20-years production is negligible after comparing the results of numerical simulation with and without region beyond fractures. In the work of Heidari et al. [22], they also did not take the region beyond fractures into account. Therefore, the contribution of the region beyond the tips of the hydraulic fractures is assumed to be negligible in this work. In this work, our analytical model is derived under the following assumptions: 1.
The reservoir is homogeneous, isopachous and isothermal in each region.

2.
Flow process is 1-D linear in each region. 3.
Flow is single gas phase. 4.
The high-velocity non-Darcy flow in the hydraulic fracture is considered. 5.
The bottom-hole pressure is constant. 6.
The impact of gravity is neglected. 7.
Gas desorption meets the Langmuir isotherm adsorption equation.

Gas Adsorption/Desorption Effect
In contrast to conventional gas reservoirs, gas adsorption is an important feature of shale gas reservoir. The Langmuir isotherm adsorption equation [34] is widely used to calculate the shale gas adsorption and its expression is as follows: where V is the volume of the adsorbed gas, and P is pressure. V L and P L stand for Langmuir volume and pressure, respectively. When considering the gas adsorption, the effect of adsorption on compressibility of the reservoir is essential. According to Bumb et al. [35], the new compressibility equation can be expressed as, where C f is rock compressibility, C w is water compressibility, C g is free gas compressibility and C gd is adsorbed gas compressibility, S w is water saturation, ρ m is matrix density, B g is gas reservoir volume factor, and φ is the porosity. Another important change is that the compressibility factor is modified by King [36], where P sc is the standard condition pressure, T sc is the standard condition temperature, z is compressibility factor, T is the reservoir temperature.

High-Velocity Non-Darcy Flow Effect
For gas flow in the hydraulic fracture, high-velocity non-Darcy effect is considered in this study. Forchheimer [37] proposed that Darcy's law is inadequate to describe high-velocity gas flow without adding an inertial effect, which is proportional to the square of the flow velocity. To account for the non-Darcy flow effect, an inertial term must be included. The Forchheimer's flow equation is given as, In order to reduce the nonlinearity, the equivalent permeability is introduced to obtain the extended Darcy's law [37], Substituting Equation (6) into Equation (5), the equivalent permeability of hydraulic fracture yields, where β = 4.1 × 10 11 (k iF ) 1.5 (8) among which, β is the non-Darcy flow coefficient, k iF is the hydraulic fracture permeability, k F eq is the equivalent permeability of hydraulic fracture, v is the flow velocity, and ρ is the gas density.

Non-Darcy Flow in the Matrix
For the nano-porous media in shale reservoir, Darcy's law has difficulty describing the actual gas behavior. The gas flow can be classified as four flow conditions, such as continuum flow, slip flow, transition flow, and free-molecular flow. According to the previous publications [38], the four flow conditions can be qualified by the Knudsen number. However, the Knudsen number varies from 10 −3 to 1 in most shale reservoirs. In order to represent the non-Darcy gas flow in matrix, the apparent permeability is presented by the following general form: where k m ∞ is the intrinsic permeability of the porous medium, f (Kn) is the correlation term expressed as a function of the Knudsen number, which is modeled as [38], in which, Meanwhile, for a capillary tube of radius, r, the intrinsic permeability can be derived [39], The Knudsen number Kn is defined as the ratio of the molecular mean free path length and the pore radius in shale matrix, The mean free path can also be calculated, Substituting Equation (14) into Equation (13) on basis of the real gas condition, we can obtain: where µ g is the gas viscosity, z is compressibility factor, r is the pore radius, P is the reservoir pressure, R is the universal gas constant, T is the reservoir temperature, and M is the gas molecular weight.

Derivation of Linearized Gas Diffusivity Equation.
For the flow of shale gas, the gas diffusivity equation will be nonlinear which makes deriving the analytical solution difficult. On one hand, with the reduction in average reservoir pressure, the gas properties like the gas viscosity (µ), gas compressibility (C t ), and gas compressibility factor (z) will change with the pressure. On the other hand, when incorporating the significant mechanisms in shale gas reservoir like gas adsorption and non-Darcy flow effect, the permeability is a variate with pressure rather than a constant.
In order to deal with this problem, the pseudo-pressure and pseudo-time instead of the pressure and time are adopted to linearize the equations. We set a general real gas diffusivity equation in three-dimensional Cartesian coordinate system as example. When the non-Darcy effect is coupled, the k(p) can be calculated by Equation (7) or Equation (9). Certainly, the c t (p) can be replaced by c t * (p) to represent the gas adsorption effect.
∂p ∂t (16) where k(p) is the pressure-dependent permeability, µ(p) is the pressure-dependent viscosity, c t (p) is the pressure-dependent compressibility.
Considering the pressure dependence of the permeability, we define a general modified pseudo-pressure transformation. Surely, k(p) can be calculated by Equation (7) or Equation (9) to couple with the non-Darcy effect. z(p) can be replaced with z * (p) calculated by Equation (4) to consider the gas adsorption.
where m(p) is the pseudo-pressure, k i is the intrinsic permeability, z(p) is the pressure-dependent compressibility factor. Substituting Equation (17) into Equation (16), the right side of the partial differential equation is still nonlinear.
where t a is the pseudo-time, µ i is the initial viscosity, c ti is the initial compressibility. After substituting the Equation (19) into Equation (18), the general linear partial differential equation is derived as, Considering the permeability and compressibility in three regions are different, therefore, the pseudo-time in three regions can be shown as follows.
where t aF , t aF , t a2 are the pseudo-time in the hydraulic fracture region, region 1, and region 2 respectively. Therefore, we can obtain the linearized gas diffusivity equation in three regions respectively. Firstly, we consider the gas adsorption and gas slippage in matrix and the high-velocity non-Darcy flow in hydraulic fracture region. Thus, we use a new total compressibility coupling the gas adsorption effect in diffusivity equation of region 2. Finally, through the linearization by the modified pseudo-pressure and pseudo-time, the non-Darcy flow effect for matrix and fracture are included in the governing linear equations in different regions.

Model Description in Matrix Region (Region 2)
The system of equations based on the conceptual model is presented as follows. For the lowpermeability matrix region (region 2), the diffusivity equation for gas flow is derived as, where m 2 (p) is the pseudo-pressure in region 2, x, y, z are the Cartesian coordinates, t a2 is the pseudo-time for gas flow, c * ti2 and µ i2 are the initial modified total compressibility and viscosity in region 2, k i2 and φ i2 are the permeability and porosity in region 2, respectively.
The initial condition for the region is that the pseudo-pressure of the region is equal to initial pseudo-pressure at t = 0.
where m(p i ) is the initial pseudo-pressure.
The boundary conditions are defined as no-flow at top and bottom of the reservoir. Additionally, the both ends of y-direction can also be regarded as no-flow boundaries.
Due to the plane of symmetry between adjacent fractures, the location x = x 2 is also a no-flow boundary.
The continuity of flux and pressure across the boundaries between the regions is assumed. In Ogunyomi's [27] and Qiu's [29] work, only one is chosen as a boundary condition. Considering the average pseudo-pressure will be adopted, the continuity of flux is chosen as the last boundary condition which can be given by,

Model Description in Stimulated Region Volume (Region 1)
For the stimulated region volume (region 1), the diffusivity equation for gas flow is also expressed as, where m 1 (p) is the pseudo-pressure in region 1, x, y, z are the Cartesian coordinates, t a1 is the pseudo-time for gas flow, c ti1 and µ i1 are the initial total compressibility and viscosity in region 1, k i1 and φ 1 are the permeability and porosity in region 1, respectively. For the whole model, the initial condition remains the same.
The flow in region 1 is also the 1D linear in x-direction, therefore, the outer boundary conditions are identical with those of region 2.
According to the continuity hypothesis, the boundary condition is expressed as,

Model Description in Hydraulic Fracture Region
For the hydraulic fracture region, the diffusivity equation for gas flow is also expressed as, where m F (p) is the pseudo-pressure in hydraulic fracture region, x, y, z are the Cartesian coordinates, t aF is the pseudo-time for gas flow, c tiF and µ iF are the initial total compressibility and viscosity in hydraulic fracture region, k iF and φ F are the permeability and porosity in hydraulic fracture region, respectively. For the whole model, the initial condition remains the same.
With the assumption of the constant bottom-hole pressure, at the location of y = 0, the pressure is the same as the bottom-hole pressure at any time.
The flow in hydraulic fracture region is also the 1D linear in the y-direction, therefore, the outer boundary conditions is expressed as, According to the continuity hypothesis, the boundary condition is expressed as,

Derivation of Analytical Solution.
The diffusivity equations in the mathematical model are all PDEs. In order to obtain the analytical solution, the first step is to transform the sets of equations into ODEs. In this work, we adopted the integral method other than the common Laplace Transform. By integrating over the spatial domain, the spatial dependence is eliminated, which is even feasible to irregular regions as will be demonstrated later. The pressure in different regions is assumed as average value in our model. Therefore, the pseudo-time is supposed to be independent with space. Integrating the equations with respect to spatial coordinates, In Equation (43), we move the pseudo-time outside the spatial integral since the pseudo-time is independent of the spatial coordinates. In order to get a simplified equation, the average pseudo-pressure and the effective pore volume is defined as, where V b is the volume of the region and V p is the pore volume of the region. Equation (43) can be rewritten as, Using the initial condition and boundary conditions, According to the equivalent Darcy's law, Define Therefore, Equation (48) can be rewritten as, where V p2 is the pore volume of region 2, q 2 is the flow rate in region 2.
For region 1, we also use the integration method to deal with the equation, After applying the initial condition and boundary conditions, Equation (51) can be rewritten as, Note that, Defining Substituting Equations (30), (45), (49), and (50) into Equation (48) results in, where V p1 is the pore volume of region 1, q 1 is the flow rate in region 1.
For hydraulic fracture region, we also use the integration method to deal with the diffusion equation, After applying the initial condition and boundary conditions, Equation (56) can be rewritten as, The flow rate in hydraulic fracture can be expressed as, Defining Substituting Equations (42), (48), (58), and (59) into Equation (57) results in, The next step is to substitute the average pseudo-pressure with the relationship between pressure and the flow rate. Since it is assumed that gas flows sequentially from region 2 to region 1 then to hydraulic fracture, a general analytical solution for one-dimensional linear gas flow is derived to solve the problem (details are shown in Appendix A), which is given by: where m(p) is the average pseudo-pressure, q Dn is the dimensionless production from the nth mode.
Combining the assumptions, the average pseudo-pressure in each region can be expressed as, We define the productivity index (J) and transmissibility (T r1F and T r21 ) as, where m 1 (p), m 2 (p) is the average pseudo-pressure in hydraulic fracture region, region 1, and region 2, respectively.
L F is the initial production rate from the hydraulic fracture region.
L 1 is the initial production rate from the region 1.
L 2 is the initial production rate from the region 2.
Substituting Equations (62) . In order to solve the system of ordinary differential equations, an approximate pseudo-time t a is introduced. Therefore, the Equations (50), (55), and (60) can be rewritten as follows.
where q DF n , q D1 n , and q D2 n are the initial production rate from the nth mode in hydraulic fracture region, region 1, and region 2, respectively. t a is an approximate pseudo-time, which is defined as the average of the t a1 , t a2 , and t aF .
Then we define three time-constant parameters as, We can rewrite this set of ODEs in matrix form, where the initial conditions for the system of equations are, After solving Equation (70), we can get the nth flow rate in combination with the initial conditions. The flow rate is the summation of them. By converting the summation to an indefinite integral, the analytical solution in real-time space can be derived as follows (details are shown in Appendix B), q F = β 3 q iF a 3 e λ 1 t a − β 2 q iF a 6 e λ 2 t a + β 1 q iF a 9 e λ 3 t a + In Equation (78), the defined parameters are, where λ 1 , λ 2 , and λ 3 are the eigenvalues of matrix in the Equation (74), a 1 ∼ a 9 are all the values in the eigenvector of matrix in the Equation (75), β 1 , β 2 , and β 3 are all the coefficients. In our model, shale gas in region 2 must firstly flow into region 1 and then into the hydraulic fracture region and finally into the wellbore. Therefore, the flow rate in hydraulic fracture is equal to that in the wellbore. From the final solution, we can find that the flow rate is related to six variables.
Through fitting the production data, these variables can be obtained and then the solution can be used for production analysis and forecasting.

Model Validation with Numerical Models
In order to verify the derived analytical solution, a numerical model is built with the commercial Eclipse reservoir simulator for comparing with the previous physical model, which is one-quarter of the volume around one hydraulic fracture. The numerical model has 27 grid cells in the x-direction, 50 grid cells in the y-direction, and only one grid cell in the z-direction. In order to capture the transient flow towards the hydraulic fracture, local grid refinement in x-direction is constructed. The top view of the model is shown in Figure 2, where the first column of grids represents half-length of hydraulic fracture, and the horizontal well located in the first row of the grids along the x-direction. The blue region represents the region 2 in which the gas adsorption plays an important role, while the red one is the region 1. Table 1 summarizes the input parameters used in the numerical models, which include the reservoir conditions, hydraulic conductivity, gas adsorption, and non-Darcy flow parameters.  In deriving the new analytical solution in this study, the pressure dependence of gas properties is considered using pseudo-pressure and pseudo-time. However, the results from the numerical models are gas rate versus real time. Therefore, the necessary step is to transform the numerical simulation results of gas rate with time into pseudo-time and then fit with the new model. Figure 3 presents the result analysis. According to the previous definition, the relationship between the pseudo-time and real time is shown in Figure 3a. The plot of 1/q vs. square root pseudo-time for regime diagnosis is shown in Figure 3b. The comparison of production rates obtained from numerical simulation and our analytical model is presented in Figure 3c. The blue dotted line represents the relationship between gas rate and pseudo-time, whereas the red one indicates that from the derived analytical solution. It can be seen that the results from the simulator and the analytical solution agree well with each other. Due to the high velocity gas flow in hydraulic facture region, the time constant in hydraulic fracture (0.001 days) is too short to be observed. Therefore, four flow regimes are identified in Figure 4. Regime 1 exhibits a half-slope straight line on the log-log plot and represents the transient linear flow in region 1. The permeability in this region is higher and hence the time constant in this region is shorter and nearly 22 days. Then an exponential curve of regime 2 indicates that the boundary of region 1 is reached, which is called the boundary-dominated flow in region 1 or inner-boundary-dominated flow. After that, the pressure propagates into the region 2. As for regime 3, it still presents an expected straight line with a half-slope signature. In our model, the permeability of region 2 is low-permeability matrix and thus the time constant is relatively longer (209 days). Regime 4 is the outer-boundary-dominated flow, which is affected by the boundary of region 2. Table 2 summarizes the four model parameters after fitting.

Parameter
Value Based on the output parameters, we can predict the values of physical parameters to further validate our model according to the following step-by-step procedure: Step 1: Calculate the productivity index.
As defined above, we can calculate the productivity index J combined with the values of initial rate q iF initial pseudo-pressure m(p i ) and pseudo-bottom-hole pressure m(p w f ).
Step 2: Calculate the transmissibility between fractures and region 1 and region 1 and region 2. Among the output parameters, we can obtain the ratio of transmissibility (T r21 /T r1F , T r1F /J). Considering the calculated value in Step 1, the transmissibility can be calculated (T r21 , T r1F ).
Step 3: Calculate the pore volume of hydraulic fracture region, region 1, and region 2. By transforming the formula that we defined in Equations (71)-(73), the pore volumes of different regions (V pF , V p1 , V p2 ) can be obtained by: Following the above steps, we calculate the physical parameters in the numerical case as shown in Table 3 to compare with the given data. It shows that our model solution is correct within the accepted error bound. Among them, V g and V c express the given volume and calculated volume, respectively.

Irregular Stimulated Region with One Hydraulic Fracture
In this section, three numerical cases with only one hydraulic fracture are designed for investigating the applicability of the analytical model for irregular stimulated region in Figure 4. For the three numerical cases, the input parameters are the same with the Case 1 except for the model dimensions. These three cases are identical with 31 grid cells in the x-direction, 51 grid cells in the y-direction, and only one grid cell in the z-direction, which represent the volume of 214 × 521 × 10 ft 3 . For Case 2, there is a rectangular stimulated region (region 1). As for Case 3, the stimulated region is irregular but symmetrical. In order to better describe the real condition, the stimulated region in Case 4 is designed to be neither regular nor symmetrical. The D-factor in hydraulic fracture is set as 0.0012 to represent the high-velocity non-Darcy flow.
In general, the results from the simulator and the analytical solution fit very well, as shown in Figure 5. There are also four flow regimes in the three cases. Comparing with Case 1 and Case 2, there are deviations for the cases with irregular stimulated region. The deviations are caused by the irregular inner boundary. Due to the irregular shapes of region 1, the time when the pressure propagates to inner boundaries changes with the distance away from the wellbore in different parts of inner boundary. However, the analytical solution is derived based on the regular inner boundary conditions. For the same reason, the variable starting time in regime 2 results in the deviation in the early time of the third flow regime. The flow time in regime 3 is long enough, thus, the curves in this regime fit well with each other at late time. For the three cases, the outer boundaries are identical and therefore the curves in the fourth regime also agree well with each other. According to the results summarized in Table 4, we can observe that the estimated results shown in Figure 5 are acceptable within engineering accuracy. In conclusion, our analytical solution is feasible for both regular and irregular region conditions.

Irregular Stimulated Regions with Several Hydraulic Fractures
In order to further validate the applicability in irregular region of the new derived model, three more numerical cases of multifractured horizontal wells are designed, as shown in Figure 6. These three cases have the identical dimension with 163 grid cells in the x-direction, 63 grid cells in the y-direction, and only one grid cell of 10 ft in the z-direction. The length of the horizontal well is 762 ft with 10 hydraulic fractures equally spaced along the x-direction. For Case 5, the stimulated regions are all the same regular and symmetrical regions. As a comparison, the stimulated regions in Case 6 are all irregular but symmetrical regions. Meanwhile, there is no interference between fractures in Case 5 and Case 6 and the length of the hydraulic fractures is 352 ft in Case 5 and 464 ft in Case 6. As for Case 7, all the stimulated regions around each hydraulic fracture are different, irregular, and asymmetric. The length of 10 hydraulic fractures ranges from 288 to 432 ft and there is interference between fractures. The other input parameters are the same with Case 1.  Figure 7 shows the comparison of gas rate versus pseudo-time obtained from numerical simulations and the new analytical solutions. It can be seen that the results agree very well. Four regimes are also identified. The linear flow time in region 1 based on the fitting results is 9, 7, and 5 days for the three cases, respectively. Then it displays the boundary flow of region 1, which lasts until about the 100th day. The time constant in region 2 is 72, 79, and 91 days, respectively. Finally, the boundary flow of region 2 lasts for hundreds of days. We can obtain that the decline rate of gas production is faster during regimes 1 and 2 and thus the regimes 3 and 4 can last for a longer time. Therefore, for the shale gas reservoir, it is crucial to enhance the volume of stimulated regions to keep a longer high production period. Meanwhile, the contribution of gas adsorption mainly reflects in regimes 3 and 4. It helps to extend the stable production period. According to the output parameters after fitting with the numerical cases, we can calculate the volume of the hydraulic fracture region, region 1, and region 2, as shown in Table 5. Considering that the relative errors meet the requirement of engineering accuracy, our new model is also suitable for the multifractured horizontal well.

Application to Field Case
The previous section demonstrated the accuracy of the derived analytical solution for production analysis. In this section, we apply the method to history matching and forecasting of field data in a shale gas reservoir. The gas well is selected for its relatively long production history and availability of pressure data. The main work flow is as follows: • Apply the given parameters and gas material balance equation to transform the time into pseudo-time and pressure into pseudo-pressure.

•
Make a diagnostic plot of production rate vs. pseudo-time.

•
Analyze the diagnostic plot to identify flow regimes.

•
Fit Equation (75) to the production data to obtain the model parameters τ F , τ 1 , τ 2 , T r21 /T r1F , T r1F /J, and q iF . • Output the model parameters until a satisfactory matching is obtained. • Following the step-by-step procedure to calculate the volume of hydraulic fracture region, region 1, and region 2.

•
Make forecast of production rate with the model parameters.
We chose a horizontal well called Well B-15 with multiple fracture stages from Barnett shale [40]. The general data of the well for analysis is listed in Table 6. The well was produced under constant pressure and even for short early variable p/q the constant pressure was also applicable [41]. Firstly, we transformed the time into pseudo-time and the relationship between pseudo-time and time is shown in Figure 8a. The plot of 1/q vs. square root pseudo-time for regime diagnosis is shown in Figure 8b. Then the log-log diagnostic plot of gas rate versus pseudo-time was obtained, which exhibits a half-slope straight line in Figure 8c. According to the previous analysis, the flow time in matrix is longer and the decline rate in region 2 is slower. Therefore, the half-slope linear flow was diagnosed as the regime 3. Following the work flow, the next step was to match our model to the production data. The results of history matching and forecasting is shown in Figure 8c. The red marks in Figure 8c represent the production data and the green ones represent the analytical solution for history matching. It indicates a good matching. The six model parameters obtained from history matching are summarized in Table 7. Based on the parameters, we forecast the gas rate represented by the black markers in Figure 8c. According to the step-by-step procedure, we calculated the volume of the hydraulic fracture region as 89 ft 3 , the volume of stimulated region as 84.2 MMscf, and the volume of region 2 as 1784.2 MMscf. It is essential for increasing high production period by extending the volume of the stimulated region and decreasing the decline rate in the stimulated region.

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 great 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.

Conclusions
In this paper, we presented a practical analytical model to study the performance of MFHWs from shale gas reservoirs. Numerical models have been used to validate the analytical solutions and an excellent agreement was obtained. The following conclusions are drawn from this work: • A simple rate versus pseudo-time relationship is presented to account for transient linear and boundary-dominated flow periods in shale gas formation.

•
Incorporating the effect of gas adsorption, non-Darcy flow, and slippage flow in the analytical model by defining the modified pseudo-pressure and pseudo-time, accuracy is improved in production forecast in shale gas reservoir.

•
Comparing to the Laplace-transform solution, our analytical model is derived in real-time space and it is unnecessary to undertake dimensionless transformation and numerical inversion. It is more applicable in field scale. • Through the model parameters obtained from history matching the field data, the production rate and cumulative production can be forecasted. In addition, the pore volume of different regions can also be calculated by step-by-step procedure, which was validated to be feasible for the irregular and asymmetric stimulated regions in multifractured horizontal wells. According to the results, the calculation accuracy is less than 10% and meets the engineering requirements.

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

Nomenclature
V the volume of the adsorbed gas, ft 3 P the reservoir pressure, Psi V L Langmuir volume P L Langmuir pressure C f rock compressibility, Psi −1 C w water compressibility, Psi −1 C g free gas compressibility, Psi −1 C gd adsorbed gas compressibility, Psi −1 C t total compressibility, Psi −1 C t * modified total compressibility, Psi −1 B g gas reservoir volume factor z * modified compressibility factor z compressibility factor Z sc standard compressibility factor T reservoir temperature, K T sc standard condition temperature, K P sc standard condition pressure, Psi p average reservoir pressure, Psi k iF hydraulic fracture permeability, mD k F eq equivalent hydraulic fracture permeability, mD k m a matrix permeability, mD k m ∞ intrinsic permeability, mD k i1 region 1 permeability, mD k i2 region 2 permeability, mD r pore radius, ft R the universal gas constant M the gas molecular weight Kn Knudsen number V b volume of the region V p pore volume of the region q iF initial production rate from the hydraulic fracture region q i1 initial production rate from region 1 q i2 initial production rate from region 2 t production time, days t a pseudo-time, days t aF pseudo-time in fracture region, days t aF pseudo-time in region 1, days t a2 pseudo-time in region 2, days t a approximate pseudo-time, days m 2 (p) pseudo-pressure in region 2 m 1 (p) pseudo-pressure in region 1 m F (p) pseudo-pressure in hydraulic region m 2 (p) average pseudo-pressure in region 2 m 1 (p) average pseudo-pressure in region 1 m F (p) average pseudo-pressure in fracture region w/2 half-width of hydraulic fracture, ft x 1 region 1 impact distance, ft x 2 half distance between fractures, ft y e half-length of macro-fracture, ft z e depth of top reservoir, ft z 0 depth of bottom reservoir, ft τ 1 region 1 constant time, days τ 2 region 2 constant time, days τ F hydraulic region constant time, days T r21 transmissibility between region 1 and region 2, STB/D/psi T r1F transmissibility between microfracture and matrix, STB/D/psi J hydraulic fracture region production index, STB/D/psi V pF pore volumes of hydraulic region V p1 pore volumes of region 1 V p2 pore volumes of region 2 α correlation parameter β non-Darcy flow coefficient φ porosity ρ m matrix density, g/cm 3 µ fluid viscosity, cp v gas flow velocity