A New Decline Curve Analysis Method for Layered Reservoirs

A New Decline-Curve-Analysis Method for Layered Reservoirs Kittiphong Jongkittinarukorn, Chulalongkorn University; Nick Last, Well Test Knowledge International; Freddy Humberto Escobar, Universidad Surcolombiana; and Kreangkrai Maneeintr, Chulalongkorn University Summary This study presents a new method to improve production forecasts and reserve estimation for a multilayer well in the early stages of production using the Arps (1945) hyperbolic decline method to model the decline rate of each layer. The method can be applied to both oil and gas wells. The new approach generates the profiles of the instantaneous decline rate (D) and instantaneous decline-curve exponent (b) from the historical flow rate (q). Because of the inherent noise in the production data, a regression technique is applied to smooth the flow-rate data, and the analysis is performed on the smoothed data. History matching is performed not only on the profile of q but also on the profiles of D and b. This results in the unique decline parameters (qi, Di, and b) for each layer. For a multilayer well, the values of D and b vary with time, which means that the well’s performance cannot be modeled using a conventional single-layer-well approach. Furthermore, the well-known nonuniqueness problem from history matching is magnified in a multilayer well: Many models can successfully match the production profile in the short-term but fail to match it in the longer term. Only the correct model can match the profiles of q, D, and b over both the short-term and the long-term. The proposed method provides the correct unique decline parameters (qi, Di, and b) for each layer, during the early stages of production, and these parameters are then valid for the life of the well. The method works well for both synthetic examples and actual field data. The novelty of the new methodology is the ability to provide the decline parameters for each layer at an early stages of production that can then be used for production forecasting in the long-term. The nonuniqueness problem from history matching is solved. Introduction Decline curve analysis (DCA) is a commonly applied petroleum engineering method used for forecasting production and estimating reserve from production data under boundary-dominated flow. For a volumetric, single-phase reservoir produced at constant bottomhole pressure, the flow rate can be empirically modeled by the Arps (1945) hyperbolic function. Note that the value of b is constant. High values of b can be observed during the infinite-acting or transient-production period (Fetkovich et al. 1990; Cheng et al. 2008), and this part of the well’s history should not be used for long-term performance predictions. During boundary-dominated flow, a single-layer well has a value of b < 0.5 (Fetkovich et al. 1996; Okuszko et al. 2007). Higher values of b (b > 0.5) can be observed for a multilayer well without crossflow (Fetkovich et al. 1996). This is caused by an early rapid-rate decline of the higher-permeability layer followed by an extended period of low-rate decline of the lower-permeability layer. Normally, a multilayer well has an unusually long production life. Arps (1945) developed the hyperbolic-rate decline empirically. However, the values of D and b could be related to fluid properties and production conditions (Fetkovich 1980; Fetkovich et al. 1996; Chen and Teufel 2002; Ayala and Ye 2013; Ye and Ayala 2013; Stumpf and Ayala 2016). Fetkovich et al. (1996) showed that the initial decline rate (Di) and b can be expressed in reservoirengineering terms. Chen and Teufel (2002) proposed a model to estimate an averaged b during boundary-dominated gas flow using pseudopressure. Fetkovich (1980) used the departure-curve method (differencing) to analyze decline-curve data from the East Side Coalinga Field, which consists of upper and lower oil sands separated by shale. At late time, the production was dominated by Layer 1. History matching during this period yielded the decline parameters (qi1 and N1) for Layer 1. The production rate for Layer 1 was extrapolated backward in time, and the difference between the actual rates and the rates from the backward extrapolation of Layer 1 was the production from Layer 2. History matching the resulting Layer 2 data yielded the decline parameters(qi2 and N2) for Layer 2. The results were consistent with the geologic description. However, the method is useful only when the well is at the late stage of production, when the production is mainly from one layer. The method cannot be applied early in the well’s life. Cheng et al. (2008) developed a method to improve reserve estimates from the DCA of tight and multilayer gas wells. The b-value is estimated during boundary-dominated flow, and history matching is performed using this b-value in a backward data-fitting fashion to determine the decline parameters (qi and Di) at the end of production. These parameters are then used for production forecasts. However, the method proposed by Cheng et al. (2008) cannot be used to estimate decline parameters for each layer. This method assumes that the performances of a multilayer well can be represented by a single-layer-well model. The production performance for a multilayer gas well is complicated because the production contribution from each layer changes constantly with time (Fetkovich 1980; Cheng et al. 2008). Several previous studies (Arps 1945; Fetkovich et al. 1996; Cheng et al. 2008) showed that DCA for multilayer gas wells could yield an Arps (1945) decline-curve exponent (b) greater than 1.0. Whereas b is constant for a single-layer well, b changes with time for a multilayer well. The nonuniqueness problem from history matching is more prominent in a multilayer system. The available methods could not provide a unique set of decline parameters (qi, Di, and b) for individual layers. The objective of this study is to propose a new method to provide valid production forecasts and reserve estimation for multilayer wells (both oil and gas) on the basis of early-stage production. The method’s innovations are to estimate the instantaneous values of D and b from the production data during boundary-dominated flow and to perform history matching on the profiles of q, D, and b to determine a set of decline parameters (qi, Di, and b) for each layer. Both a synthetic example and field data are used to validate the C 2020 Society of Petroleum Engineers Copyright V This paper (SPE 195085) was accepted for presentation at the SPE Middle East Oil and Gas Show and Conference, Manama, Bahrain, 18–21 March 2019, and revised for publication. Original manuscript received for review 30 May 2019. Revised manuscript received for review 10 August 2019. Paper peer approved 2 September 2019. 2020 SPE Journal 1 proposed method. The nonuniqueness problem from history matching is solved, and the best decline models for individual layers are identified. Consequently, better production forecasts and more accurate reserve estimation can be achieved for both the short- and the long-term. A New Methodology To best illustrate the principles of the method, this study considers a two-layer system. Layer 1 has higher productivity than Layer 2. The main assumptions are constant flowing wellbore pressure, single-phase flow, no aquifer support, no reservoir crossflow, no wellbore crossflow, and no changes in completion and operating conditions. For commingled production at a constant wellbore pressure, production from any layer is not dependent on the production from other layers. Therefore, the well flow rate is the sum of the flow rates from each layer (Russell and Prats 1962). The flow rate of the well for a two-layer system is qðtÞ ¼ q1 ðtÞ þ q2 ðtÞ: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð1Þ The production rate for layer j (Arps 1945) is qij qj ðtÞ ¼ ð1 þ bj Dij tÞ1=bj for j ¼ 1; 2: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð2Þ At t ¼ 0, the initial flow rate of a well is qi ¼ qi1 þ qi2 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð3Þ Following Fetkovich (1980), the decline-curve dimensionless rate of a well is defined as qDd ¼ qðtÞ ¼ qi     qi1 q qDd1 þ i2 qDd2 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð4Þ qi qi The decline-curve dimensionless rate is a weighted average of the decline-curve dimensionless rates of Layers 1 and 2. The weights of Layers 1 and 2 are the initial production contribution of Layers 1 and 2, respectively. Note that the weights are constant. Then the D of the well is     1 dqðtÞ q ðtÞ q ðtÞ ¼ 1 D1 ðtÞ þ 2 D2 ðtÞ; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð5Þ DðtÞ ¼  qðtÞ dt qðtÞ qðtÞ where the D for each layer (Yu 2011; Blasingame and Rushing 2005) is Dj ðtÞ ¼ Dij for j ¼ 1; 2: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð6Þ 1 þ bj Dij t The D of the well is a weighted average of the decline rates of Layers 1 and 2. The weights of Layers 1 and 2 are the instantaneous production contributions of Layers 1 and 2, respectively. These weights are not constant. The b of the well is bðtÞ ¼     dD1 ðtÞ q ðtÞ D21 ðtÞ q2 ðtÞ D22 ðtÞ ¼ 1 ð1 þ b Þ þ ð1 þ b2 Þ  1: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð7Þ 1 dt qðtÞ D2 ðtÞ qðtÞ D2 ðtÞ Whereas Eq. 5 is similar to the result from Spivey et al. (2001), Eq. 7 is different. However, the result from Eq. 7 is perfectly matched with Fig. 2 from Spivey et al. (2001). Fetkovich et al. (1996) showed that the average value of b(t) is not less than the values of b1 and b2. To understand Eqs. 1 through 7, consider a two-layer gas well with the following properties for Layers 1 and 2, respectively: • Initial flow rate (qi) ¼ 5.0 and 2.0 MMscf/D • Di ¼ 0.010 and 0.001 d1 • b ¼ 0.0 and 0.50 During the boundary-dominated flow period, the production profile and the production forecasts based on of the production data of 100 days and 1, 2, and 4 years are shown in Fig. 1. The first forecasting method is from Cheng et al. (2008) (Fig. 1a). The two key components ofs their method are: • The b-value is the averaged value during the production period. • Other decline parameters (qi and Di) are the values at the end of the production period. This set of decline parameters is used for the production forecast. Note that some values of b are greater than 1.0. The forecasts can match the production profile only in the short-term and fail to do so in the longer term. If the value of b is restricted to be less than or equal to 1.0 (b  1), the production forecast is improved significantly. Fetkovich et al. (1990) recommended a b-value of 0.6, which further improves the forecast. The errors in reserve estimations at the end of each production period are also shown in Fig. 1. With complex production, the error percentage does not decrease monotonically with time. The reserve are underestimated during the early stage of production and tend to be overestimated during the late stage of production. Yu (2011) pointed out the slope change (flattening and then increasing again) in the log-log plot. This is caused by a large permeability contrast between layers. With a long production life, the layer with lower permeability dominates the well production in the long-term. It is apparent that we cannot predict the performance of a multilayer system with a significant permeability contrast using a single-layer-well approach. 2 2020 SPE Journal 10 q (MMscf/D) q (MMscf/D) 10 1 Well Forecast at 100 days Forecast at 1 year Forecast at 2 years Forecast at 4 years 0.1 10 Well Forecast at 100 days Forecast at 1 year Forecast at 2 years Forecast at 4 years 1 0.1 100 1,000 0 10,000 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 t (days) t (days) (a) b is calculated from the history matching (b = 0.46 for 100 days, 1.30 for 1 year, 1.61 for 2 years, and 1.10 for 4 years) 10 q (MMscf/D) q (MMscf/D) 10 1 Well Forecast at 100 days Forecast at 1 year Forecast at 2 years Forecast at 4 years 0.1 10 Well Forecast at 100 days Forecast at 1 year Forecast at 2 years Forecast at 4 years 1 0.1 100 1,000 10,000 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 t (days) t (days) (b) b < 1 (b = 0.46 for 100 days, 1.00 for 1 year, 2 years, and 4 years) 10 q (MMscf/D) q (MMscf/D) 10 1 Well Forecast at 100 days Forecast at 1 year Forecast at 2 years Forecast at 4 years 0.1 10 100 10,000 1,000 t (days) Well Forecast at 100 days Forecast at 1 year Forecast at 2 years Forecast at 4 years 1 0.1 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 8,000 t (days) (c) b = 0.6 for all forecasts Production Data EUR Error (%) Historical b b<1 100 days –61 –61 –52 1 year 98 33 –25 2 years 202 70 7 4 years 67 53 7 b = 0.6 (d) Estimated ultimate recovery (EUR) errors Fig. 1—Flow rate (q) and production forecasts for a two-layer gas system. The profiles of q, the decline-curve dimensionless rate (q/qi), D, and b are illustrated in Fig. 2. We can make the following observations: • The production contribution from each layer varies with time. The higher-permeability layer has a relatively high early-time production rate with a steeper decline. The lower-permeability layer declines less rapidly and dominates the later-time production. The production contribution from Layer 2 is increasing with time. This phenomenon was clearly explained by previous studies (Fetkovich et al. 1990; Cheng et al. 2008; Yu 2011). For a system with a large permeability contrast, the early-rate decline is large. Then, the production is dominated by Layer 2, with a much-flatter rate decline (Yu 2011). This characteristic can be identified by the log-log plot of flow rate vs. time. 2020 SPE Journal 3 • As described in Eq. 4, the q/qi of a well is a weighted average of the q/qi of Layers 1 and 2. The weights are constant, not a function of time. As a result, on a plot of q/qi vs. t, the respective distances between the well and each layer are always equal for any point in time, regardless of the instantaneous production contribution from each layer. • D is a weighted average of the decline rate of Layers 1 and 2 (D1 and D2). The weights are not constant. The value of D(t) is always between the values of D1(t) and D2(t). As the contribution of Layer 1 decreases, the value of D(t) approaches the value for D2(t). • Eq. 7 explains the b of a well. The value of b is increasing and then decreasing. The maximum value of b can be greater than 1.0. As the contribution of Layer 1 diminishes, the value of b(t) approaches the value of b2 (Fetkovich et al. 1990). • Fig. 2 clearly shows that the values of D and b are not constant but vary with time. This is the reason the multilayer system performance is much more complicated than that of a single-layer system. 0.012 10 Well Layer 1 Layer 2 0.010 Well Layer 1 Layer 2 D, d 1 q (MMscf/D) 0.008 0.006 1 0.004 0.002 0.1 10 0 100 1,000 10,000 0 200 400 t (days) 600 800 1,000 1,200 1,400 t (days) 3.0 1.00 Well Layer 1 Layer 2 0.75 Well Layer 1 Layer 2 2.5 0.50 α b q /qi 2.0 q1(t ) q(t ) 1.5 1.0 0.25 α 0 0 200 400 0.5 q2(t ) q(t ) 600 800 t (days) 1,000 1,200 1,400 0 0 200 400 600 800 1,000 1,200 1,400 t (days) Fig. 2—Flow rate (q), decline-curve dimensionless rate (q/qi), decline rate (D), and decline-curve exponent (b). The methodology of this study improves production forecasts and reserve estimation from DCA for a multilayer well. The available production data are flow rate [q(t)] and cumulative production [Gp(t)]. The unknown decline parameters for Layers 1 and 2 are • Initial flow rates (qi1, qi2) • Initial decline rates (Di1, Di2) • Decline-curve exponents (b1, b2) To find models for Layers 1 and 2, history matching is performed to match the q(t), D(t), and b(t) of the well. Additional history matching on D(t) and b(t) provides the decline parameters for each layer. Therefore, the new method will solve the nonuniqueness problem. The model will yield good production forecasts and reserve estimation at early stages of production. The proposed method can be used for both oil and gas wells. Applications Example 1—Nonuniqueness Problem from History Matching. The nonuniqueness problem from history matching has been well-documented in the literature. It becomes highly problematic when we try to match performance from a multilayer well with a single-layer system. The objective of this example is to illustrate the nonuniqueness problem for a multilayer system. Let’s consider a two-layer gas well with the properties listed in the “Well” column in Table 1. Layer 1 has a higher initial-production rate, a higher initial-decline rate, and a lower decline-curve exponent than Layer 2. The details of the calculation are in Appendix A. On the basis of these parameters, Layer 1 will be depleted faster than Layer 2. In this case, we assume that we have only 400 days of production data. Four different cases of history matching are considered in this example: • Case I: The initial flow was incorrect, and Layer 2 has the higher initial-production rate. • Case II: The initial decline rate for Layer 1 is overestimated by 12%. • Case III: The initial decline rate for Layer 1 is underestimated by 12%. • Case IV: The decline-curve exponent was incorrect, and Layer 1 has a higher decline-curve exponent. The decline parameters (qi, Di, and b) and the errors in reserve estimation from history matching are reported in Table 1. It is apparent that the value of Di has a significant effect on the estimated reserve. The results of history matching are shown in Figs. 3 through 8. Fig. 3 plots profiles of q for the different cases. It is apparent that these four models all can match the 400 days of production data fairly well. This is the well-known nonuniqueness problem from history matching. However, the four models fail to match well performance 4 2020 SPE Journal in the long-term. Plots of q vs. Gp for the different cases are illustrated in Fig. 4. These four models predict significantly different ultimate recovery from the actual performance regardless of apparently good history matching. History matching on q and Gp leads to the nonuniqueness problem. Case I: qi1 < qi2 Well Case II: Di1 is 12% Overestimated Case IV: b1< b2 Case III: Di1 is 12% Underestimated Parameters Layer 1 Layer 2 Layer 1 Layer 2 Layer 1 Layer 2 Layer 1 Layer 2 Layer 1 Layer 2 qi 5.0 2.0 2.0 5.0 5.0 2.0 5.0 2.0 5.0 2.0 Di 0.0045 0.0010 0.0025 0.0040 0.0050 0.0005 0.0040 0.0018 0.0029 0.0052 b 0.0 0.4 0.0 0.4 0.0 0.4 0.0 0.4 0.4 0.0 Reserve (MMscf) 2,883 Reserve error (%) — — 1,305 6,096 1,557 1,728 –55% 111% –46% –40% Table 1—Decline parameters (qi, Di, and b) from history matching for Example 1. 10 10 Case III q (MMscf/D) q (MMscf/D) Case I 1 1 Well Well History matching History matching Forecast Forecast 0.1 10 100 t (days) 0.1 10 10,000 1,000 100 t (days) 1,000 10 10 Case IV q (MMscf/D) q (MMscf/D) Case II 1 1 Well Well History matching History matching Forecast 0.1 10 10,000 Forecast 100 t (days) 1,000 10,000 0.1 10 100 t (days) 1,000 10,000 Fig. 3—Flow rate (q) for Example 1. The profiles of q/qi for different cases are illustrated in Fig. 5. The models for all cases match the q/qi of a well but fail to match the q/qi of individual layers. Note that in the actual DCA, we cannot see the actual profile of the q/qi of individual layers. The profiles of D, b, and (q2/q) are illustrated in Figs. 6 through 8, respectively. Note that the profile of q2/q can be obtained when production-logging data are available. Different models yield different profiles of these parameters. Therefore, we could use this set of parameters to solve the nonuniqueness problem and to identify the best model for production forecasting and for reserve estimation. The results of this example illustrate that history matching on only q and Gp might yield misleading results. Several different models can have apparently good history matching on the short-term production data but yield different production forecasts and reserve estimations. The profiles of D and b have a greater ability to identify the valid model than the profile of q. Example 2—An Oil Well from Coalinga Oil Field. The data in this example is from the East Side Coalinga Field, California, USA. There are isolated (no crossflow) upper and lower oil sands. The departure-curve method (differencing) is applied to analyze the production data (Fetkovich 1980). Note that the departure-curve method requires the last portion of the production data to differentiate the performances of upper and lower layers. This is the period when only one layer significantly contributes to the production. Therefore, (departure-curve method) cannot be applied to analyze production data during the early stages of well life when both layers significantly contribute to the production. The result is the following oil rate: qðtÞ ¼ 58;824 50;000 þ 0:535t ; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð8Þ e0:200t e where q(t) is the total production rate from both layers in bbl/yr and t is time in years. 2020 SPE Journal 5 10 10 Case III q (MMscf/D) q (MMscf/D) Case I 1 1 Well History matching Forecast Well History matching Forecast 0.1 0.1 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 Gp (MMscf) Gp (MMscf) 10 10 Case IV q (MMscf/D) q (MMscf/D) Case II 1 1 Well Well History matching Forecast History matching Forecast 0.1 0.1 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 Gp (MMscf) Gp (MMscf) Fig. 4—Flow rate (q) vs. cumulative gas production (Gp) for Example 1. 1.0 1.0 Case I 0.8 0.7 0.7 0.6 0.6 0.5 Layer 1 Layer 2 Well Layer 1, matching Layer 2, matching Well, matching 0.4 0.3 0.2 0.1 0 50 100 150 Case III 0.9 0.8 q /q i q /q i 0.9 0.5 Layer 1 Layer 2 Well Layer 1, matching Layer 2, matching Well, matching 0.4 0.3 0.2 0.1 200 250 300 350 0 400 50 100 150 t (days) 250 300 350 400 1.0 1.0 Case II 0.9 0.8 0.7 0.7 0.6 0.6 0.5 Layer 1 Layer 2 Well Layer 1, matching Layer 2, matching Well, matching 0.4 0.3 0.2 0.1 0 50 100 150 Case IV 0.9 0.8 q /qi q /qi 200 t (days) 0.5 Layer 1 Layer 2 Well Layer 1, matching Layer 2, matching Well, matching 0.4 0.3 0.2 0.1 200 t (days) 250 300 350 400 0 50 100 150 200 250 300 350 400 t (days) Fig. 5—Decline-curve dimensionless rate (q/qi) for Example 1. 6 2020 SPE Journal 0.004 0.004 Well History matching Forecast 0.003 D (days–1) 0.003 D (days–1) Case I Well History matching Forecast 0.002 0.001 0.002 0.001 0 0 0 400 800 1,200 1,600 2,000 2,400 2,800 3,200 3,600 4,000 0 400 800 1,200 1,600 2,000 2,400 2,800 3,200 3,600 4,000 t (days) t (days) 0.004 0.004 Well History matching Forecast Case II Well History matching Forecast 0.003 D (days–1) 0.003 D (days–1) Case III 0.002 Case IV 0.002 0.001 0.001 0 0 0 400 800 1,200 1,600 2,000 2,400 2,800 3,200 3,600 4,000 0 400 800 1,200 1,600 2,000 2,400 2,800 3,200 3,600 4,000 t (days) t (days) Fig. 6—Decline rate (D) for Example 1. 3.0 3.0 Case III 2.5 2.5 2.0 2.0 1.5 b b Case I Well History matching Forecast 1.0 1.5 Well History matching Forecast 1.0 0.5 0.5 0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 t (days) t (days) 3.0 3.0 Case II 2.5 2.0 2.0 1.5 b b Case IV 2.5 Well History matching Forecast 1.0 1.5 Well History matching Forecast 1.0 0.5 0.5 0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 t (days) t (days) Fig. 7—Decline-curve exponent (b) for Example 1. Fetkovich (1980) mentioned that a single-layer model with a b-value of 0.2 can match nearly all of the data points but cannot be explained by any of the drive mechanisms in this field. We tried history matching with a single layer with the following result: qðtÞ ¼ 108;824 ½1 þ 0:2  0:353 t 1=0:2 2020 SPE Journal : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð9Þ 7 1.0 1.0 Case III 0.8 0.8 0.6 0.6 q2 /q q2 /q Case I 0.4 Well History matching Forecast 0.2 500 Well History matching Forecast 0.2 0 0 0.4 1,000 1,500 2,000 2,500 3,000 3,500 4,000 0 0 500 t (days) t (days) 1.0 1.0 Case II 0.8 0.6 q2 /q q2 /q Case IV 0.8 0.6 0.4 Well History matching Forecast 0.2 0 1,000 1,500 2,000 2,500 3,000 3,500 4,000 0.4 Well History matching Forecast 0.2 0 0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 t (days) 0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 t (days) Fig. 8—Production contribution from Layer 2 (q2/q) for Example 1. The proposed method is also applied to analyze the production from this field. The details of the calculation are in Appendix B. It is apparent that the production data are very noisy. Therefore, we perform regression on q and t, yielding a profile for the oil rate. This profile is then used to establish the profiles for D and b. History matching on q, D, and b simultaneously yields the following model: qðtÞ ¼ 61;732 42;083 þ 0:181t : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð10Þ e0:462t e The profiles of q, D, and b and the production contribution from Layer 2 (q2/q) are illustrated in Fig. 9. Note that the profile of b is slightly distorted by the polynomial function used in the regression. It is clearly seen that the model in this study can better fit the profiles of q, D, and b than the ones from the departure-curve method and from a single layer. In addition, for the model from the departure-curve method, Layer 2 has a lower initial contribution (qi2 < qi1) but was depleted faster (Di2 > Di1). Layer 1 with higher productivity dominated the production during the late-production period. However, the model in this study indicates that Layer 2 with lower productivity should be depleted more slowly (Di2 < Di1) and should dominate the production during the late period. Discussion For a Multilayer Well. Fetkovich (1980) recommended that two or more layers can be perfectly represented by a single layer only when these layers share the same value of the decline-curve dimensionless rate (qDd). 8 > ð1 þ bDi tÞ1=b for 0 < b  1 qðtÞ < qDd ¼ ¼ > qi : Di t e for b ¼ 0:                                                 ð11Þ In other words, these layers will be depleted at the same rate when they have the same values of Di and b. The q/qi is more sensitive to the value of Di than to the value of b. For practical purposes, we might concentrate on Di and ignore the effect of b. Fetkovich et al. (1996) recommended reducing a multilayer system to a two-layer system by combining layers with similar values of qi/G. Many studies recommend no more than two layers for history matching because that would complicate the matching process and magnify the nonuniqueness problem. Noise on Production Data. The proposed method uses the parameters D and b, which involve the derivative of production data. The effect of inherent noise in the production data is magnified by taking the derivative (Ilk et al. 2011; Varma et al. 2018). Therefore, data smoothing and/or modern differentiation algorithms are required to improve the derivative calculations (Kupchenko et al. 2008; Varma et al. 2018). The rate-cumulative profile is always smoother than the rate-time profile. Ilk et al. (2008) recommended the following alternative to compute the values of D and b for noisy production data: D¼ b¼q 8 dq ; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð12Þ dG dD1 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð13Þ dG 2020 SPE Journal 0.40 100,000 Production Departure-curve method Single-layer model Model of this study D (years–1) q (BOPY) 0.35 10,000 Production Departure-curve method Single-layer model Model of this study 0.30 0.25 0.20 0.15 1,000 0 4 2 6 8 10 12 14 16 0 2 4 6 t (years) 0.30 10 12 14 16 1.00 Production Departure-curve method Single-layer model Model of this study 0.25 8 t (years) 0.75 q2 /q b 0.20 0.15 Model Departure-curve method 0.50 0.10 0.25 0.05 0 0 0 2 4 6 8 10 12 14 16 t (years) 0 2 4 6 8 10 12 14 16 t (years) Fig. 9—Flow rate (q), decline rate (D), decline-curve exponent (b), and production contribution from Layer 2 (q2/q) for Example 2. Only Data in the Hyperbolic Window Can Be Used for History Matching. The decline-curve exponent (b) in the Arps (1945) hyperbolic decline is assumed to be constant throughout a well’s life. Stumpf and Ayala (2016) showed that b is approximately constant only in the hyperbolic window. This window occurs during the early boundary-dominated flow period. Furthermore, the value of b approaches zero during the late boundary-dominated flow period. Therefore, only the data in the hyperbolic window can be used for history matching. Prediction outside this window yields invalid matching results because the requirement that b is approximately constant is not honored. Production Logging (PL). PL is a casedhole logging technique that can be used to allocate production on a layer-by-layer basis. In addition to defining zonal flow contributions at a given point in time, PL can define the inflow performance relationship and average drainage area pressure separately for each zone on the basis of the zones’ individual flow contributions at different total-well flow rates. PL surveys can be repeated multiple times and used to define the variation of layer-flow contributions with time. If the magnitude of inherent noise in the production data is too large and the profiles of D and b cannot be established, PL can be applied to solve the nonuniqueness problem. Conclusions 1. Both D and b of a multilayer system are not constant but vary with time. The production profile of this system cannot be modeled by a one-layer system (with a constant value of b), especially when there is a significant contrast in permeability. Consequently, the conventional DCA method yields unreliable production forecasts and reserve estimation. 2. The nonuniqueness problem of history matching for a multilayer system is illustrated. Many different one-layer models can match production from a multilayer system, especially when a well is at an early stage of production. The model with the correct decline parameters (qi, Di, and b) in each layer could not be identified from the methods available in the literature. 3. This study proposes a new method to improve production forecasts and reserve estimation and to solve the nonuniqueness problem for a well in a multilayer system. With this new methodology, history matching is performed not only on q but also on D and b. This method yields the set of decline parameters (qi, Di, and b) of each individual layer at early stages of production. 4. The proposed method is successfully validated using several synthetic examples and an actual field data set. 5. Improvements in production forecasting and reserve estimation, especially during early stages of production, can be achieved by applying the proposed method. Nomenclature b ¼ decline-curve exponent, dimensionless bE ¼ equivalent or averaged decline-curve exponent, dimensionless cg ¼ isothermal gas compressibility, psi1 D ¼ decline rate, day–1 G ¼ ultimate gas recovery, MMscf Gp ¼ cumulative gas production, MMscf 2020 SPE Journal 9 n N p q qDd t UR z lg ¼ ¼ ¼ ¼ ¼ ¼ ¼ ¼ ¼ backpressure curve exponent, dimensionless ultimate oil recovery, STB pressure, psia flow rate, MMscf/D or bbl/yr decline-curve dimensionless rate, dimensionless time, days ultimate recoverable oil or gas, STB or MMscf gas compressibility factor, dimensionless gas viscosity, cp Subscripts 1 ¼ Layer 1 2 ¼ Layer 2 i ¼ initial condition wf ¼ well flowing References Arps, J. J. 1945. Analysis of Decline Curves. Trans AIME 160 (1): 228–247. SPE-945228-G. https://doi.org/10.2118/945228-G. Ayala, L. F. and Ye, P. 2013. Density-Based Decline Performance Analysis of Natural Gas Reservoirs Using a Universal Type-Curve. J Energy Resour Technol 135 (4): 042701. JERT-12-1268. https://doi.org/10.1115/1.4023867. Blasingame, T. A. and Rushing, J. A. 2005. A Production-Based Method for Direct Estimation of Gas-in-Place and Reserve. Paper presented at the SPE Eastern Regional Meeting, Morgantown, West Virginia, USA, 14–16 September. SPE-98042-MS. https://doi.org/10.2118/98042-MS. Chen, H.-Y. and Teufel, L. W. 2002. Estimating Gas Decline-Exponent Before Decline-Curve Analysis. Paper presented at the SPE Gas Technology Symposium, Calgary, Alberta, Canada, 30 April–2 May. SPE-75693-MS. https://doi.org/10.2118/75693-MS. Cheng, Y., Lee, W. J., and McVay, D. A. 2008. Improving Reserve Estimates from Decline-Curve Analysis of Tight and Multilayer Gas Wells. SPE Res Eval & Eng 11 (5): 912–920. SPE-108176-PA. https://doi.org/10.2118/108176-PA. Fetkovich, M. J. 1980. Decline Curve Analysis Using Type Curves. J Pet Tech 32 (6): 1065–1077. SPE-4629-PA. https://doi.org/10.2118/4629-PA. Fetkovich, M. J., Bradley, M. D., Works, A. M. et al. 1990. Depletion Performance of Layered Reservoirs Without Crossflow. SPE Form Eval 5 (3): 310–318. SPE-18266-PA. https://doi.org/10.2118/18266-PA. Fetkovich, M. J., Fetkovich, E. J., and Fetkovich, M. D. 1996. Useful Concepts for Decline-Curve Forecasting, Reserve Estimation, and Analysis. SPE Res Eng 11 (1): 13–22. SPE-28628-PA. https://doi.org/10.2118/28628-PA. Ilk, D., Jenkins, C. D., and Blasingame, T. A. 2011. Production Analysis in Unconventional Reservoirs—Diagnostics, Challenges, and Methodologies. Paper presented at the North American Unconventional Gas Conference and Exhibition, The Woodlands, Texas, USA, 14–16 June. SPE-144376-MS. https://doi.org/10.2118/144376-MS. Ilk, D., Rushing, J. A., Perego, A. D. et al. 2008. Exponential vs. Hyperbolic Decline in Tight Gas Sands: Understanding the Origin and Implications for Reserve Estimates Using Arps’ Decline Curve. Paper presented at the SPE Annual Technical Conference and Exhibition, Denver, Colorado, USA, 21–24 September. SPE-116731-MS. https://doi.org/10.2118/116731-MS. Johnson, R. H. and Bollens, A. L. 1927. The Loss Ratio Method of Extrapolating Oil Well Decline Curves. Trans AIME 77 (1): 771–778. SPE-927771-G. https://doi.org/10.2118/927771-G. Kupchenko, C. L., Gault, B. W., and Mattar, L. 2008. Tight Gas Production Performance Using Decline Curves. Paper presented at the CIPC/SPE Gas Technology Symposium Joint Conference, Calgary, Alberta, Canada, 16–19 June. SPE-114991-MS. https://doi.org/10.2118/114991-MS. Lee, W. J., Rollins, J. B., and Spivey, J. P. 2003. Pressure Transient Testing, Vol. 9. Richardson, Texas, USA: Textbook Series, Society of Petroleum Engineers. Okuszko, K. E., Gault, B. W., and Mattar, L. 2007. Production Decline Performance of CBM Wells. Paper presented at the Canadian International Petroleum Conference, Calgary, Alberta, Canada, 12–14 June. PETSOC-2007-078. https://doi.org/10.2118/2007-078. Russell, D. G. and Prats, M. 1962. Performance of Layered Reservoirs with Crossflow—Single-Compressible Fluid Case. SPE J. 2 (1): 53–67. SPE-99-PA. https://doi.org/10.2118/99-PA. Spivey, J. P., Frantz, J. H., Williamson, J. R. et al. 2001. Applications of the Transient Hyperbolic Exponent. Paper presented at the SPE Rocky Mountain Petroleum Technology Conference, Keystone, Colorado, USA, 21–23 May. SPE-71038-MS. https://doi.org/10.2118/71038-MS. Stumpf, T. N. and Ayala, L. F. 2016. Rigorous and Explicit Determination of Reserve and Hyperbolic Exponents in Gas-Well Decline Analysis. SPE J. 21 (5): 1843–1857. SPE-180909-PA. https://doi.org/10.2118/180909-PA. Varma, S., Tabatabaie, S. H., Ewert, J. et al. 2018. Variation of Hyperbolic-b-Parameter for Unconventional Reservoirs, and 3-Segment Hyperbolic Decline. Paper presented at the SPE/AAPG/SEG Unconventional Resources Technology Conference, Houston, Texas, USA, 23–25 July. URTEC2892966-MS. https://doi.org/10.15530/URTEC-2018-2892966. Ye, P. and Ayala, L. F. 2013. Straight-Line Analysis of Flow Rate vs. Cumulative-Production Data for the Explicit Determination of Gas Reserve. J Can Pet Technol 52 (4): 296–305. SPE-165583-PA. https://doi.org/10.2118/165583-PA. Yu, S. 2011. An Improved Methodology to Obtain the Arps’ Decline Curve Exponent (b) for Tight/Stacked Gas Reservoirs. Paper presented at the North American Unconventional Gas Conference and Exhibition, The Woodlands, Texas, USA, 14–16 June. SPE-143907-MS. https://doi.org/10.2118/143907-MS. Appendix A—Details of Calculation for Example 1 The production data in Example 1 are a synthetic data set that is based on the following model: q1 ðtÞ ¼ q2 ðtÞ ¼ 5 expð0:0045tÞ 2 ð1 þ 0:4  0:001tÞ1=0:4 5 2 þ qðtÞ ¼ : expð0:0045tÞ ð1 þ 0:4  0:001tÞ1=0:4 10                                               ðA-1Þ 2020 SPE Journal The D for Layers 1 and 2 and wells, respectively, is D1 ðtÞ ¼ 0:0045; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðA-2Þ 0:001 ; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðA-3Þ 1 þ 0:0004 t     q ðtÞ 0:001 q2 ðtÞ DðtÞ ¼ 0:0045 1 þ : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðA-4Þ qðtÞ 1 þ 0:0004 t qðtÞ D2 ðtÞ ¼ The instantaneous decline-curve exponent is bðtÞ ¼     q1 ðtÞ D21 ðtÞ q2 ðtÞ D22 ðtÞ þ 1:4  1: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðA-5Þ qðtÞ D2 ðtÞ qðtÞ D2 ðtÞ Because there is no noise in this production data, there is no need to use the polynomial function to represent the data. We assume that there are only 400 days of production data. The objective of this example is to illustrate the nonuniqueness problem. We can find at least four different models (Cases I, II, III, and IV), that can apparently match the production data, q(t). These models are qI ðtÞ ¼ 2 5 þ expð0:0025 tÞ ð1 þ 0:4  0:004 tÞ 1=0:4 qII ðtÞ ¼ 5 2 þ expð0:005 tÞ ð1 þ 0:4  0:0005 tÞ 1=0:4 qIII ðtÞ ¼ 5 2 þ expð0:004 tÞ ð1 þ 0:4  0:0018 tÞ 1=0:4 qIV ðtÞ ¼ 5 ð1 þ 0:4  0:0029 tÞ 1=0:4 þ 2 ðA-6Þ expð0:0052 tÞ                                               For each case, D can be calculated as DðtÞ ¼ Di : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðA-7Þ 1 þ bDi t The b for each case is defined as bðtÞ ¼ dD1 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðA-8Þ dt Appendix B—Details of Calculation for Example 2 The production data in Table B-1 is from the East Side Coalinga Field, California, USA. Years bbl/yr 0.5 90,000 1.5 64,000 2.5 48,000 3.5 36,000 4.5 27,500 5.5 21,250 6.5 16,250 7.5 13,000 8.5 10,500 9.5 8,500 10.5 6,900 11.5 5,600 12.5 4,550 13.5 3,800 14.5 3,200 15.5 2,750 Table B-1—Production data from the East Side Coalinga Field. 2020 SPE Journal 11 Fetkovich (1980) applied the departure-curve method (differencing) to analyze the production data and yielded the following model: qðtÞ ¼ 58;824 50;000 þ 0:535t : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðB-1Þ e0:200t e We tried history matching with a single layer [b ¼ 0.2, recommended by Fetkovich (1980)]. The result was qðtÞ ¼ 108;824 ð1 þ 0:2  0:353 tÞ1=0:2 : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðB-2Þ Because of the noisy production, this study recommends using the polynomial function to represent the actual production data. The regression yields the following result (Fig. B-1): lnðqÞ ¼ 0:000000878t6  0:0000389t5 þ 0:000660t4  0:00555t3 þ 0:0307t2  0:375t þ 11:6 : . . . . . . . . . . . . . . . . . . . ðB-3Þ 12.0 Production Polynomial (production) 11.5 11.0 10.5 ln(q ) 10.0 9.5 9.0 8.5 8.0 y = 8.78×10–7x 6 – 3.89×10–5x 5 + 6.60×10–4x 4 – 5.55×10–3x 3 + 3.07×10–2x 2 – 3.75×10–1x + 1.16×101 R 2 = 1.00×100 7.5 7.0 0 2 4 6 8 10 12 14 16 t (years) Fig. B-1—Profile of ln(q). With this smooth production profile, the profile of D(t) can be easily calculated as dln½qðtÞ dt ¼ 6  0:000000878t5 þ 5  0:0000389t4  4  0:000660t3 þ 3  0:00555t2  2  0:0307t þ 0:375:         ðB-4Þ DðtÞ ¼  Then [1/D(t)] is calculated and plotted in Fig. B-2. 6.000 5.000 1/D 4.000 3.000 y = 0.000050576x 4 – 0.002554848x 3 + 0.033381547x 2 + 0.044476090x + 3.102211202 R 2 = 0.998423945 2.000 1.000 0 2 4 6 8 t 10 12 14 16 18 Fig. B-2—Profile of (1/D). 12 2020 SPE Journal The result of the regression on (1/D) is 1 ¼ 0:0000506t4  0:00255t3 þ 0:0334t2 þ 0:0445t þ 3:10: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ðB-5Þ D The profile can be calculated as dð1=DÞ dt ¼ 4  0:0000506t3  3  0:00255t2 þ 2  0:0334t þ 0:0445: bðtÞ ¼                                   ðB-6Þ Kittiphong Jongkittinarukorn is a lecturer in the Department of Mining and Petroleum Engineering, Chulalongkorn University, Thailand. His research interests are decline curve analysis, material balance, and commingled production. Jongkittinarukorn holds a PhD degree in petroleum engineering from the University of Oklahoma and a DBA in finance from Chulalongkorn University, Thailand. He is a member of SPE. Nick Last is a consulting engineer specializing in the analysis of dynamic well and reservoir data, especially well tests and production logs. He began his career with Schlumberger in the North Sea and Southeast Asia, and has subsequently delivered consulting projects for numerous major operators and government bodies worldwide. Last is a founder of Well Test Knowledge International (WTKI) and has taught numerous classes under the WTKI banner. He holds a BSc degree (with honors) in physics from Imperial College, London. Freddy Humberto Escobar is a professor at Universidad Surcolombiana. His field of research is well-test analysis. Escobar holds MSc and PhD degrees, both in petroleum engineering from the University Surcolombiana. He joined SPE in 1987 and served as a member of the SPE Formation Evaluation Award Committee. Escobar is the faculty sponsor of the SPE student chapter at Universidad Surcolombiana. Kreangkrai Maneeintr is a lecturer in the Department of Mining and Petroleum Engineering, Chulalongkorn University, Thailand. His areas of interest are energy and environment, especially carbon capture, storage, and utilization; coal waste management; and heavy-oil production. Maneeintr holds a PhD degree from the University of Regina, Saskatchewan, Canada. 2020 SPE Journal 13