Approach to shape factor determination in estimating horizontal wells productivity for pressure transient model of Mutalik et al.

Horizontal well productivity (Jh) equation at pseudo-steady state (PSS) of Mutalik et al. (in: Presented at SPE 63rd annual technical conference and exhibition, Houston, Texas, 1988) has long and widely been known to estimate well performance. One key parameter of that equation is the so-called shape factor (CA) that characterizes the effect of the pressure transient behavior at PSS depending on well and reservoir configuration. For limited well and reservoir configurations, CA tables had been generated decades ago, but at unknown accuracy they are still used widely via tedious correlation or interpolation from those reported CA tables to estimate Jh for different configurations even today. The innovation points in this study are as follows: (1) The pressure transient model of Mutalik et al. (in: Presented at SPE 63rd annual technical conference and exhibition, Houston, Texas, 1988) (PTM) is improved to generalized pressure transient model (GPTM) and verified by converting several dimensional parameters to dimensionless forms. (2) An efficient algorithm is developed, and accurate CA and Jh are obtained for any well and reservoir configuration. (3) The accuracy of the reported CA tables is quantified, and the consequences are determined when they are used to estimate Jh. (4) New CA tables and their correspondence dimensionless time based on drainage area at PSS ((TDA)PSS) tables are generated for a wider range of well and reservoir configurations. Although horizontal wells usually have been drilled with high dimensionless well length (LD) industrially (i.e., LD > 10), it is still a significant error in Jh by (< 7%) that has been determined when using the reported CA of PTM for rectangular reservoir aspect ratio (xe/ye) = (1) and (2) in the both assumed input data, while for (xe/ye = 5), significant error in Jh starts from (46.43%) to (1037%) for LD = (10) to (100), respectively, in the first assumed input data [variable, dimensionless wellbore radius (rWD)] and expected the same in the second assumed input data (constant rWD).


Introduction
J h equation of Mutalik et al. (1988) is well known and presented by many authors (Duda et al. (1991), Economides et al. (1991b), Economides et al. (1991a), Frick and Economides (1993), Thomas et al. (1998), Penmatcha et al. (1999), Saavedra and Joshi (2002), Cho (2003), Jackson et al. (2011), Adesina et al. (2011), Ouyang (2015), Cetkovic et al. (2016), Wu et al. (2018), , and Zhang et al. (2020)). In reservoir simulators, J h equation of Mutalik et al. (1988) has been used to compute well productivity index (Krawchuk et al. 2006;Pinzon et al. 2007;Huebsch et al. 2008;Haidar et al. 2008;Brown and Tiwari 2010;Wood et al. 2010;and Leone et al. 2015). The J h equation is defined as follows: where J h is the horizontal well productivity (STB/day.psi), k is the reservoir permeability (md), h is the reservoir thickness (ft), µ o is the oil viscosity (cp), B o is the oil formation volume factor (bbl/STB), A is the drainage area (ft 2 ), r w is the well radius (ft), L is the length of horizontal well (ft), x e (1) (2) C A = 8.9832x e y e L 2 exp 2 2 T DA PSS − P D T PSS is the reservoir length in the X-direction (ft), the horizontal well is placed along the X-direction, y e is the reservoir length in the Y-direction (ft), and (P D T) PSS is the total dimensionless pressure at PSS. The J h equation [Eq. (1)] is for rectangular bounded reservoir at depletion stage (i.e., PSS). In C A equation [i.e., Eq. (2)], the (P D T) PSS was obtained by the PTM: where P D T is the total dimensionless pressure, T DA is the dimensionless time based on drainage area, x wD , y wD , and z wD are the dimensionless well location coordinates, x D is the X-coordinate for a point in the porous media. (x D =0) at the horizontal well center, and (x D = 1 and −1) at the well tips. yD is the Y-coordinate for a point in the porous media. (y D =r WD ) if the point is at the well wall, z D is the Z-coordinate for a point in the porous media. (z D = z WD ) if the point is at the well wall PTM [i.e., Eq. (3)] is based on 3D instantaneous source functions with Newman product method under homogenous and isotropic rectangular bounded reservoir at constant well flow rate and uniform flux wellbore condition of Gringarten and Ramey (1973). It is assumed that horizontal well is analogous to a vertical well of a vertical fracture. The equivalent infinite conductivity wellbore condition was obtained by solving P D T at x D equal to (0.732) (Mutalik et al. (1988) and Gringarten et al. (1974)). "Appendix 1" presents a summary of PTM.
When P D T slope with time becomes constant at PSS and is equal to (2π), both the (T DA ) PSS and (P D T) PSS have been obtained to use for C A calculation. C A at PSS also becomes constant for a given well and reservoir configuration.
PTM and their reported C A tables have identified three issues. Firstly, PTM as defined in Eq. (3) contains several dimensional parameters (i.e., x e , y e , h, and L) and this makes the solution not only awkward, but also difficult to use PTM to expand the reported C A tables to any other well and reservoir configurations. Secondly, PTM has not been verified in terms of P D T which may lead to errors in the reported C A tables. The third issue is the discrepancy in the assumed input data of the reported C A tables. The first assumption of the input data (3) exp −4n 2 2 x e y e h 2 T DA cos n z wD cos n z D dT DA (i.e., r w = 0.225 ft and y e = 933.33 ft) produces variable r WD where (r WD = 2r w /L) depends on penetration ratio (PR) and (x e /y e ). The second assumption of the input data shows that r wD is constant and is equal to (10 -4 ) for any well and reservoir configurations. These three limitations cast a doubt on the representativeness and usefulness of the reported C A tables for their respective configurations or other configurations for which simply interpolation or extrapolation of those tables is required as the correlative models previously presented by Bahadori (2012), Bahadori et al. (2013), and Ahmadi et al. (2015).
The main purposes of this study are in four parts.
(1) To overcome the first weakness in PTM (i.e., the dimensional parameters) by recasting PTM into a fully dimensionless form (i.e., GPTM) and to verify GPTM in terms of P D T. (2) To develop efficient algorithm to compute C A and J h accurately for any well and reservoir configuration. Because C A depends on the all-configuration parameters, unlimited tables of C A can be generated; therefore, the MATLAB code is provided. (3) To quantify the accuracy of the reported C A tables in their both assumed input data (i.e., constant and variable r WD ), and to determine the consequences of any identified errors in the reported C A tables when they are used to estimate J h . (4) To generate new C A tables for a wider spectrum of well and reservoir configurations.

The generalized pressure transient model (GPTM)
The PTM is improved to fully dimensionless form (i.e., to GPTM) by applying the following dimensionless parameters into Eq. (3): where x eD is the dimensionless reservoir length in the X-direction and y eD is the dimensionless reservoir length in the Y-direction For anisotropic reservoir solution, the dimensionless parameters are presented in "Appendix 2." The following meaningful mathematical relationship should be noted in the well and reservoir configuration: Used for centrally and off-centrally well loca- Used for centrally and OCWL in the X-direction. z wD = z w h Used for centrally and OCWL in the Z-direction. (6) y eD = 2y e L (7) exp −4n 2 2 x eD y eD L 2 D T DA cos n z wD cos n z D dT DA Note that OCWL is used to indicate the offset of the well relative to the central well location (i.e., OCWL = 1) at a respective coordinate direction for the purpose of classifying the configurations only.
The efficient algorithm to solve Eq. (7) is developed and implemented in MATLAB code and is given in "Appendix 3." P D T by GPTM of this study [i.e., Eq. (7)] is verified against its counterpart by the horizontal well pressure transient model of  for a square bounded reservoir (i.e., x e /y e = 1). L D is equal to (10), and PR = (0.2) and (1), for a well that is centered in x-, y-, and z-directions [i.e., (x wD /x eD = 1), (y wD /y eD = 1), and (z wD = 0.5), or (OCWL = 1)]. The wellbore condition is uniform flux (x D = 0), and r wD is equal to (2 × 10 -4 ). Figure 1 shows the identical results of P D T by both the models.
For further verification on GPTM, P D T of this study is compared with P D T of the pressure transient model of Gringarten et al. (1974) for vertical fracture well under square bounded reservoir. (L D = 10 3 ) (i.e., to make the thickness of the fracture to be zero), PR = (0.5) and (1), (OCWL = 1),

Fig. 1
Comparison of P D T by GPTM and by the pressure transient model of  for a centrally well located in a square reservoir (x D = 0), and (r wD = 10 -4 ). Figure 2 shows the identical results of P D T by both the models.
The C A equation [i.e., Eq.
In what follows, the reported C A tables by PTM for selected well and reservoir configurations are examined, and their well and reservoir configurations can be classified into four cases in terms of OCWL. Table 2 lists the four cases of the reported C A tables. The parameters in Cases (1) to (3) correspond to tables (4), (7), and (8) in Mutalik et al. (1988) paper, respectively. Note that parameters in Case (4) are deduced from Fig. 12 in Mutalik et al. (1988) paper.

The accuracy of the reported C A tables by PTM (case 1)
Although the accuracy of the reported C A tables by PTM for other cases has been examined, here only the results for the Case 1 are presented for both assumed input data (i.e., the discrepancy issue). For each C A by PTM, the counterpart C A Fig. 2 Comparison of P D T by GPTM and by the pressure transient model of Gringarten et al. (1974) for vertical fracture well under square reservoir Table 1 Comparison of C A by this study, Earlougher (1977) and Hagoort (2009) Earlougher (1977) FTLC of this study Hagoort (2009) FTLC of this study 10 -5 10 -9 10 -5 10 -9 for the both assumed input data. Tables 3 and 4 show the relative errors in the reported C A tables by PTM for the variable r wD and the constant r wD , respectively. For the first assumed input data (variable r wD ), Table 3 shows that in x e /y e = (1) and (2), the C A errors of any PR are in reverse order with L D , while in x e /y e = (5), they are in non-order with L D . The maximum C A error is (50%) at (PR = 1), (L D = 1), and (x e /y e = 2). For the second assumed input data (constant r wD = 10 -4 ), Table 4 shows that in x e /y e = (1) and (2), the C A errors of any PR are in reverse order with L D , while in x e /y e = (5), they are in non-order with L D . The maximum C A error is (1215%) at (PR = 0.2), (L D = 1), and (x e /y e = 1). For any assumed input data, (T DA ) PSS by GPTM is noticeably close to its counterpart by PTM and both are independent of L D .
The errors in the reported C A tables are also found in the cases of off-center well location in x-and y-directions for both assumed input data but are not reported here.
To determine the consequences of the identified errors in the reported C A tables of Mutalik et al. (1988) when they are used to estimate J h , an example problem is applied. The reservoir of (y e = 933.33 ft) and (r w = 0.225 ft) is for the first assumed input data, and the reservoir of (y e = 2640 ft) and (r wD = 10 -4 ) is for the second assumed input data, and all are at (k = 1 md), (µ o = 0.5 cp), and (B o = 1.2 rbbl/STB). The J h is calculated by using Eq. (1), and the C A (s) by PTM and GPTM is applied by using Tables 3 and 4. For the first assumed input data, PR = (0.4), (1), and (0.8) for x e /y e = (1), (2), and (5), respectively. For the second assumed input data, PR = (0.2) and (0.8) for (x e /y e = 1), and for x e /y e = (2) and (5) they are not reported here. (OCWL = 1) and L D is from (1) to (100). Figures 3 and 4 show a comparison of those J h (s) and C A (s) for the both assumed input data, respectively. Figure 3 shows that the reported C A tables of PTM lead to underestimation in the J h for x e /y e = (1) and (2) and lead to overestimation in the J h for x e /y e = (5). The relative error in the J h for each selected x e /y e = (1), (2), and (5) follows the same trend as those in C A table. The maximum underestimation in J h is (8.29%) corresponding to the error in C A at (x e /y e = 1) and (L D = 1), while the maximum overestimation in J h is (1037%) corresponding to the error in C A at (x e /y e = 5) and (L D = 100) . For both x e /y e = (1) and (2), although the maximum error in J h (i.e., underestimation) at (L D = 1), a rather unrealistic case in most real fields, from Fig. 3 the error in J h can still reach (6.38%) underestimation for even a realistic configuration of L D equal to (10). For (x e /y e = 5), the error in J h (i.e., overestimation) is much more significant and starts from (46.43%) to (1037%) overestimation for L D = (10) to (100), respectively. Figure 4 shows that the reported C A tables of PTM lead to overestimation in the J h for both PR = (0.2) and (0.8). The relative error in the J h for each PR follows the same trend as those in C A table. The maximum overestimation in J h is (32%) corresponding to the error in C A at (L D = 1) and (PR = 0.2). Although the maximum error in J h (i.e., overestimation) is at (L D = 1), a rather unrealistic case in most real fields, from Fig. 4 the error in J h can still reach (7%) overestimation for even a realistic configuration of L D equal to (10).

The source of errors in the reported C A tables
Given that both PTM and GPTM solutions are derived identically, except re-casting some variables into dimensionless forms in the latter model, the source of the errors in the reported C A tables may be originated from the evaluation of the P D T which was not verified by Mutalik et al. (1988). Figure 5 shows the effect of using FTLC = (0-10 -9 ) and (0-10 -5 ) on the C A accuracy. GPTM is used to calculate C A for square reservoir, (r WD = 10 -4 ), (OCWL = 1), PR = (0.1), (0.5), and (1), and L D from (1) to (10 3 ). The maximum error in C A is (3630%) which appears at the lowest PR (i.e., PR = 0.1) and L D (i.e., L D = 1). Notice that the error in C A becomes small and not significant at (L D > 100) and suggests that the FTLC of (0-10 -5 ) would be sufficient for C A accuracy when (L D > 100) for any well and reservoir configuration. This has also been shown previously in Table 1 for vertical fracture well case when (L D > 10 3 ).
The pressure drop that occurs at the very early time in the P D T solution, especially when L D is small (i.e., when the reservoir thickness is large), has significant effects on the C A accuracy. The evaluation of P D T in that period requires  to use a much shorter FTLC [i.e., FTLC = (0-10 -9 )] to gain sufficient accuracy. Due to the absence of the information about the selected FTLC in the calculation of the reported C A tables by PTM, C A tables have errors in reverse order with L D as shown previously for x e /y e = (1) and (2) in Tables 3 and 4. Therefore, the source of errors in the reported C A tables by PTM could be originated and justified from the using of long FTLC in the P D T evaluation. About (x e /y e = 5), no analysis can be achieved. Figure 6 shows the effect of OCWL in the Z-direction on C A for square reservoir, (r WD = 10 -4 ) and PR is equal to (0.1) and (0.9). The effect of z wD (i.e., z D = z wD ) on the C A that is calculated by GPTM becomes insignificant when (L D > 35) (C A % error < %5) which is inconsistent with their C A counterparts by PTM [i.e., Fig. 12 in Mutalik et al. (1988) paper] because z wD will be insignificant when (L D > 5). Therefore, z wD should be considered as a major parameter for any correlation. Table 5 which is presented at the end of the paper shows the effect of the wellbore radius on C A for square reservoir and PR = (0.1) and (1). When r wD increases two times (i.e., from (5 × 10 -5 ) to (10 -4 )), a significant error occurs when (L D < 15), and when r wD increases (10) to (20) times [i.e., from (5 × 10 -4 ) to (10 -3 )], a significant error occurs when (L D < 50). This is also inconsistent with their C A counterparts by PTM (i.e., table (9) in Mutalik et al. (1988) paper), the significant error will occur when (L D < 3) if r wD increases 10 times. Therefore, r wD should be considered as a major parameter for any correlation.  Table 4 Comparison between C A and (T DA ) PSS of this study and Mutalik et al. (1988) (when constant r wD = 10 -4 ) PR = 0.2, Central well, x eD /y eD = 1, r wD = 10 -4 PR = 0.4, Central well, x eD /y eD = 1, r wD = 10 -4 L D Mutalik et al. (1988) This study C A % Err L D Mutalik et al. (1988) This study C A % Err PR = 0.6, Central well, x eD /y eD = 1, r wD = 10 -4 PR = 0.8, Central well, x eD /y eD = 1, r wD = 10 -4 L D Mutalik et al. (1988) This study C A % Err L D Mutalik et al. (1988) This study PR = 1, Central well, x eD /y eD = 1, r wD = 10 -4 PR = 0.2, Central well, x eD /y eD = 2, r wD = 10 -4 L D Mutalik et al. (1988) This study C A % Err L D Mutalik et al. (1988) This study PR = 0.4, Central well, x eD /y eD = 2, r wD = 10 -4 PR = 0.6, Central well, x eD /y eD = 2, r wD = 10 -4 L D Mutalik et al. (1988) This study C A % Err L D Mutalik et al. (1988) This study C A % Err  PR = 0.2, Central well, x eD /y eD = 5, r wD = 10 -4 PR = 0.4, Central well, x eD /y eD = 5, r wD = 10 -4 L D Mutalik et al. (1988) This study C A % Err L D Mutalik et al. (1988) This study C A % Err PR = 0.6, Central well, x eD /y eD = 5, r wD = 10 -4 PR = 0.8, Central well, x eD /y eD = 5, r wD = 10 -4 L D Mutalik et al. (1988) This study C A % Err L D Mutalik et al. (1988) This study

Conclusions
1. The PTM was improved to GPTM and verified by converting several dimensional parameters to dimensionless forms. 2. An efficient algorithm was developed, and accurate C A and J h were obtained for any well and reservoir configuration. 3. The accuracy of the reported C A tables was quantified, and the consequences were determined when they were used to estimate J h .
4. New C A tables and their (T DA ) PSS tables were generated for a wider range of well and reservoir configurations. The tables are provided in the Supplementary Material. 5. Because C A depends on the all-configuration parameters, unlimited tables of C A can be generated; therefore, the MATLAB code was provided. 6. Although horizontal wells usually have been drilled with high L D industrially (i.e., L D > 10), it is still a significant error in J h by (< 7%) that has been determined when using the reported C A of PTM for (x e /y e ) = (1) and (2) in the both assumed input data, while for (x e /y e = 5), a significant error in J h starts from (46.43%) to (1037%) for L D = (10) to (100), respectively, in the first assumed input data (i.e., variable r WD ) and expected the same in the second assumed input data (i.e., constant r WD ).

Appendix 1
Summary of pressure transient model of Mutalik et al. (1988) The mathematical expression of the dimensionless parameters of PTM (Eq. (3)) and the geometrical figure of PTM is given as follows: 12) x D = 2x L Fig. 6 The effect of OCWL in the Z-direction on C A

Appendix 2 The dimensionless parameters in anisotropic reservoir
The mathematical expression of the dimensionless parameters for anisotropic reservoir by Ozkan et al. (1987),  and Besson (1990) in (') symbol is given as follows: (13) y D = 2y L = r wD = 2r w L
Funding No funding supported this study.

Conflict of interest
The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.