Transient Behavior of Vertical Commingled Well in Vertical Non-Uniform Boundary Radii Reservoir

: Transient behavior analysis, including rate transient analysis (RTA) and pressure transient analysis (PTA), is a powerful tool to investigate the production performance of the vertical commingled well from long-term production data and capture the formation parameters of the multilayer reservoir from transient well testing data. Current transient behavior analysis models hardly consider the e ﬀ ect of vertical non-uniform boundary radii (VNBR) on transient performances (rate decline and pressure response). The VNBR may cause an obvious novel radial ﬂow regime and rate decline type, which can easily be mistaken as a radial composite reservoir. In this paper, we present a VNBR transient behavior analysis model, the extended vertical uniform boundary radii (VUBR), to investigate the rate decline behavior and pressure response characteristics through diagnostic type curves. Results indicate that the dimensionless production integral derivative curve or pressure derivative curve can magnify the di ﬀ erences so that we can diagnose the outer-convex shape and size of the VNBR. Therefore, it is signiﬁcant to incorporate the e ﬀ ects of VNBR into the transient behavior analysis models of the vertical commingled well, and the model proposed in this paper enables us to better evaluate well performance, capture formation characteristics and diagnose ﬂow regimes based on rate / pressure transient data.


Introduction
With the depletion of conventional oil and gas resources, more and more complex, multilayer oil and gas reservoirs have been discovered with the improvement of exploration technology. The well production performance and reservoir evaluation are of great significance for enhancing recovery in multilayer oil and gas field development [1][2][3]. Both production decline analysis and well testing analysis are important dynamic technologies for petroleum engineers to evaluate well performance [4][5][6]. During the past years, many researchers have investigated the production decline behavior and pressure response characteristics of the vertical commingled well. Their main focus was four types of analysis: transient behavior analysis (including pressure transient analysis (PTA) [7][8][9][10][11][12] and rate transient analysis (RTA) [12][13][14][15]), vertical-heterogeneity analysis [16][17][18] and boundary effect analysis [19][20][21], for a commingled well in a multilayer reservoir.
For the PTA and RTA for the vertical commingled production well, Satman used the discontinuous inner boundary radii to characterize the front of the displacement fluid during the enhanced recovery processes, such as steam flooding and gravity override effect [7,8]. Tariq and Ramey presented a PTA model considering the wellbore storage and skin factor in the multilayer reservoir [9]. Anbarci et al. presented a PTA model to determine the front location in a multilayer reservoir [10]. Ehlig-Economides and Joseph validated a general analytical solution of PTA for an n-layer, composite, crossflow reservoir system [11]. Fetkovich et al. conducted studies on the RTA and PTA for a two-layer non-communicating reservoir, and the pressure response curves showed that these layers can be combined into one layer if the layers have the same physical properties [12]. El-Banbi and Wattenbarger developed an RTA model to match the gas well production data for a three-layer tight gas reservoir [13]. Villanueva-Triana and Civan developed an RTA model for three-layer commingled reservoirs involving formation crossflow and general external flow boundaries [14]. Milad simulated a commingled well RTA model for multilayer shale gas reservoirs, in which the pressures and flow rates are not constant in various zones [15].
For the vertical-heterogeneity of the multilayer reservoir, Aly et al. examined the effect of unequal initial layer pressure on RTA and PTA models [16]. However, their objective was the permeability and skin factor rather than boundary radii of an individual layer. Jordan and Mattar compared the PTA curves of the composite reservoir with that of the multilayer reservoir [17]. His modified composite reservoir demonstrated that the composite reservoir pressure behavior would be like that of a multilayer reservoir if the vertical permeability is very small. Lodon et al. provided a two-layer semi-analytical PTA model of a commingled reservoir with different boundary distances of each layer, and their model shows that the PTA curve of the two-layer reservoir is located between that of the two-layer reservoir with the same minimum and maximum boundary distance [18]. However, their research is limited to the two-layer reservoir and cannot provide a quantitative relation of pressure behavior between the total reservoir and individual layers.
For the boundary effect of the reservoir, Larsen generated an exact Laplace-transformed solution including a great variety of boundary configurations. Results show the PTA of a square boundary reservoir is the same as that of a circle boundary reservoir [19]. For a well near a liner leaky boundary in the two-layer reservoir, Boussila et al. revealed that there is a quantitative relationship between the pressure derivative and the mobility ratio during the boundary flow stage [20]. For the tidal flat layered carbonate reservoir in Sichuan Basin, Shi et al. characterized the pressure transient behavior by a layered reservoir with the vertical inhomogeneous closed boundary [21]. Their results showed that the shape of the closed boundary has an intuitive inverse proportional relationship on the pressure derivative curve during the radial flow regime.
The existing research mainly focuses on the analysis and acquisition of the parameters of layered reservoirs from the production data, and few works discuss the equivalence and substitutability of the boundary effect of multilayer reservoirs [22]. The objective of this paper is to investigate the influence of vertical non-uniform boundary radii (VNBR) on transient behavior and explore the corresponding relation between them. This paper extends the boundaries issues of the multilayer reservoir by a VNBR transient behavior analysis method. Some interesting and valuable phenomena, such as the "equivalent seepage volume effect" and the "equivalent composite zone effect", are identified, which give a new perspective on the flow regimes of the vertical commingled well in a VNBR reservoir and provide a quick approach to capturing the formation parameters from transient well test data. The transient behavior analysis model presented in this paper provides a tool to characterize the production performance and rate decline type from the long-term production data. Figure 1a shows a vertical commingled well in the VNBR reservoir. The following assumptions are considered to yield the mathematical model:

1.
A vertical production well is located at the center of a multilayer reservoir formation with constant initial reservoir pressure (p i ) in any layer (j). The formation has a circular-closed boundary with non-uniform boundary radii in the vertical direction, but any layer j has an individual constant thickness (h j ) and boundary radius (r ej ).

2.
The formation is considered to be isotropic and homogeneous. Any layer has individual physical parameters, including the permeability (k j ), porosity (ϕ j ) and total compressibility (ctj).

3.
The formation is fully penetrated by the vertical production well with height (ht) and radius (rw). During the production stage, both the pressure (p j ) and rate (q j ) of any layer are different. The bottom-hole pressure (p w ) and wellbore rate (q) are variable with the production time.

4.
The porosity of the formation is filled with a single-phase fluid with constant viscosity (µ). Radial flow regime is assumed in any layer (see Figure 1b,c). The isothermal volume factor (B) is descript the volume change of the natural gas from the wellbore to the wellhead. 5.
The gravity and temperature effects are ignored.
Energies 2020, 13, x FOR PEER REVIEW 3 of 14 Figure 1a shows a vertical commingled well in the VNBR reservoir. The following assumptions are considered to yield the mathematical model: 1. A vertical production well is located at the center of a multilayer reservoir formation with constant initial reservoir pressure (pi) in any layer (j). The formation has a circular-closed boundary with non-uniform boundary radii in the vertical direction, but any layer j has an individual constant thickness (hj) and boundary radius (rej). 2. The formation is considered to be isotropic and homogeneous. Any layer has individual physical parameters, including the permeability (kj), porosity (φj) and total compressibility (ctj). 3. The formation is fully penetrated by the vertical production well with height (ht) and radius (rw).

Physical Model
During the production stage, both the pressure (pj) and rate (qj) of any layer are different. The bottom-hole pressure (pw) and wellbore rate (q) are variable with the production time. 4. The porosity of the formation is filled with a single-phase fluid with constant viscosity (μ). Radial flow regime is assumed in any layer (see Figure 1b,c). The isothermal volume factor (B) is descript the volume change of the natural gas from the wellbore to the wellhead. 5. The gravity and temperature effects are ignored.

Mathematical Model
The dimensionless flow equations of any layer j in Laplace domain are (see the detailed derivation in Appendix A): The definition of the parameters is provided in Table A1, and the boundary conditions are presented in Appendix A.
Substituting Aj & Bj, obtained by solving Equation (A11) into Equation (A12), the dimensionless rate solution of RTA model is obtained as: where the parameters are presented in Appendix B.

Mathematical Model
The dimensionless flow equations of any layer j in Laplace domain are (see the detailed derivation in Appendix A): The definition of the parameters is provided in Table A1, and the boundary conditions are presented in Appendix A.
Substituting A j & B j, obtained by solving Equation (A11) into Equation (A12), the dimensionless rate solution of RTA model is obtained as: where the parameters are presented in Appendix B.
Energies 2020, 13, 2305 4 of 13 Stehfest provides the mathematical relationship of a dimensionless solution between Laplace domain and real-time domain [23]. Blasingame et al. defined three parameters to characterize dimensionless transient rate as follows [24]: where qDd denotes the normalized dimensionless decline rate, qDdi represents the normalized dimensionless decline rate integral and qDdid is the normalized dimensionless decline rate integral derivative. tDd is the normalized dimensionless decline time, which can be calculated by: According to the research results of Van-Everdingen and Hurst, the dimensionless pressure solution of PTA mode can be obtained through the dimensionless rate solution of the RTA model in Laplace domain [25]: Based on Duhamel's superposition principle, the dimensionless pressure solution of the PTA model in Laplace domain with the effects of skin factor and wellbore storage is [25]: The dimensionless pressure solution of the PTA model can be obtained in the real-time domain by Stehfest inversion [23]. In the double logarithmic coordinate system, the dimensionless pressure derivative can zoom the characteristics of the dimensionless pressure features, and the dimensionless pressure derivative is:

Solution Validation
The dimensionless pressure solution of the PTA model in a two-layer reservoir with the same boundary radii is [26]: where Setting n = 2 and r e1D = r e2D = 1000, other parameters are listed in Table 1, the n-layer VNBR reservoir becomes a two-layer vertical uniform boundary radii (VUBR) reservoir. Figure 2 shows that the simplified VNBR PTA model matches that of the two-layer VUBR reservoir perfectly, which proves our model is reliable.

VNBR VUBR
Dimensionless wellbore storage (C D ) Setting n = 2 and re1D = re2D = 1000, other parameters are listed in Table 1, the n-layer VNBR reservoir becomes a two-layer vertical uniform boundary radii (VUBR) reservoir. Figure 2 shows that the simplified VNBR PTA model matches that of the two-layer VUBR reservoir perfectly, which proves our model is reliable.

Combined Type Curve
To characterize the effect of vertical non-uniform boundary radii (VNBR) on rate decline behavior and pressure response, we developed the combined type curves of VNBR and equivalent vertical uniform boundary radii (VUBR). The parameters of VNBR and equivalent VUBR for the combined type curves can be found in Table 2. The outer-convex VNBR, parabolic interpolation from 100 to 1000, is only a typical case of analysis. Based on the equivalent seepage volume (ESV) [21], the value of the equivalent VUBR can be calculated by the following equation: As shown in Figure 3a, the RTA combined type curves (including the normalized dimensionless decline rate curve (qDd), the normalized dimensionless decline rate integral curve (qDdi) and the normalized dimensionless decline rate integral derivative curve (qDdid)) show the PTA combined type curves of VNBR and VUBR are different, especially before the boundary-dominated flow. The RTA combined curves are much higher than the VUBR case even if the parameters are the same. The deployment of the dimensionless production integral derivative curve magnifies the differences so that we can capture the effect of VUBR. As shown in Figure 3b, the PTA combined type curves (including the dimensionless pressure curve (pwD) and dimensionless pressure derivative curve (p'wD)) show the flow of the VNBR reservoir can be divided into six flow regimes. Flow regime I is caused by wellbore storage. Flow regime II is induced by the skin effect near the wellbore. The pressure derivative is constant (L0 = 0.5) in flow regime III. In the early stage of flow regime IV, the pressure derivative begins to increase due to the pressure drop having reached the closest boundary. In the late stage of flow regime IV, the pressure derivative begins to decrease as the pressure drop propagates far through other layers. The pressure derivative curve appears a second constant value (L) during flow stage V. The pressure and its derivative increase and eventually coincide into the unit-slope line in flow regime VI.

Sensitivity Analysis
The β is introduced to represent the outer-convex relative thickness. In Case 1 (see Figure 4a), the outer-convex thickness is increasing but the maximum radius is constant. In Case 2 (see Figure  4b), the outer-convex number increases but each outer-convex thickness is the same. In Case 3 (see Figure 4c), the outer-convex positions are different but the each outer-convex is the same.
(a)  As shown in Figure 3a, the RTA combined type curves (including the normalized dimensionless decline rate curve (qDd), the normalized dimensionless decline rate integral curve (qDdi) and the normalized dimensionless decline rate integral derivative curve (qDdid)) show the PTA combined type curves of VNBR and VUBR are different, especially before the boundary-dominated flow. The RTA combined curves are much higher than the VUBR case even if the parameters are the same. The deployment of the dimensionless production integral derivative curve magnifies the differences so that we can capture the effect of VUBR. As shown in Figure 3b, the PTA combined type curves (including the dimensionless pressure curve (pwD) and dimensionless pressure derivative curve (p'wD)) show the flow of the VNBR reservoir can be divided into six flow regimes. Flow regime I is caused by wellbore storage. Flow regime II is induced by the skin effect near the wellbore. The pressure derivative is constant (L0 = 0.5) in flow regime III. In the early stage of flow regime IV, the pressure derivative begins to increase due to the pressure drop having reached the closest boundary. In the late stage of flow regime IV, the pressure derivative begins to decrease as the pressure drop propagates far through other layers. The pressure derivative curve appears a second constant value (L) during flow stage V. The pressure and its derivative increase and eventually coincide into the unit-slope line in flow regime VI.

Sensitivity Analysis
The β is introduced to represent the outer-convex relative thickness. In Case 1 (see Figure 4a), the outer-convex thickness is increasing but the maximum radius is constant. In Case 2 (see Figure  4b), the outer-convex number increases but each outer-convex thickness is the same. In Case 3 (see Figure 4c), the outer-convex positions are different but the each outer-convex is the same. (a)

Outer-Convex Thickness
As shown in Figure 5a, the RTA combined type curves of the VNBR move down and closer to that of the VUBR as the β increases. Therefore, the outer-convex thickness can be captured through the combined type curves. As shown in Figure 5b, L approaches L0 as the β increases, and the mathematical relationship is L = L0/β. This mathematical relationship clarifies the seepage physical meanings of β and provides a fast way to obtain the outer-convex thickness just from the pressure derivative curve.

Outer-Convex Thickness
As shown in Figure 5a, the RTA combined type curves of the VNBR move down and closer to that of the VUBR as the β increases. Therefore, the outer-convex thickness can be captured through the combined type curves. As shown in Figure 5b, L approaches L0 as the β increases, and the mathematical relationship is L = L0/β. This mathematical relationship clarifies the seepage physical meanings of β and provides a fast way to obtain the outer-convex thickness just from the pressure derivative curve.  Energies 2020, 13, 2305 7 of 13 curves move down due to the total thickness of outer-convex increasing from single to triple outerconvex as the outer-convex number is increasing. As shown in Figure 6b, the PTA combined type curves move down due to the outer-convex total thicknesses increasing. On the one hand, the start time of the VNBR transition flow regime is the same, as the minimum VNBR is unchanged. On the other hand, the start time of the VNBR boundary flow regime is later, since the equivalent VUBR becomes larger.

Outer-Convex Position
The outer-convex position refers to the outer-convex located at the top, middle and bottom of the reservoir (see Figure 4c). The ESV and β are unchanged due to the outer-convex thickness and outer-convex number being constant. Therefore, the rate decline curves (see Figure 7a) or pressure response curves (see Figure 7b) are completely coincident, regardless of the outer-convex position. The phenomenon indicates that the outer-convex position cannot be captured only from long-term production data or transient well test data.

Conclusions
This paper presents the extended RTA and PTA models of a vertical well in the VNBR reservoir by Laplace transform, Shehfest inverse and Duhamel's superposition principles, and further develops

Outer-Convex Position
The outer-convex position refers to the outer-convex located at the top, middle and bottom of the reservoir (see Figure 4c). The ESV and β are unchanged due to the outer-convex thickness and outer-convex number being constant. Therefore, the rate decline curves (see Figure 7a) or pressure response curves (see Figure 7b) are completely coincident, regardless of the outer-convex position. The phenomenon indicates that the outer-convex position cannot be captured only from long-term production data or transient well test data.

Conclusions
This paper presents the extended RTA and PTA models of a vertical well in the VNBR reservoir by Laplace transform, Shehfest inverse and Duhamel's superposition principles, and further develops

Combined Type Curve
To characterize the effect of vertical non-uniform boundary radii (VNBR) on rate decline behavior and pressure response, we developed the combined type curves of VNBR and equivalent vertical uniform boundary radii (VUBR). The parameters of VNBR and equivalent VUBR for the combined type curves can be found in Table 2. The outer-convex VNBR, parabolic interpolation from 100 to 1000, is only a typical case of analysis. Based on the equivalent seepage volume (ESV) [21], the value of the equivalent VUBR can be calculated by the following equation: As shown in Figure 3a, the RTA combined type curves (including the normalized dimensionless decline rate curve (qDd), the normalized dimensionless decline rate integral curve (qDdi) and the Energies 2020, 13, 2305 8 of 13 normalized dimensionless decline rate integral derivative curve (qDdid)) show the PTA combined type curves of VNBR and VUBR are different, especially before the boundary-dominated flow. The RTA combined curves are much higher than the VUBR case even if the parameters are the same. The deployment of the dimensionless production integral derivative curve magnifies the differences so that we can capture the effect of VUBR. As shown in Figure 3b, the PTA combined type curves (including the dimensionless pressure curve (pwD) and dimensionless pressure derivative curve (p'wD)) show the flow of the VNBR reservoir can be divided into six flow regimes. Flow regime I is caused by wellbore storage. Flow regime II is induced by the skin effect near the wellbore. The pressure derivative is constant (L 0 = 0.5) in flow regime III. In the early stage of flow regime IV, the pressure derivative begins to increase due to the pressure drop having reached the closest boundary. In the late stage of flow regime IV, the pressure derivative begins to decrease as the pressure drop propagates far through other layers. The pressure derivative curve appears a second constant value (L) during flow stage V. The pressure and its derivative increase and eventually coincide into the unit-slope line in flow regime VI.

Sensitivity Analysis
The β is introduced to represent the outer-convex relative thickness. In Case 1 (see Figure 4a), the outer-convex thickness is increasing but the maximum radius is constant. In Case 2 (see Figure 4b), the outer-convex number increases but each outer-convex thickness is the same. In Case 3 (see Figure 4c), the outer-convex positions are different but the each outer-convex is the same.

Outer-Convex Thickness
As shown in Figure 5a, the RTA combined type curves of the VNBR move down and closer to that of the VUBR as the β increases. Therefore, the outer-convex thickness can be captured through the combined type curves. As shown in Figure 5b, L approaches L 0 as the β increases, and the mathematical relationship is L = L 0 /β. This mathematical relationship clarifies the seepage physical meanings of β and provides a fast way to obtain the outer-convex thickness just from the pressure derivative curve.

Outer-Convex Number
As shown in Figure 4b, the β is 3/15 when a single outer-convex is located at the top, the β is 6/15 when the double outer-convex is distributed in the top and middle and the β is 9/15 when the triple outer-convex is distributed in the top, middle and bottom. As shown in Figure 6a, the RTA combined curves move down due to the total thickness of outer-convex increasing from single to triple outer-convex as the outer-convex number is increasing. As shown in Figure 6b, the PTA combined type curves move down due to the outer-convex total thicknesses increasing. On the one hand, the start time of the VNBR transition flow regime is the same, as the minimum VNBR is unchanged. On the other hand, the start time of the VNBR boundary flow regime is later, since the equivalent VUBR becomes larger.

Outer-Convex Position
The outer-convex position refers to the outer-convex located at the top, middle and bottom of the reservoir (see Figure 4c). The ESV and β are unchanged due to the outer-convex thickness and outer-convex number being constant. Therefore, the rate decline curves (see Figure 7a) or pressure response curves (see Figure 7b) are completely coincident, regardless of the outer-convex position. The phenomenon indicates that the outer-convex position cannot be captured only from long-term production data or transient well test data.

Conclusions
This paper presents the extended RTA and PTA models of a vertical well in the VNBR reservoir by Laplace transform, Shehfest inverse and Duhamel's superposition principles, and further develops the rate and pressure combined type curves to capture the production decline behavior and flow regimes' characteristics. Several conclusions and suggestions are obtained from this work.
(1) The production decline behavior of the VNBR reservoir is different from that of the VUBR reservoir, and the normalized dimensionless decline rate integral derivative can magnify the differences in the RTA combined type curves. (2) A new radial flow regime is identified in the VNBR reservoir that is captured by the PTA combined type curves. This behavior can be easily misinterpreted by a radial composite model as it has two stabilizations of the pressure derivative. Thus, the well testing analyst should be aware of the "equivalent composite zone effect". (3) The outer-convex thickness and outer-convex number of the VNBR can be quickly obtained from the second level of pressure derivative in the PTA combined type curve. The equivalent VUBR can be obtained from the boundary radius corresponding to the boundary flow regime based on the "equivalent seepage volume effect". (4) The actual exact outer-convex position of the VNBR, especially the VNBR with multiple outer-convexes, can neither be identified by RTA nor PTA combined type curves when the outer-convex total thickness and equivalent distance of the VNBR are unchanged.
In conclusion, the proposed RTA model enables petroleum engineers to better evaluate the production performance of a vertical well in the VNBR reservoir and interpret reservoir and well parameters more accurately based on the long-term production data. The PTA model can be used to analyze the pressure response of the VNBR reservoir, and quickly capture the flow regime characteristics from transient well test data. Acknowledgments: The achievements of this paper are supported by Sinopec's "Shitiaolong" special project (P18062). The corresponding author (CSC NO. 201906440196) would like to thank the China Scholarship Council for supporting this research at the University of Leeds, UK.

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

PTA
Pressure transient analysis, abbreviation RTA Rate transient analysis, abbreviation VNBR Vertical non-uniform boundary radii, abbreviation VUBR Vertical uniform boundary radii, abbreviation ∆ The thickness of the inner zone, m ∆p The pressure difference between initial formation and any position, Pa ∆p w The pressure difference between initial formation and wellbore, Pa Before production (t = 0), the pressure of any layer j is uniform and equals the initial formation pressure, which can be written as: At the wellbore position (r = r w ), the pressure of layer j is the same and equal to the bottom hole pressure. Therefore, the pressure at the inner boundary is: At the wellbore position (r = r w ), the wellbore rate is the summation of all the layers. Hence, the rate at the inner boundary is expressed as: There is no flow at the reservoir circle-closed boundary position (r = r ej ), which can be expressed as: To obtain the dimensionless mathematical model, first, define the following variables as shown in Table A1. Mobility thickness product ratios in layer j χ j = k j h jD /k a Storability ratios in layer j ω j = h jD (ϕc t ) j /(ϕc t ) a Based on the dimensionless parameters and dimensionless variables, the dimensionless equations of the RTA mathematical model in the real-time domain are:  where σ j = uω j /χ j .

Appendix B
In the Laplace domain, the pressure solution of layer j in Equations (A1) and (A2) can be written as [7,9]: p jD (r D , u) = A j I 0 (r D σ j ) + B j K 0 (r D σ j ) (A8) For any layer j, there are four unknown coefficients (A j , B j ) in Equation (A8). Since the range of j is from 1 to n, the number of unknown coefficients in the PTA model is 2n. Substituting Equation (8) and its derivative into the boundary conditions of Equation (7) yields: where the elements of Equation (A9) are: The A j and B j are obtained by solving Equation (9). Substituting A j and B j into Equation (8), the dimensionless pressure solution of any layer j can be obtained. Likewise, the dimensionless rate of layer j is obtained from Equation (A4) as follows: From Equation (A4) and Equation (A11), the dimensionless rate solution is obtained as: q D (u) = n j=1 q jD (u) (A12)