A practical analytical model for performance prediction in unconventional gas reservoir

Ninety percentage of newly proven natural gas reservoirs of China are mainly unconventional resources, which can be typically developed by multi-stage fracturing horizontal well technology. Two regions typically occur near each fracture after hydraulic fracturing, in which the lower permeability region is considered as the storage source and the higher permeability region as the flow channel. The existed analytical models so far are mainly derived by Laplace transforming. In this paper, an improved practical analytical solution is derived for unconventional gas reservoirs bypassing the Laplace transform and numerical inversion. Through solving material balance equation and adopting the integration, a rate vs. pseudo-time solution in real-time domain can be directly obtained. Five numerical cases are created to verify the accuracy of the proposed analytical solution and the ratio of regular/irregular region pore volume is also proved to be derived reversely by the output parameters, which is significant for the field engineers to evaluate the effect of hydraulic fracturing. Moreover, a field example of a multi-fractured horizontal well in a tight gas reservoir is provided for demonstrate application.


Introduction
Unconventional resources including tight gas and shale gas have become a significant source of global hydrocarbon supply due to advanced reservoir evaluation and drilling/completion/ fracturing technology (Clarkson, 2013;Dejam et al., 2018). As for the ultra-low permeability in unconventional reservoirs, the volume fracturing technology is adopted extensively to improve reservoir properties and achieve economic exploitation. Due to the dense distribution of the preexisting natural fractures, the complex fracture geometry with high permeability will be created around the hydraulic fractures after fracturing operation (Wei et al., 2019;Xia et al., 2019;Liu et al., 2021), which is significant for fluid flow and resource extraction.
In order to make accurate reservoir evaluation, a large number of studies have been conducted for production data analysis. The popular techniques for quantitative production data analysis including empirical methods, numerical methods and analytical models. Firstly, the typical empirical methods, firstly proposed by Arps (1945), have been evolved by many researchers (Fetkovich, 1980;Blasingame et al., 1991;Duong, 2011). These decline curve analysis methods remain commonly used because of their simple and fast empirical regression. However, these methods will result in significant errors especially when applied to unconventional reservoirs. Numerical simulation for unconventional reservoirs has received extensive attention (Cipolla et al., 2010;Ding et al., 2014;Rao et al., 2022). However, it is a time-consuming process and depends on numerous original physical parameters. Compared with the numerical simulation, the analytical methods are models that derive analytically solutions to mathematical models including reservoir behavior as well as flow behavior associated with complex fracture geometry. Therefore, many scholars have focused on analytical models for production data analysis in unconventional reservoirs. Ci-qun (1983) proposed the analytical solution considering the twodimensional pseudo steady state flow from two matrix blocks to the fracture. Lee and Brockenbrough, (1983) derived a tri-linear model to represent the transient flow behavior within the single fracture in an infinite homogeneous reservoir. Wattenbarger et al. (1998) represented an analytical solution accounting for linear flow in fractures with infinite conductivity applied in tight gas wells. Brown et al. (2011) presented a three-region analytical model considering the hydraulic fracture region, the stimulated region including dense branch fractures as well as the outer region containing low permeability matrix beyond the tip of hydraulic fractures. Stalgorova and Mattar, (2012) also introduced a three-region analytical model to consider complex fracture geometry. The outer region in their model adjacent to the stimulated region, which is the only difference from Brown et al. (2011) model. Stalgorova and Mattar, (2013) extended the three-region model to five-region model. Three outer regions are considered in their model. And it notes that the contribution from the outer regions may be negligible except for very late in the production life of the well. Abbassi et al. (2019) put forward a three-region model considering the transient flow among matrix, vugs and fractures. The analytical relationship for calculating the shape factor between different regions is also presented.
The analytical solutions mentioned above were derived analytically based on Laplace-transform technique and the final solution in real-time domain can be obtained by Stehfest numerical inversion (Stehfest, 1970). Development of a straightforward solution for analytical model considering the complex fracture geometry is still lacking. Ogunyomi et al. (2016) and Qiu andLi, (2018) presented the approximate analytical solutions to the dualporosity model and triple-porosity model in real-time domain separately by adopting the integration and average pressure replacement. However, their solutions are applicable for tight oil reservoirs because the pressure-dependent properties in their model are not taken into consideration.
In this paper, our goal is to derive a straightforward solution in real-time domain for unconventional gas reservoirs considering complex fractures. Firstly, we define two regions and the high permeability regions contains hydraulic fractures and complex branch fractures as well as the low permeability region which adjacent to the high permeability region contains the matrix. Then the pseudo-variables (Anderson and Mattar, 2005) are introduced to eliminate the non-linearity caused by the pressuredependent properties of gas. The derivation is based on the integration in real-time domain bypassing Laplace transform and numerical inversion. Lastly, the accuracy of the analytical solution is validated with equivalent numerical model.

Model description
In order to improve gas production in unconventional reservoirs, it is necessary to create as much connection between the unconventional gas reservoir and artificial or natural fractures as possible. Therefore, multistage hydraulic fracturing technology is applied to generate the hydraulic fractures and complex fracture networks to improve gas recovery. As shown in Figure 1A, there is a network of branch fractures forming the flow channel within the tight matrix. By virtue of the highly connected and high conductivity of the network of fractures, the oil/gas stored in ultra-tight matrix can flow into the wellbore. The black arrows represent the flow directions. Obviously, the flow in fractures network is complex and it is impossible to describe each flow process using mathematical models. Therefore, we simplified the complex flow system into two-region system, in which the defined high permeability region (Region 1) includes all the fractures and the low permeability region (Region 2) is the aggregated volume of the matrix directly connected to the high permeability region. Because the contribution of regions beyond hydraulic fractures is validated to be negligible after comparing the results of numerical simulation with and without region beyond fractures (Stalgorova and Mattar, 2012;Abbassi et al., 2019). Figure 1B is a 3D schematic of the reservoir with a multi-stage horizontal well, which is the conceptual model for constructing analytical model.
In this paper, our analytical model is derived analytically based on the following assumption.
1. The reservoir is homogeneous, and isothermal. 2. Flow process is linear in each region. 3. Flow is single gas phase. 4. Flow within the hydraulic fracture is neglected for the high-speed gas flow rate. 5. Flow process is under constant bottom-hole pressure. 6. Pressure at the interface between two regions is constant. 7. Gravity effect is neglected.
Once the gas flow reaches the boundary and reservoir average pressure will decrease, and meanwhile the gas properties such as gas viscosity (μ), gas compressibility (C t ) and gas compressibility (Z) is varying with reservoir pressure. Consequently, the gas diffusivity equation will be non-linear which is impossible to derive the analytical solution directly. In order to deal with this problem, the pseudo-pressure and pseudo-time (Anderson and Mattar, 2005) are adopted to linearize the equation. The equations for pseudo-pressure and pseudo-time can be expressed as, 3 Model development The system of equations based on the conceptual model is presented below. In the low-permeability region (Region 2), the governing equation for gas flow after substituting the Eqs 1, 2 is expressed as: Where m 2 (p) is the pseudo-pressure in region 2, x, y and z are Cartesian coordinates, t a is the pseudo-time for gas flow, k 2 and Φ is the permeability and porosity, respectively, of region 2. The initial condition for the region 2 is equal to initial reservoir pseudo-pressure before the production start.
Where m (p i ) is the initial reservoir pseudo-pressure.
The closed boundary is defined for the top and bottom of the reservoir. Therefore, the outer boundary conditions are regarded as no-flow boundaries.
Considering the symmetry between adjacent hydraulic fractures, the location of x = x 2 is also regarded as no-flow boundary.
Continuity of the flux and pressure across the boundaries between two regions is assumed. Therefore, the inner boundary condition can be written as, Similarly, the governing equation for gas flow in region 1 can be also expressed as, Where m 1 (p) is the pseudo-pressure in region 1, x, y and z are Cartesian coordinates, t a is the pseudo-time for gas flow, k 1 and Φ is the permeability and porosity in region 1 respectively. In regard to the whole flow process, the initial condition of two regions is identical. And the constant bottom-hole pressure is assumed. Thus, the pressure is equal to the bottom-hole pressure at the location of x = x 0 .
Where m (p wf ) is the bottom-hole pseudo-pressure.
The flow in two regions is both 1D linear flow in the x-direction. Therefore, the outer boundary conditions are identical to those of region 2.
zm 1 p zy y y wf ,ye 0 (12) The inner boundary condition based on the continuity of flux at the location of x=x 1 is expressed as: 4 Analytical solution The equations mentioned above used for mathematical model are partial differential equations (PDEs). To solve mathematical model analytically, the system of PDEs must be transformed into the system of ordinary differential equations (ODEs). Considering the Laplace transform and numerical inversion are timeconsuming, the integration method is adopted in this paper to eliminate the spatial dependences and obtain the ODEs directly. Integrating the governing equation with respect to spatial coordinates for region 2: The pseudo-time can be moved outside the spatial integral because the pseudo-time is independent of the spatial coordinates. To obtain a simplified equation, we define the average pseudo-pressure and the effective pore volume as: Where V b is the volume of the region and V p is the pore volume of the region. Substituting the Eqs 15-17 can be rewritten as: With the use of the initial and boundary conditions, Eq. 18 can be simplified as: According to the Darcy's law, the following applies: Therefore, Eq. 19 can be rewritten as: Where V p2 is the pore volume of region 2 and q 2 is the flow rate in region 2.
In region 1, we also apply the integration method to address the diffusion equation as: After applying the initial and boundary conditions, Eq. 22 can be rewritten as, Note that the following applies: Substituting Eqs 9, 16, 34, 38 into Eq. 24 results in: Where V p1 is the pore volume of region 1, q 1 is the flow rate in region 1. Obviously, the next step involves the substitution of the average pseudo-pressure with the relationship between the pressure and flow rate. Since it is assumed that gas flows sequentially from region 2 into region 1, a general analytical solution for 1D linear gas flow is derived to solve the problem (Details are provided in Supplementary Appendix SA), which is given by: Where m(p) is the average pseudo-pressure, q Dn is the dimensionless production of the nth mode. Based on the model assumptions, the average pseudo-pressure in each region can be expressed as: We define the productivity index (J) and transmissibility (T r ) between region 1 and region 2 as: Where m 1 (p) is the average pseudo-pressure in region 1, q i and q i ′ is the initial production rate from the region 1 and region 2 respectively.
Substituting Eqs 27, 28 into Eqs 21, 25 results in: Where q D1n and q D2n is the initial production rate from the nth mode in region 1 and region 2 respectively. We can rewrite this set of ODEs in the following matrix form: Where the initial conditions to solve the system of equations are,

FIGURE 2
Reservoir grid of the numerical model. The blue grids represent the high-permeability region and the gray grids represent the lowpermeability region.
Frontiers in Earth Science frontiersin.org After solving Eq. 33, we can obtain the nth flow rate in combination with the initial conditions. The production rate is the summation of all production rate terms. By converting the summation to an indefinite integral, the analytical solution in realtime domain can be derived below (Details are provided in Supplementary Appendix SB), Where the defined parameters are expressed as: Where λ 1 and λ 2 are both the eigenvalue in the matrix, a 1 and a 2 are both the factor in the eigenvector in the matrix, τ 1 and τ 2 is the time constant in region 1 and region 2 respectively. Based on the final solution in Eq. 36, it reveals that the production rate is related to four variables (productivity index (J), transmissibility between region 1 and region 2 (T r ), time constant in region 1 (τ 1 ) and time constant in region 2 (τ 2 )). By fitting production data, four variables can be determined and the analytical solution can be applied for performance prediction.

Validation against numerical cases
To verify the derived analytical solution, one numerical model contains 35 grid cells in the x-direction, 50 grid cells in the y-direction and only one grid cell in the z-direction is built with Eclipse reservoir simulator and the model parameters are summarized in Table 1. Considering the symmetry of the conceptual model and computational convenience, one-quarter of the region 1 and region 2 around one hydraulic fracture exhibited in Figure 1A is extracted for comparison. A top view of the numerical model is shown in Figure 2, where the first column of grids represents the half-length of the hydraulic fracture and the first row of grids in the x-direction represents the horizontal wellbore only connected with the region 1. The gray region represents region 2, while the blue region represents region 1.
Pressure-dependent properties is considered and linearized using the pseudo-pressure and pseudo-time during the derivation of the new analytical solution. It is essential to transform the results obtained with numerical models into production rate over pseudo-time and then fit the transformed results with the derived analytical solution. A comparison of the production rates obtained with the numerical simulation and our analytical model is shown in Figure 3. The blue dotted line indicates the relationship between the gas rate and pseudo-  time whereas the red dotted line indicates that determined from the derived analytical solution. The simulation and analytical results agree well. As shown in Figure 3, four flow regimes are identified. Regime 1 exhibits a half-slope straight line in a log-log plot which represents the transient linear flow in region 1. The permeability in this region is relatively high and hence the time constant in this region is relatively short and equals nearly 4 days. Then, the exponential curve of Regime 2 represents the boundary of region 1 is reached, which is referred to as boundary-dominated flow in region 1 or inner-boundary-dominated flow. Regime 3 exhibits the expected straight line with a half-slope signature. In our model, the permeability of region 2 is low and thus the time constant is relatively long (246 days). Regime 4 experiences  Frontiers in Earth Science frontiersin.org outer-boundary-dominated flow, which is caused by the boundary of region 2. Table 2 summarizes the four variables after fitting.

Validation for the pore volume ratio
According to Eqs 37, 38 the time constants in region 1 and region 2 are defined to relate strongly with the pore volumes of regions. When Eq. 37 is divided by Eq. 38, we can derive the pore volume ratio of two regions as below: Substituting the model parameters in Table 2 into Eq. 43, the pore volume ratio can be obtained. We compare it to the given value from the numerical model which is proved correct within the accepted error bound.
The calculated pore volume ratio is nearly equal to the given value in Case 1. To further verify the accuracy of Eq. 43, four more numerical cases are conducted based on the previous numerical model parameters. Table 3 summarizes all the output values from these four numerical cases. Through comparing the given pore volume ratio to the calculated results, the comparison reveals that the accuracy of the calculated results is higher than 95% and thus the model parameters from our new analytical solution can be applied to estimate the pore volume ratio.

Application to field example
In this section, we apply the practical analytical solution to field data. The Mexico gas well is a 'near textbook' quality tight gas well with long production history and high quality data, which is operated under a constant bottom-hole pressure of 920 Psi. Several analyses have shown that the reservoir has an ultra-low permeability less than 0.001 md and there is only one well in this field. The original production data is from Blasingame et al. (2007). First, we transform time into the pseudo-time and then the log-log diagnostic plot of gas rate versus pseudo-time is obtained, which exhibit a half-slope straight line and a nearly unit-slope line in Figure 4A indicating linear and boundarydominated flow. Since the linear flow lasts for nearly 10,000 days, this linear flow regime is identified as Regime 3 and the boundary-dominated flow represents Regime 4. The next step involved matching our model to the production data.
The history matching and forecasting results is shown in Figure 4B. The red marks in Figure 4B indicate the production data and the green marks indicate the analytical solution of the history match. A suitable matching degree is revealed. The four model parameters obtained from the history matching process is summarized in Table 4. Based on these output parameters, we forecast the gas rate until 30,000 days, indicated by the black markers in Figure 4B. The more intuitive results are shown in Figure 4C by Cartesian plots. Meanwhile, we can calculate the ratio of region pore volume by Eq. 43 as V p2 /V p1 =104. The ratio can also be used to explain why only two flow regions can be observed rather than four regimes as seen in the previous numerical models. From the ratio, we can obtain that the volume of high permeability region is far smaller than the matrix region and thus the flow in region 1 with high permeability is too fast to be observed.

Conclusion
In this paper, a practical analytical model is presented for performance prediction in unconventional gas reservoirs under constant bottomhole pressure. Numerical models were employed to verify the accuracy of derived analytical solution and an excellent agreement was revealed. The following conclusions are drawn: (1) The proposed analytical solution is derived in real-time domain bypassing Laplace transform and numerical inversion, which is highly suitable for field applications.
(2) Pressure-dependent properties is considered and linearized using the pseudo-pressure and pseudo-time in the analytical model and it is applicable for performance prediction in unconventional gas reservoirs. (3) The pore volume ratio of different regions can be calculated reversely and the relative error is less than 10%, which is helpful for rapidly evaluating the effect of hydraulic fracturing.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions
Conceptualization; methodology and investigation; writing-original draft preparation; project administration; funding acquisition-KQ.

Funding
This study was supported by Basic Research Project from Jiangmen Science and Technology Bureau (Grant No. 2220002000356).