Transient pressure behavior of horizontal well in gas reservoirs with arbitrary boundary

Transient pressure analysis is a crucial tool to forecast the production performance during the exploration and production process in gas reservoirs. Usually, a regular shaped outer boundary is assumed in previous studies for well-testing analysis, which is just a simplification of practical cases and cannot reflect the actual boundaries of reservoirs. In this paper, a mathematical model is established to analyze the transient pressure behaviors of a horizontal well in an arbitrarily shaped gas reservoir. Dimensionless treatment, Laplace transformation, and boundary element method are applied in solving the model, which is verified by comparing with the results from the source function method. Based on the Stehfest numerical inversion method, the models of single-porosity media and dual-porosity media are solved respectively. Then, the time-domain curves of pseudo pressure and its derivative are obtained, and the flow regimes are identified. Finally, the impacts of some critical parameters on pressure transient behaviors are analyzed, including storativity ratio, interporosity coefficient, well length, and well orientation. This paper presents an effective way to handle complex external boundary problems in gas reservoirs.


Introduction
The exploration and development of natural gas plays an irreplaceable role in building a green and clean energy system (Zou et al., 2018). Horizontal wells can significantly accelerate production and increase recovery efficiency of gas reservoirs, especially in low permeability formations (Biryukov and Kuchuk, 2015;He et al., 2020;Zheng et al., 2020). The application of this technique belongs to the innovation of drilling technology and the improvement of engineering techniques (Al Rbeawi and Tiab, 2011). However, corresponding theoretical studies of pressure transient behaviors in gas reservoirs are far behind practices, which restrict the efficient development of these reservoir resources. The outlines of the outer boundary are usually simplified as circle or rectangle in well testing for closed outer boundary reservoirs. As we know, the practical reservoir shape cannot be very regular as a circle or a rectangle, but with odd outlines. Thus, it is necessary to study the effects of boundary geometry on transient pressure behaviors, which is helpful to increase the precision of well testing interpretation and enhance the producing efficiency in gas reservoirs.
The source function method is one of the most common methods to solve percolation issues of oil and gas reservoirs and analyze transient behaviors of wells (Ozkan and Raghavan, 1991;Pan et al., 2007;Qin et al., 2019;Raghavan et al., 1997). Chu and Shank (1993) used the source function method to analyze the pressure performance of the fractured wells with finite conductivity and uniform flow in the radial composite reservoir. Chen and Raghavan (1997) analyzed transient pressure behaviors of horizontal wells in a rectangular oil reservoir by the source function method in the Laplace domain. Wan and Aziz (1999Aziz ( , 2002) established a new 3-D seepage model for a horizontal well and utilized Fourier transformation to convert 3-D seepage problems into 2-D seepage problems. Later, this approach was widely adopted in the well-testing analysis of complex structural wells. Medeiros et al. (2007) studied the well testing model of a horizontal well in rectangular composite formation, and then they obtained its semi-analytical solution. Wang et al. (2012) established the inflow model of highly deviated wells in anisotropic reservoir according to the coupling of the reservoir and wellbore. And the inflow performance in an infinite slab anisotropic reservoir with impermeable top and bottom boundaries were studied. Wang and Yi (2018) presented a coupling mathematical model for a fractured horizontal well in triple media carbonate reservoir with three kinds of outer boundaries by conceptualizing vugs as spherical shapes. Then, a semi-analytical solution was obtained in the Laplace domain, and the transient pressure responses of fractured horizontal wells were analyzed.
However, the source function method can only deal with seepage problems in oil and gas reservoirs with regular-shaped external boundaries, such as infinite, circular, or rectangular external boundaries. Still, it is tough to use the source function method to study transient behaviors for reservoirs with complex outer boundaries. As a result, the numerical methods need to be employed to obtain the characteristics of transient pressure for horizontal wells in gas reservoirs with arbitrary boundaries. The boundary element method (BEM), as one of the numerical methods, has been proved to be suitable to solve the transient flow model of a complicated well and boundary models, especially for the reservoirs with irregularly shaped boundaries. Compared with the finite difference method, the BEM does not need to divide the reservoir into many grids and spend much time to solve the pressure of all grids at each time step. The BEM just need to divide the reservoir boundary into some elements, and treat every element as a plane source. Therefore, the BEM just need to solve the pressure solutions of the boundary elements. Kikani and Horne (1989) pointed out that the numerical divergence and grid orientation problems of the BEM can be ignored, due to the fundamental solution of the BEM is analytical. Zhang and Zeng (1992) analyzed the transient pressure behavior of a vertical well in arbitrarily shaped oil reservoirs by the BEM. Jongkittinarukom and Tiab (1998) combined the methods of the numerical well testing technique and the BEM to study the pressure transient behaviors in heterogeneous oil reservoirs. Their study showed that the pressure derivative values obtained by the BEM were more accurate than those from the finite difference method. Wang and Zhang (2009) studied the transient seepage problem of a vertical well in a gas reservoir with irregular form. Zhao et al. (2016) used the BEM to analyze pressure transient of multiple fractured horizontal wells (MFHWs) in the gas reservoirs with arbitrary outer boundary, and the solutions obtained by the BEM fit well with the solutions derived from the semi-analytical method and well testing software, respectively. Zhang et al. (2017) applied the BEM for transient pressure analysis of a fractured vertical well in a rectangular composite coal bed gas reservoir. Besides, other scholars presented BEM solutions for a well located in geometrically complex reservoirs (Idorenyin and Shirif, 2017;Zhang et al., 2018;Zhao et al., 2017). Zhao et al. (2018) firstly utilized the BEM to analyze the pressure behavior of the MFHW with the consideration of stimulated reservoir volume (SRV). In their model, the SRV and the outer boundary of the reservoir are treated as rectangles. And later, the same method for wells in the composite reservoir is also analyzed by Wu et al. (2018) and Chen et al. (2018).
As we can see from the literature review, the source function method can only deal with seepage problems of wells in conventional shaped reservoirs, which cannot deal with the seepage problems of wells in arbitrarily shaped reservoirs. Given that the formation geology is very complex, practical reservoir shape cannot be regular. Although the BEM can handle seepage problems in complex reservoir geometries, the relative studies are mainly focused on vertical wells or vertically fractured wells. In contrast, transient pressure research on horizontal wells in oil or gas reservoir is rare. Since the long horizontal well technique is a necessary measure for commercial development of a gas reservoir, it is essential to understand the transient seepage mechanism of long horizontal well in arbitrarily shaped reservoirs.
In this paper, a seepage model of a long horizontal well in irregular shaped gas reservoirs was developed, based on gas continuity equation, motion equation, and state equation. The mathematical model is solved by the BEM method, and the results are verified by the source function method. Then, different flow patterns were distinguished, and sensitivity analysis of parameters is conducted through the proposed mathematical model.

Physical model and assumptions
A horizontal well in an irregular gas reservoir is considered in this study. As shown in Figure 1, a horizontal well is located in the reservoir, whose external boundaries can have arbitrary shapes. In order to develop the mathematical model, some assumptions are listed as 1. The external boundary of the gas reservoir is impermeable with an arbitrary shape. 2. The upper and lower boundaries of the formation are also impermeable, with a uniform formation thickness h.
3. The effects of gravity, capillary force, and rock compressibility are ignored. Meanwhile, the seepage problem is deemed as the isothermal process. 4. The initial reservoir pressure is uniform, with the value of p i . The effect of the wellbore storage and skin factor is considered. 5. The horizontal wellbore of the well has infinite conductibility.

Governing equation
According to mass conservation theory, the governing equation of gas flow in a reservoir can be obtained by combining the gas diffusivity equation, the motion equation of Darcy's law, and the real gas state equation. The expression is as follows, Where k is permeability, m 2 ; l is viscosity, PaÁs; p is pressure, Pa; Z is gas deviation factor, dimensionless; / is porosity, dimensionless; T is reservoir temperature, K; p sc is the pressure at standard conditions, Pa; T sc is the temperature at standard conditions, K; q l is line source strength, m/s; Q is any point in the field; Q 0 is point source location in the field and d is Dirac function, which has the following properties Introducing pseudo pressure (Al-Hussainy et al., 1966) to simplify the above nonlinear diffusivity equation Then, the governing equation of equation (1) can be expressed as Where c g is gas compression coefficient, which can be calculated by According to assumption 4, the initial condition can be given as The impermeable outer boundary can be defined by Neumann boundary condition as Where n is the normal direction at the outer boundary C. By using the dimensionless variables in Table 1, the seepage governing equation can be dimensionalized. Meanwhile, due to formation is fully penetrated by vertical fractures, the 3-D seepage problems can be simplified into 2-D problems. Thus, the governing equation can be transformed into Where q sct is total production rate at standard conditions, m 3 /s; q v is volume flux from point source per unit time; r D and r 0 D are the locations of an arbitrary point and point source in the 2-D field, respectively.
Imposing Laplace transformation to equation (9), we can obtain

Dimensionless variables Expression
The above equation (10) is the governing equation of gas seepage behavior in a homogeneous formation. If the dual-porosity model is considered, the parameter s needs to be replaced by f(s), which has the following form, For homogeneous reservoirs, For dual-media reservoirs (Al-Rbeawi, 2019),

Model solution by boundary element method
To solve the governing equation of equation (11), a boundary element method was employed in this paper, which discretized the arbitrary shaped outer boundary into small segments. Before its application, it's necessary to give the fundamental solution, E(P, Q, s), which satisfies equation (14) in the 2-D domain, @ 2 EðP; Q; sÞ @x 2 D þ @ 2 EðP; Q; sÞ @y 2 D À fðsÞEðP; Q; sÞ þ 2pdðP; QÞ ¼ 0 The fundamental solution is expressed as follows, Any point in the formation satisfies the governing equation of equation (11), Multiplying equation (14) (18) Subtracting equation (18) from equation (17), we have 2p m D ðP; sÞdðP; QÞ À 2p fðsÞ q D EðP; Q; sÞdðP; QÞ Integrating equation (19) over the reservoir domain X, we have With the characteristic of the Dirac function of equation (3) and Green's second identity, equation (20) To make equation (21) more general, a parameter h is added in equation (21) Where h equals to 1/2 for point Q on smooth segments of the boundary C; h equals to a/2p for point Q on non-smooth segments of boundary C; h equals to 1 for point Q inside the domain X; a is the angle of the tangent of the boundary segment at point Q.
As can be seen above, the integral boundary equation has been derived in the Laplace domain for both homogeneous and dual media reservoirs. If we change the observation point on the horizontal well, the dimensionless bottom-hole pressure m wDN can be obtained, regardless of wellbore storage and skin effects. The sequence of discretized segments satisfied the right-hand rule, which can be seen in Figure 2. The solution procedure in detail can refer to Wang and Zhang (2009) and Zhao et al. (2017).
Duhamel's principle is employed to take wellbore storage and skin effects into accounts, which is Where m wD and m wDN are bottom hole pressure with and without considering wellbore storage and skin effects; S skin is skin factor; C D is wellbore storage factor.

Numerical inversion algorithm
All the solutions derived previously are in the Laplace space. To plot their type curves, they need to convert into the real space. For a complex format, it is difficult to get a primitive function through the Laplace inverse transformation or contour integral. In engineering applications, what people usually care about are the characteristics of the pressure type curves in different flow stages, and the real data can be explained through matching with the type curves. The difficulty in inversion through a contour integral and the complexity of image functions result in difficulty in the acquisition of their primitive functions. This constrains an application of the Laplace transformation in the engineering field. For a long time, people had been seeking a satisfactory numerical inversion method to avoid shortages of analytical inversion. Until the 1970s, various numerical inversion algorithms were developed, which accelerated the application of the Laplace transformation in engineering. In petroleum engineering, there are two main methods in numerical inversion. One is the Stehfest inversion algorithm (Stehfest, 1970), which is based on a function probability density theory, and the other one is the Crump inversion algorithm, which is based on the Fourier series theory. Due to its simplicity and programming friendliness, the Stehfest inversion algorithm is the most common one in engineering applications.
For the Stehfest numerical inversion algorithm, assuming that F(s) is the Laplace transform of a time function f(t), then, Therefore, the Stehfest numerical inversion algorithm to invert a solution from the Laplace space to the real space solution is where s j ¼ ln2 t j, and According to the above inversion method, the image function f(t) can be calculated by equation (26) when N and t are known. It is worthy of noting that N must be an even number, and its value can significantly affect the calculation accuracy. The N value can be decided for different types of functions through practical calculations. In most cases, the N values of 8, 10, and 12 are appropriate. In this paper, the N is set to be 8.

Model verification
To verify the model with irregular outer boundary, a simplified case is adopted here to compare with the results from different methods. As mentioned before, the source function method is a good tool to deal with seepage problems in oil or gas reservoirs, and the results from the source function method are very reliable. However, this method cannot handle seepage problems in reservoirs with complex outer boundaries, such as the arbitrary shaped outer boundary. Thus, in this part, the transient pressure behavior of a horizontal well in a gas reservoir with a rectangular shaped outer boundary (seen in Figure 2) is used to verify our model. The length of the horizontal well is taken as the value of 200 m. The length of the gas reservoir is 1000 m, and the width of the reservoir is 400 m. The wellbore storage coefficient and skin factor are respectively 0.001 and 0.1. Meanwhile, the formation is deemed as singe media, where the single porosity model was employed in this case. After calculation, the results from the source function method and our boundary element method are compared.
According to the souse function method proposed by Ozkan and Raghavan (1991), the continuous point source function for a gas reservoir in an infinite outer boundary is as follows.
where n n ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi L D 2 b n 2 þ sfðsÞ p , b n ¼ np, n ¼ 0,1,2,3Á Á Á. As shown in Figure 3, the calculation results of transient pressure behaviors of a horizontal well in box gas reservoir are in good agreement between the traditional source function method and our new boundary element method (See Appendix 1 for the original data). The good agreement between the results from two different methods reflects the correctness of our model. In the next part, our model will be used to conduct the identification of flow regimes and sensitivity analysis.

Flow regimes identification
In this section, a 250 m horizontal well in a circular gas reservoir with a radius of 12000 m. Other parameters are shown in Table 2, which were used in all the following cases, and the physical model is shown in Figure 4. As can be seen from Figure 5, there are seven main flow regimes for transient pressure behaviors of a horizontal well in a gas reservoir. ‹ The first flow regime is the wellbore storage flow regime. In this flow regime, the curves of pseudo pressure and its derivative are overlapped lines of a unit slope. › The second flow regime is the transition flow period, which is controlled by the skin factor. fi The third is the early linear flow period, where the pseudo pressure derivative curve appears a straight line with 0.5 slopes. The schematic of the early linear flow for a horizontal well is showed in Figure 6. fl The forth flow regime is the early radial flow period, where there is a horizontal line in the pseudo pressure derivative curve. The position of the horizontal line is related to the length of the horizontal well.°T he fifth flow regime is the interporosity flow period. In this period, the pressure in fractures is much lower than that in the matrix, resulting in fluid supply from the matrix system to the fracture system. -The sixth flow regime is the later radial flow period, where there is a horizontal stage with a value of 0.5 in the pseudo pressure derivative curve. † The last flow regime is the boundary dominated flow period, where the curves of pseudo pressure and its derivative are overlapped straight lines with a unit slope in the type curves. Storativity ratio of the fracture system The effects of the storativity ratio on transient pressure behaviors of a horizontal well in a gas reservoir can be seen in Figure 7. The storativity ratio mainly affects the shape of the type curve  in the interporosity flow period. With the storativity ratio increasing, the storage capacity in the fracture system becomes stronger, resulting in the hump in the type curve becomes lower. That is because a bigger storativity ratio, which means a stronger capacity of fluid supply from the fracture system, causes smaller pressure depletion for constant well production.

Interporosity coefficient
The interporosity coefficient can affect the beginning time of the interporosity flow period, as shown in Figure 8. The interporosity coefficient determines the location of the concave in pseudo pressure derivative curves. The bigger the interporosity coefficient is, the easier the interporosity flow could occur, and the earlier the supply from the matrix system to fracture system occurs. Therefore, a bigger interporosity coefficient brings an earlier interporosity flow period. Contrarily, a smaller interporosity coefficient may cover the appearance of the later radial flow period in pseudo pressure derivative curves.

Length of horizontal well
As we can see in Figure 9, the impact of horizontal well length on the pressure behaviors mainly take place in the wellbore storage effect stage and later radial flow stage. The longer the horizontal well is, the harder the later radial flow period could be observed for constant rate production in the reservoirs with the same size.

Horizontal well orientation
In this part, the effects of horizontal well orientation on transient pressure behaviors are considered. As can be seen from Figure 10, the angle between the horizontal well and the long axis of an elliptical-shaped reservoir are 0 , 30 , 60 for case 1, case 2, and case 3, respectively. As we can see in Figure 11

Transient pressure behavior in an arbitrarily shaped reservoir
Subsequently, the transient pressure behavior of a horizontal well in an irregular reservoir is analyzed. The outer boundary shape of the gas reservoir can be seen in Figure 12. As mentioned, conventional means like Green function and source function methods cannot solve the transient pressure behavior of a horizontal well in such arbitrary shaped reservoirs. Thus, we calculated the pseudo pressure and its derivatives using our proposed model. As we can see in Figure 13, even though the external boundaries of the reservoirs are odd, our model can still obtain pressure transient behaviors conveniently.

Conclusions
In this paper, a mathematical model is established to analyze the transient pressure behaviors of a horizontal well in an arbitrarily shaped gas reservoir. The model is solved by the boundary element method, and type curves are plotted and further discussed. Through the above discussion, we can obtain the following conclusions: 1. The boundary element method can be applied to solving the mathematical model of a horizontal well in a gas reservoir with arbitrary shaped outer boundaries, which was verified by comparing it with the source function method. 2. Seven flow regimes can be identified for transient pressure behaviors of a horizontal well in a gas reservoir. They are wellbore storage regime, transitional flow period, early liner flow stage, early radial flow stage, interporosity flow stage, later radial flow stage, and boundary dominated stage. 3. The storativity ratio mainly determines the width and depth of concaves in dimensionless pseudo-pressure derivative curves during the interporosity flow period. The smaller the storativity ratio is, the wider and deeper the concave will be. The interporosity coefficient determines the beginning time of the interporosity flow period. The larger the interporosity coefficient is, the earlier the interposity flow period will happen. 4. The influences of horizontal well length and horizontal well orientation on transient pressure behaviors are obtained and analyzed. According to the sensitivity analysis, the closer the well locates from the boundary, the deviation will be more serious compared to the conventional model. For real situations, if a horizontal well drilled in a blocked reservoir with an arbitrary outer boundary, and then the model proposed in this paper will be more accurate. 3.42Eþ04 3.43Eþ04 2.51Eþ07 3.94Eþ03 3.95Eþ03 3.63Eþ08 5.68Eþ04 5.67Eþ04 3.98Eþ07 6.24Eþ03 6.25Eþ03 6.03Eþ08 9.41Eþ04 9.39Eþ04 6.31Eþ07 9.89Eþ03 9.90Eþ03 1.00Eþ09 1.56Eþ05 1.55Eþ05 1.00Eþ08 1.57Eþ04 1.57Eþ04