Performance forecasting of hydraulically fractured horizontal wells using dynamic-drainage-area concept with correct distance of investigation

Performance forecasting of multi-fractured horizontal wells (MFHWs) completed in tight/shale oil or gas reservoirs is of great significance in the development of unconventional resources. Dynamic drainage area (DDA) concept has emerged as a production forecasting methodology for unconventional reservoirs. This paper investigated the DDA method and derived correct expressions of distance of investigation (DOI) for both constant production rate and constant bottom-hole pressure cases by integrating material balance and deliverability equation within the range of DOI. The modified DDA method with correct DOI coefficients provided in this paper permits direct calculation of dynamic performance at arbitrary time steps before the end of transient linear flow. At the same time smooth production forecasting of MFHWs from transient linear to boundary-dominated flow is realized without modification by extra coefficient. Hybrid models such as DDA plus dual-exponential and DDA plus hyperbolic are presented, which can be applied quickly and easily to MFHWs in unconventional reservoirs as alternatives to complex numerical simulation. Meanwhile, average pressure in the range of DDA can be readily obtained with correct DOI coefficients, avoiding complex iterative calculations. The reliability and practicability of this solution have been demonstrated by synthetic and field cases in this work.

Distance of investigation, ft y invD Dimensionless distance of investigation, y invD = y inv ∕x f

Introduction
Over the past two decades of years, multi-fractured horizontal wells have become the primary and fundamental recovery technology for tight/shale oil and gas reservoirs. Besides complex numerical simulation, different analytical, semianalytical, and empirical methods (Clarkson 2013(Clarkson , 2014 have been developed for history matching and production forecasting of MFHWs completed in these unconventional resources. Due to the ultra-low permeability associated with unconventional reservoirs, the major flow regime usually exhibits a very long transient linear flow with or without boundary-dominated flow depending on the property of the reservoir and characteristics of fractures. The rigorous reservoir numerical simulation technique usually cannot be applied to every well due to a lack of supporting data and time-cost problems. Empirical methods, such as power-law exponential (PLE) (Ilk et al. 2008), stretched exponential production diagnostic (SEPD) (Valko and Lee 2010), and the Duong model (Duong Anh 2014), even though remain popular, do not rigorously account for flow regime and fracture geometries. Shahamat et al. (2014) introduced a concept of continuous succession of pseudo-steady states (SPSS), in which the flow regime in the region of investigation is assumed to be pseudo-steady state. He used this method to predict reservoir performance by combining material balance equation, boundary-dominated flow and distance of investigation. A capacitance-resistance methodology (CRM) was presented, in which the capacitance is expressed as the product of the total system compressibility and the reservoir volume of DOI and resistance is equivalent to the inverse of productivity index in deliverability equation.  extended the idea of SPSS and introduced dynamic drainage area (DDA) concept, which is actually a continuous succession of steady state. As a production-forecasting method, the DDA concept has attracted much more attention because of its simplicity and physicalmodel-based nature of application. The DDA method has been extended to cases of multi-wells and multi-phases with complex fracture geometry (Clarkson andQanbari 2015, 2016a, b;Qanbari and Clarkson 2016;Shahamat and Clarkson 2018;Ahmadi et al. 2021). The DDA method combines the deliverability equation, material balance equation and distance of investigation together and requires an iterative calculation procedure to determine the average pressure.
The calculation of DOI is critical in determining the end of linear flow, as well as the flow continuity in the linear-toboundary (LTB) model. Some "hybrid" approach has been used for modeling the LTB system (Nobakht and Clarkson 2011;Nobakht andClarkson 2011, 2012), in which different models are adopted separately for transient-linear-flow and boundary-dominated flow periods, such as square-root-oftime plot and hyperbolic decline model. The key problem here is how to determine the end of linear flow ( t elf ).
In this work, we modified and developed the DDA method by introducing a derived DOI equation, avoiding the iterative calculation and maintaining a good continuity between transient linear and boundary-dominated flow. Average pressure in the range of DOI can be directly calculated without iteration as previously required. This new method is robust and flexible enough for performance forecasting of MFHWs under various conditions. The new modified DDA approach has been validated using synthetic simulation and field examples, showing reasonable match both in transient and boundary dominated flow periods.
In the following, the detailed development of the new method is presented for constant production rate and constant bottom-hole pressure cases. Applications on performance forecasting of MFHWs and comparisons with numerical simulation case are also provided to demonstrate the practical applicability of this solution.

Theory and model development
The calculation of distance of investigation is essentially required in both SPSS and DDA methods. The following equation is generally used. (1) The coefficient α must be carefully selected, which covers a large range depending on the criterion used for defining the DOI. Many values of α have been suggested based on the type of production (constant-rate or constant-pressure) or selective criteria (gauge resolution, rate ratio and cut-off values). Both Shahamat et al. (2014) and ,  adopted the empirical values suggested by Wattenbarger et al. (1998), which are α = 0.113 for the case of constant production rate and α = 0.159 for constant production pressure. When expressed in dimensionless form, Eq. (1) becomes The constant 0.159 of α is equivalent to D = 2 in Eq. (2). When DOI reaches the boundaries of formation or fractures inferences start, the value of y inv equals to y e , marking the start of pseudo steady state (PSS) and the time equals to the start of boundary-dominated flow ( t BDF ). In the iterative calculation of DDA method, an estimated value of t BDF must be provided in advance based on DOI equation. There is no supporting proof provided by Nobakht and Clarkson (2012), Shahamat et al. (2014) and ,  for the selection of 0.159. On the other hand, because of the inappropriate value of coefficient D in DOI equation as Eq.
(2), an extra coefficient α o was forcibly incorporated to keep continuity of production rate at the end of transient linear flow and start of boundary-dominated flow (Clarkson and Qanbari 2015).
In this work, we reinvestigated and simplified the DDA method using a derived value of coefficient D, which can avoid the complex iterative algorithm for average pressure and remove the artificial continuity coefficient α o . We demonstrated that in DDA concept the 0.159 value is just suitable for constant production rate transient radial flow case (Appendix B). The correct D values for transient linear flow at different production-control conditions are derived in the following sections.

Basic model
The basic physical model of MFHW is depicted in Fig. 1. For simplicity we just select one element of MFHW as shown in Fig. 1b, similar to that of Wattenbarger et al. (1998) and Nobakht and Clarkson (2012), in which a hydraulic fracture lies in the center traversing the reservoir. For the model in Fig. 1b, the reservoir exhibits transient linear flow until the investigated distance equals to the element length in y-direction and thereafter boundary-dominated flow prevails (and fractures interference starts).
In the next two sections, we derived the expressions of DOI for DDA method at constant production rate and constant bottom-hole pressure cases. Even though the performance study of MFHWs involves transient linear flow only, we still incorporated the derivation of D-value in transient radial flow (which is shown in Appendix B). This is for theoretical completeness and to prove that the D-value of 2 (and the corresponding α = 0.159) is only suitable for transient radial flow. The DOI equation for transient radial flow cannot be used for transient linear flow.

Constant production rate case
The key algorithm of DDA is to find an average pressure to calculate oil and gas rate using an iterative algorithm. Assuming p y inv , t is the average pressure in the drainage area within the distance of investigation. The following expression holds unconditionally.
The right side of Eq. (3) implies that the total pressure drop of the well can be separated into two terms. The first pressure drop from initial pressure to average pressure accounts for reservoir depletion, described by material balance equation. The second pressure drop is due to the flow from reservoir into the well. For constant-production-rate case in transient linear flow, we can give the following material balance equation for the model described in Fig. 1b.
Define a dimensionless material-balance pressure at constant-rate production, Rearranging Eq. (4) and substituting Eq. (2) yields As previously mentioned, the DDA method assumes a succession of steady state flow. The productivity equation (suggested by Clarkson and Qanbari 2015, and revised here) is Based on Eq. (7), we define a steady-state dimensionless pressure for constant production rate case as implied in DDA concept.
Multiplying kh∕ (141.2q B) on both sides of Eq. (3), we have From solution of transient linear diffusivity equation (Wattenbarger et al. 1998), we have Substituting Eq. (6), Eq. (8), and Eq. (10) into Eq. (9), we obtain Solving Eq. (11) yields the value of D, which is 1.651. This value of D is corresponding to α value of 0.131 in (1). This value of α is greater than 0.113 from unit impulse method and less than 0.141 from method of intersection suggested by Behmanesh et al. (2015). This derived D value implies that the average pressures in material balance and deliverability equation are equivalent to each other. Most importantly, this D value avoids the iterative calculation of average pressure suggested by Qanbari 2015 2016a, b;Qanbari and Clarkson 2016

Constant production pressure case
In constant-production-pressure situation, the production rate varies with time. Dividing (p i − p wf ) on both sides of Eq. (3) gives Like constant-production-rate case, the first term on the right side of Eq. (12) is material balance contribution due to reservoir pressure depletion. The second term can be treated as dimensionless pressure drop caused by fluid flow toward the hydraulic fracture.
Because production rate varies with time in this case, we modify Eq. (4) as Define dimensionless cumulative production and dimensionless material-balance-pressure at constant-production-pressure case, as Eq. 14 and Eq. 15 From Eq. (13), we have The production rate in transient linear flow is given by We note that dimensionless cumulative production can be alternatively derived from the integration of dimensionless production rate.

Equation (12) becomes
From the definition of dimensionless production rate and the deliverability equation of Eq. (7) (1), which is greater than 0.194 from unit impulse method and 0.180 from method of intersection suggested by Behmanesh et al. (2015). This value of D for constantproduction-pressure case also implies that we can avoid the iterative algorithm to calculate the average pressure in material balance and deliverability equation. Meanwhile, this derived value of D makes it easier to estimate the end of linear flow or the start of boundary-dominated flow based on equation of DOI. The match of dimensionless production rate under different D values is illustrated in Fig. 3. It is noted that dimensionless production rate q D under D = 2.84 is in excellent agreement with the theoretical solution. Lower value of D than 2.84 will overestimate production rate, while higher D value than 2.84 will underestimate production rate. Both flow regime and production condition can affect the selection of correct DOI coefficient. For constant-production-pressure case, different values of D are corresponding to different dimensionless front-bottomhole pressure ratio (FBPR, which define the pressure front for constant BHP case), p Dr . According to the definition of FBPR (Fan 2021), the D value of 2.84 is related to a FBPR of 0.045. Comparison of p Dr for different D values is listed in Table 2.

Average pressure calculation
Average reservoir pressure (p av ) is one of the fundamental parameters in reservoir engineering calculations. It is essential to evaluate fluid properties in material balance applications and to calculate pseudo-time for pressure-sensitive reservoirs. Traditionally the direct method of estimating p av is based on the relation of flowing bottom-hole pressure and average reservoir pressure during pseudo-steady-state (PSS) of a constant production rate system. However, in unconventional reservoirs, it usually takes very long time to reach PSS flow and it is almost impossible to reach PSS for constant bottom-hole pressure case. Anderson and Mattar (2007) proposed that during transient flow the pseudo-time should be evaluated at the average pressure of the region of influence rather than the average reservoir pressure. For DDA and SPSS methods, the average pressure within the area of distance of investigation (DOI) is an essential parameter that can only be solved by iterative methods. Alternatively, the average pressure can also be evaluated by the volumetric average method, which is highly dependent on the definition of distance of investigation. For transient linear flow of MFHWs, the dimensionless average pressure in the influence region can be expressed by From the transient linear flow solution of constant-production-rate case (Wattenbarger et al. 1998;Behmanesh et al. 2015), it has Similarly, from the solution of transient linear flow at constant production pressure case (Wattenbarger et al. 1998;Behmanesh et al. 2015), the dimensionless average pressure in the region of influence for constant BHP is Based on the derived value of D = 2.84 for constant production pressure, it gives From the definition of dimensionless pressure in constant BHP case, we can obtain another direct calculation method of average pressure from initial pressure and bottom-hole pressure in dimensional form.
We note that the average pressure for constant production rate is a function of square-root-of-time, while the average pressure for constant production pressure case is a constant value, as shown in Fig. 4 and Fig. 5. Coefficients of these two relations depend on the selection of DOI equation. Equation (28) and Eq. (31) indicate that average pressure in the drainage area is a unique function of initial pressure and flowing wellbore pressure, which allows for a direct calculation of pseudo-time and avoids an iterative process in current practice. Roadifer (2015) derived different average pressure expressions for the above two cases, which are p = 0.6 √ t D and p = 0.4411 , respectively. The reason for the differences is that he used the DOI equations from Behmanesh (2015), which (25) are D = √ 2 for constant production rate and D= √ 6 for constant production pressure case. In fact, these two values are not suitable for DDA concepts. For instance, if we use D= √ 6 for constant pressure case, we can derive p mbD,p = 0.23 and p ssD,p = 0.69 , which makes the equation p mbD,p + p ssD,p = 1 not hold. The derived coefficients in this work actually synchronize the change of average pressure in material balance and deliverability equations within the dynamic drainage area, leading to a direct calculation of dynamic performance at arbitrary time step before the end of transient linear flow.
After the end of linear flow, we can derive the expression of average pressure from material balance equation (Roadifer and Kalaei 2015).
In which, t elf is the time at end of linear flow, p elf is the average pressure at end of linear flow and N pelf is the cumulative production at end of linear flow. Wattenbarger et al. (1998) presented the analytical solution to hydraulically fractured wells under constant bottom-hole condition, as in Eq. (34). However, the use of such an infinite-series function is not necessarily desirable for production data analysis. Using DDA method, we can simplify the production forecasting procedure. As previously described, under correct DOI the iterative calculation of average pressure can be avoided. We can determine the time at the end of linear flow or start of boundarydominated flow under correct value of D coefficient. This is crucial for selection of predicting equations.

Production forecasting
Substituting Eq. (18) and Eq. (2) into material balance Eq. (13) yields From deliverability equation Eq. (7) and the definition of dimensionless production rate, we can obtain From the summation of the above two equations, we can derive an expression of dimensionless production rate, as in Eq. (37).
This equation is equivalent to the linear approximation of Eq. (34), which can be used to forecast production before the boundary is reached or the fractures interference starts. The time at end of linear flow or the start of boundary-dominated flow is easy to derive from Eq. (2).
When fracture half-length is not known explicitly, the slope of the inverse of production rate vs. square-root-oftime plot m L can be used to estimate its value. Fig. 4 Pressure profile and average pressure at constant-rate case After the end of linear flow, the long-term exponential approximation is generally suggested for pressure or flow rate dynamic analysis. However, we adopt a dual-exponential model here, which is derived by taking the first two terms of the infinite-summation series in Wattenbarger's solution.
The advantage of this dual-exponential model is that a smooth continuity of flow rate is obtained, without introducing a connecting coefficient suggested by Clarkson. This is illustrated in Fig. 6. The blue circle series (DDA + EXP D = 2) is from DDA method suggested by Clarkson (Clarkson and Qanbari 2015, 2016a, b, Qanbari and Clarkson 2016Ahmadi et al. 2021), which is not continuous between transient linear flow and boundarydominated flow without introducing additional connecting coefficient. The reason of this discontinuity is due to inappropriate D-value in DOI equation.
It is obvious that late-time data will eventually deviate from early straight line and marks the end of transient linear flow. In many cases, however, late-time data does not follow a strict theoretical exponential decline. The reasons may include, properties change due to pressure drop, contribution of unstimulated region, heterogeneities of fracture and reservoir properties, or transitional flow regimes (elliptical/radial) before a complete boundarydominated flow regime is reached. From reservoir practices, hyperbolic model (Nobakht and Clarkson 2012) or SEPD model (Valko and Lee 2010) are suggested after the end of linear flow. (39) The hyperbolic model is (Nobahkt et al. 2010) In which, q Delf and D elf are the production rate and decline rate at the end of linear flow, respectively, which can be derived based on Eq. (37). The constant b is the decline exponent, which varies from 0 to 1.
The SEPD model is (Valko and Lee 2010) In which, n and D are constants, which can be determined from group average or statistical experience.

Synthetic case
To demonstrate the applicability of the modified DDA method in this paper, we generated synthetic production data using reservoir simulation software. The model is an element (as shown in Fig. 1b) of hydraulically fractured well from a tight oil reservoir. The bottom-hole pressure is set to ensure pure single-phase linear flow and the producing time is set long enough to generate both transient and boundarydominated flow. The parameters are listed in Table 3.
The reservoir starts with a transient linear flow and the pressure distribution after two months is shown in Fig. 7. After outer boundary is reached, the production does not follow a strict exponential decline, as illustrated in Fig. 8. Before the end of transient linear flow, DDA method is used to match the production rate. After the end of linear flow, both hyperbolic model and SEPD can match the simulation results. It shows that the SEPD model can provide a much better match than hyperbolic decline.

Field case
The field case is a multi-fractured horizontal well completed in DL tight oil reservoir of Shengli oilfield, East China. The number of hydraulic fractures is 10 and the average fracture spacing is 300ft. The oil production is shown in Fig. 9. From the log-log plot of production rate versus time, a transient linear flow period is observed. The late-time data exhibits a deviation from early straight line, indicating fracture interference and boundary effect. The fracture half-length is 256ft, determined from the slope of the early straight line. The modified DDA method is used to match the early transient linear flow period. After the end of linear flow, a hyperbolic decline is adopted to match the production decline, in which D elf is 1.79 and the value of b equals to 1.

Conclusions
1. This paper investigated and developed dynamic drainage area (DDA) method for performance forecasting of multi-fractured horizontal wells (MFHWs) completed in unconventional reservoirs. Correct expressions of distance of investigation (DOI) for DDA method are derived for both constant production rate and constant bottom-hole pressure cases. The DOI coefficients of D-value are 1.651 and 2.84, respectively. 2. Simple and direct expressions of average pressure for performance forecasting of MFHWs are presented to avoid complex iterative calculation in DDA method. 3. Hybrid models, such as DDA plus dual-exponential, DDA plus stretched exponential production diagnostic model, and DDA plus hyperbolic model are presented for performance forecasting of MFHWs in unconventional reservoirs as alternatives to complex numerical simulation, bridging the gap between analytical and empirical methods. 4. This modified DDA method can be extended to MFHWs in unconventional gas reservoirs by adopting pseudopressure. Reliable forecasts can be obtained without using pseudo-time, reducing complexities and iterative procedures.  Production match of field case based on modified DDA method 5. Even though the results in this work are aimed at performance forecasting of MFHWs, the introduced technique is readily adaptable for transient radial flow and other geometries. This study will be of interest to those petroleum engineers who desire to find a simple tool for performance forecasting of MFHWs in unconventional reservoirs.