Shale technology: Bayesian variable pressure decline-curve analysis for shale gas wells

Decline-curve analysis is the most popular technique in the oil and gas industry to forecast the production and estimate reserves from wells. Popular decline-curve models include the Arps hyperbolic model, as well as the more recent proliferation of rate-time models developed specifically for unconventional wells, such as the power-law exponential, stretched exponential and the logistic growth models. This latter category also includes physics-based models. 

All these models have one thing in common: the assumption of a constant BHP. While this assumption is often satisfactory for many conventional oil and gas wells, it is not suitable for many unconventional reservoirs. The application of DCA without properly correcting for variable BHP conditions might lead to several inaccuracies, including—but not limited to—grossly erroneous EURs and incorrect flow regime characterization. 

There are two main ways to account for variable BHP conditions: rate-pressure deconvolution and time superposition. The first approach transforms a variable-pressure rate into a constant-pressure rate response. Once we have the constant-pressure response, we fit the rate-time model to it and then apply convolution to include back pressure variations to history-match and forecast production. Kuchuk (2005) presented an inverse scheme that uses an explicit positivity constraint transforming the linear inversion into non-linear problem. Ilk (2007) introduced a rate-pressure deconvolution approach, based on the cumulative production. They propose an inversion scheme, using b-spline basis functions with regularization to obtain the constant unit-pressure rate solution.  

Deconvolution approaches have an important limitation. Errors in the data (rate, initial reservoir pressure, and BHP) lead to a convolution matrix that is usually rank-deficient and ill-conditioned, making the constant-pressure rate response highly oscillatory and unstable. The second approach uses a decline-curve model, along with the time superposition principle. Ilk and Blasingame (2013) and later Collins (2014) propose applying time superposition, along with decline-curve models, to history-match and forecast production of unconventional wells subject to variable BHP. The authors note the limitations of the method by stating that the superposition principle is only applicable under consistent rate and pressure data. Errors in rate and pressure data lead to poor production history matches and forecasts. 

Bayesian inference uses probability to model uncertainty and variation. The basic ideas behind Bayesian inference are the following: a) it treats the parameters of a given model as random variables and estimates their probability distributions; and b) it incorporates our beliefs and prior information, such as expert knowledge or analogs. 

Gong (2014) was the first to apply Bayesian inference to unconventional wells of the Eagle Ford and Barnett shales, using Arps hyperbolic model. Subsequently, Gonzalez (2012, 2013) showed the application of Bayesian inference to decline analysis with several empirical models: modified Arps power-law exponential, stretched exponential, Duong and logistic growth. Fulford (2016) used the transient hyperbolic relation, a semi-empirical model, to generate distributions of production forecasts of wells in the Elm Coulee field of the Bakken formation.  

Holanda (2018) applies the Bayesian inference, using the analytic solution of the single-phase slightly compressible pressure diffusivity equation to quantify the uncertainty of future production forecasts of Barnett shale wells. The mentioned Bayesian DCA studies all used the Metropolis-Hastings (M-H) MCMC algorithm to sample the posterior distribution of the decline-curve models’ parameters. However, the M-H MCMC technique converges very slowly to stationary posterior distributions when the model’s parameters are correlated. 

GOAL OF WORK 

The objectives of this work are: a) to extend the development and application of the Ruiz Maraggi (2023) variable pressure DCA method for tight-oil reservoirs to shale gas wells in a fast, simple and robust way; and b) to quantify the uncertainty in the model’s parameters and future production predictions. The technique is divided into two stages: I) variable pressure DCA; and II) probabilistic variable pressure DCA. Stage I sequentially estimates: 1) the decline-curve model parameters, 2) the BHP; and 3) the initial reservoir pressure to accurately history-match the production of shale gas wells subject to variable BHP. Stage II computes the posterior distributions of the model’s parameters and its future production prediction. This work uses an adaptive M-H MCMC algorithm that provides fast convergence to stationary posterior distributions. 

This article is organized as follows: First, we detail the steps of the proposed approach. Second, we validate the technique with a synthetic case that has errors in the initial reservoir pressure and BHP history. Finally, this work compares the results of the application of the probabilistic variable pressure drop DCA and traditional DCA for shale gas wells. We illustrate these results, using a physics-based real gas numerical decline-curve model. 

METHODS

Figure 1 is a schematic illustrating the main idea of the variable pressure DCA method, and it depicts two pressure domains: constant and variable pressure. Figure 1a depicts the constant-pressure domain, where the decline-curve model’s rate history is customarily displayed in terms of the unit-pressure rate, q_up(t;x), which is a function of time t and parameters x.

Figure 1b shows the variable-pressure domain, where the actual rate history is displayed. The goal of this step is to map a decline-curve model from a constant pressure drop domain to variable pressure drop domain while yielding a good match of the actual rate history. 

Fig. 1. Schematic illustrating the main concept of the variable pressure DCA method. The goal is to map a decline-curve model from (a) a constant (pseudo) pressure drop domain to (b) variable (pseudo) pressure conditions to history-match and forecast the production of shale gas wells.

Figure 2 illustrates the problems that arise when applying a decline-curve model from a constant-pressure domain to a variable-pressure domain for a real shale gas well example. The model’s history-match (dashed blue lines) to production data (red dots) is extremely poor because of errors in the data. Incorrect estimates of either the initial reservoir pressure P_i and the BHP produce errors in the rate response. These errors propagate forward in time, leading to large oscillations in the variable-pressure rate response of the decline-curve model, as shown by the dashed blue curve in Fig. 2.  

Fig. 2. Gas rate production history of a shale gas well (red dots) and production history-match, using a decline-curve model (dashed blue lines). The model’s history-match (dashed blue lines) to production data (red dots) is poor because of errors in the data. Incorrect estimates in the initial reservoir pressure Pi and the BHP produce errors in the rate response. These errors propagate forward in time leading to large oscillations in the variable-pressure rate response of the decline-curve model.

To avoid these problems, we propose an alternative approach, which was introduced for tight-oil wells in our companion paper, Ruiz Maraggi (2023). This technique provides a fast and simple solution to incorporate variable BHP into any decline-curve model while addressing the following problems (data issues): 

1. Incorrect estimate of initial reservoir pressure P_i
2. Incorrect values of bottomhole flowing pressures P_wf
3. Rate and bottomhole flowing pressure inconsistencies (e.g. flowing pressure increases and gas rate increases). 

There are four main goals of the present method: 

1. History-match the gas production, using a decline-curve model while accounting for variable pressure drop effects. 
2. Estimate the BHP history, consistent with the decline-curve model history-match. 
3. Estimate the initial reservoir pressure, consistent with the decline-curve model history-match. 
4. Quantify the uncertainty in the model’s estimates: estimate distributions of parameters and EUR. 

WORKFLOW

We divide the workflow into two stages: I) variable pressure DCA; and II) probabilistic variable pressure DCA. Ruiz Maraggi (2023) introduces the variable pressure DCA method for the case of a slightly compressible fluid. This study adapts and extends this method for the case of a compressible gas. 

Stage I: Variable pressure DCA.This stage consists of three main steps: 1) production history-match; 2) (pseudo) pressure drop correction; and 3) initial reservoir (pseudo) pressure correction. 

Step 1: Initialization. The first step selects a decline-curve model q_up for the analysis of the gas production. The decline-curve model should be defined in terms of rate per unit pressure drop (dimensions of rate divided by pressure). In addition, the method requires the preliminary estimates (l=0) of the initial reservoir pressure P⁽⁰⁾_i and the bottomhole flowing pressure P⁽⁰⁾_wf. The algorithm requires a reasonable estimate for the BHP history. The bottomhole pressures are usually estimated from wellhead pressures, using flow correlations. 

Step 2: Convert pressures to pseudo-pressures. We must convert the initial and bottomhole flowing pressures into pseudo-pressures. Equation 1 gives the relationship between pressures and pseudo-pressures. The following definition of pseudo-pressure is convenient, since it has units of pressure. 

This study uses the Peng-Robinson 78 EOS using Van der Waals mixing rule to compute the compressibility factor of the natural gas mixture (Sandler 2017), and the Lee-Gonzalez correlation to compute the viscosity of the mixture. After we calculate the pseudo-pressures, we compute the pseudo-pressure drop for each time using Equation 2: 

Step 3: G_p history matching. The third step involves history-matching the gas production, using the time superposition equation with the cumulative production of the decline-curve model G_pup. The cumulative production monotonically increases with time. In addition, it is smoother than the rate history and provides a better history-match while dampening the effects of potential errors in the pressure drop estimates.

Equation 3 assumes that the system behaves linearly with time. The variations of the gas compressibility and viscosity with pressure might introduce non-linearities with time. Agarwal (1979) introduces a pseudo-time transformation to account for this nonlinearity. However, the definition of a pseudo-time requires prior knowledge of the average reservoir pressure at each production time, complicated by the fact that the average reservoir pressure changes with the flow regime(s) present in the reservoir. For these reasons, the application of pseudo-time transformation is not straightforward. 

Equation 4 is the least-squares regression of the cumulative production. 

The goal of this step is to estimate the decline-curve model parameters x^(l+1) and to use them in next step to correct/estimate the pseudo-pressure drop. 

Step 4: Pseudo-pressure drop correction. This step applies rate-(pseudo) pressure deconvolution (Kuchuk et al. 2010) to estimate the pseudo-pressure drop ΔP_wf=ψi−ψ_wf using the decline-curve model q_up and its parameters x^(l+1) from step 3. Equation 5 is an integrated form of the convolution integral. 

Considering a constant pseudo-pressure drop for each time interval, Eq. 5 is discretized as follows. 

where G_p and Δψ_wf are the cumulative production and the pseudo-pressure drop vectors, respectively. The matrix B^(l+1)  contains the time integrals of the decline-curve model q_up. 

The goal of this step is to estimate the pseudo-pressure drop Δψ_wf^(l+1) . This estimate is then used in step 5 to correct/estimate the initial reservoir pseudo-pressure. 

Step 5: Initial reservoir pseudo-pressure correction. This step estimates the initial reservoir pseudo-pressure, using the following least-squares optimization. 

Once Equation 9 is solved, we return to step 3 with the estimated initial reservoir pseudo-pressure, and we repeat steps 3 to 5 until: 

The final outputs of the algorithm are: (a) the model’s parameters x^(l+1), (b) the initial reservoir pseudo-pressure ψ_i^(l+1), and (c) the pseudo-pressure drop Δψ_wf^(l+1). 

Step 6: Compute rate history. This step history-matches the gas rate production, using Equation 1 with the estimated model’s parameters, the initial reservoir pseudo-pressure, and the pseudo-pressure drop. We convert the pseudo-pressure drops into bottomhole flowing pseudo-pressures as follows:

ψ_wf^(l+1)=−Δψ_wf^(l+1)+ψ_i^(l+1)

Stage II: Probabilistic variable pressure DCA. The goal of stage two is to estimate probability distributions of the model’s parameters to generate probabilistic production history-matches and forecasts of shale gas wells. Bayesian inference computes the probability distributions of the model’s parameters given the data p(x|q) (Gelman et al. 2014). 

Prior distribution. π(x) is the distribution of the model’s parameters x, before observing the data. This term incorporates our knowledge or belief into the computation of the posterior probability. 

Likelihood function. The likelihood function f(q|x) is the probability of generating the data q if the parameters of the model were equal to x. This work assumes that the observations of the data q_k follow a Gaussian distribution that is centered around the model value q_k with constant variance σ² (homoscedasticity). 

This work assumes that the sample standard deviation of the residuals ε∗ between the actual production and the calculated production has a normal distribution N(0,σ). 

where σ is the sample standard deviation of the residuals of the least-squares fit (Eq. 14) and ε_k∗ is the residual between the oil rate and the model evaluated at the values of the proposed parameters at the production time k:ε_k∗ =q_k−qˆ_k(x∗). 

where d is the number of parameters of the decline-curve model. The summations in Eq. 13 and Eq. 14 extend over all the production data points k. The superscript I denotes the parameter’s values determined by the optimization routine of stage I. 

Evidence term. The denominator in Eq. 11 is the marginal distribution of the data. In our case, f(q). This term ensures closure; namely, the summation of all probabilities adds to one. 

Posterior distribution. Bayesian inference obtains the posterior distribution p(x|q) of the model’s parameters x after observing the data q. It is not possible to determine analytically the evidence term in Equation 12. However, since the evidence term (denominator of Eq. 11) is a constant, the posterior, likelihood, and prior can be written as a proportionality: 

MCMC methods work by exploiting the product of the prior and likelihood to estimate the posterior distribution of the model’s parameters. They build a sequence of random states x₁→x₂→…→x_s that uses a conditional probability function [Equation] to determine the next state x_s conditional to the previous state x_s−1. The sequence is called a Markov chain, if the conditional distribution of xs given x₀, x₁, …, x_s−1 depends only on the previous state x_s−1 (Markov property): 

MCMC works designing a Markov chain, such that the stationary distribution is the posterior distribution p(x |q). Below we detail the steps of the adaptive Metropolis-Hastings (M-H) MCMC algorithm. 

Step 1: Initialize Markov chains and parameter’s vector x₀. Multiple Markov chains are necessary to perform MCMC. The standard procedure is to run at least three (C=3) parallel chains to assess convergence using a metric, such as the Rˆ diagnostic. In addition, we set the initial starting simulation value x₀ and total number of simulations S for each chain c. 

Step 2: Sample parameter vector x∗ Since the posterior distribution is unknown, we draw samples from another distribution, this is called the proposal distribution [Equation] to generate a candidate parameter vector x∗ (Monte Carlo sampling) based on the previous state of the sampler x_s−1. This work multivariate Gaussian distribution for the proposal distribution r. 

The multivariate normal distribution with dimensionality d (number of model’s parameters) has the following joint density: 

The covariance matrix Σ incorporates the correlation among the model’s parameters. We update this matrix in Step 6. 

Step 3: Compute acceptance ratio. The acceptance ratio is a measure of the likelihood of the sampled parameter vector x∗ to lie on a high-density region of the posterior distribution p(x|q) compared to the previous accepted move x_s−1. 

Step 4: Generate uniform random number. We sample a uniform random number u to decide on whether to accept or reject the proposal x∗. 

Step 5: Update x_s. If u≤ α�≤ 𝛼, we accept the proposal and set the next state equal to the proposal: x_s=x^∗. If u> α, we reject the proposal and set the next state equal to the previous state: x_s=x_s−1. 

Step 6: Update Covariance Matrix Σ

 The acceptance rate ζ is the fraction of accepted moves. The optimal acceptance rate ζ_opt is 0.441 for the dimension d=1 (one parameter), which decreases as the number of dimensions (number of model’s parameters) increases, reaching an asymptotic value of 0.234 for 𝑑→∞ (Brooks et al. 2011). 

The adaptive M-H algorithm updates the covariance matrix Σ of the proposal distribution, based on past simulation draws to ensure that the acceptance rate ζ becomes close to the optimal acceptance ratio ζ_opt. (Haario 2001). 

Step 7: Compute probabilistic rate history. Once we have the posterior distribution of the model’s parameters for each chain x_c=[x₁,x₂,…,x_S], we use these values to compute the probabilistic production history-matches, using Equation 10: 

where x_sc is a posterior simulated draw s of the model’s parameters for a given Markov chain c and qˆ_sc is the rate estimate corresponding to that posterior parameter value. Figure 3 is a flowchart illustrating the steps of the Bayesian variable pressure DCA method. 

Fig. 3 Flowchart illustrating the steps of the probabilistic variable pressure drop DCA method. The method is divided into two stages: I) variable pressure DCA and II) probabilistic variable pressure DCA.

Computational performance. The design of the algorithm meets three requirements: high stability, robustness, and speed. Currently, a typical analysis for an average well requires about 10 sec for the variable pressure DCA method (stage I) and 1 min. for the probabilistic variable-pressure DCA technique (stage II) on a typical laptop computer. We expect much shorter execution times on more sophisticated hardware. Further improvements will shorten the computational time. 

Validation with synthetic case. To validate the workflow, we ran a numerical simulation subject to variable BHP conditions. The numerical simulation represents an idealized hydraulically fractured well that consists of nf  evenly spaced, transverse planar hydraulic fractures that stimulate a rock volume around them. These stimulated rock volumes are assumed to be equal for all hydro-fractures, constituting the stimulated rock volume (SRV) of the reservoir. The fractures have infinite conductivity and thus, fluid that reaches them is instantly produced. In addition, we invoke the following model assumptions: (a) planar fractures of constant length 2xf  and height hf  and two producing faces, (b) flow is one-dimensional and perpendicular to the faces of the hydraulic fractures, (c) natural gas mixture density and compressibility follow the Peng-Robinson 78 EOS (Robinson and Peng 1978) using Van der Waals mixing rule (Sandler 2017), (d) natural gas mixture viscosity is modeled using the Lee-Gonzalez correlation (Lee et al. 1966), (e) Darcy’s law applies, (f) gravity and capillary forces are negligible, (g) constant water saturation at the irreducible level Swi, (h) an isothermal reservoir, (i) constant initial reservoir pressure Pi, and (j) there is a no-flow boundary at half the distance L between adjacent fractures. The natural gas composition for the simulated case is shown in Table 1. This is representative of a well with dry gas composition in the Haynesville shale. 

Table 2 defines the input reservoir properties and completion characteristics for the simulated case under study. 

These original gas in place OGIP_std and the characteristic time τ_i allow to answer two important questions: how much gas is present in the system and how fast we can recover it. 

Using the gas composition in Table 1 and reservoir properties in Table 2, Equation 23 and Equation 24 give the following values for OGIP_std and τ_i for the simulation case

• OGIP_std=10 BSCF
• τ_i=3 years

Figure 4 shows the input BHP and the simulated gas rate histories. We assume a BHP history with periodic step changes (shown by the cyan solid lines). However, for the variable-pressure DCA method, we will purposely include both an incorrect initial reservoir pressure (dashed red line) and noise on the input BHP (shown by the brown diamonds) to effectively simulate uncertainty and error in these two variables. The errors in the initial reservoir pressure and BHP are 15% and up to 45%, respectively. 

Fig. 4. Variable pressure gas rate history (red dots) generated using an initial reservoir pressure of 7,000 psi and a bottomhole flowing pressure history with step changes represented by the cyan lines. We input incorrect initial reservoir pressure and BHP (dashed red line and the brown diamonds, to test the method’s ability to accurately history-match the production and to estimate both the true initial reservoir pressure and the BHP. Errors in the initial reservoir pressure and BHP are 15% and up to 45%.

Decline-curve model. We use the 1-D constant-pressure scaled numerical solution of the diffusivity equation for a real gas mixture of an idealized hydro-fractured well. This model has two fitting parameters: OGIP_std and τ_i: 

Figure 5 illustrates the dimensionless rate solution q_D with dimensionless time t_D on a log-log plot for an initial reservoir pressure of 8,000 psi and different constant BHP conditions: 5,000 psi, 2,000 psi, and 500 psi. For each solution plotted in Fig. 5, we evidence two distinct behaviors. At early dimensionless times (t_D≪1), production follows a constant slope of -1/2 corresponding to the transient flow regime (no-boundary effects). For late dimensionless times (t_D>1), there is a period where the slope is monotonically and continually decreasing; this represents a period of boundary-dominated flow. In addition, as we lower the BHP, the dimensionless rate solution is stretched in dimensionless time.  

Fig. 5. 1-D single-phase real gas mixture dimensionless rate solution qD � � as a function of dimensionless time tD � � for an initial reservoir pressure of 8,000 psi and different constant BHP conditions: 5,000 psi, 2,000 psi, and 500 psi. As the BHP decreases, the dimensionless rate solution is stretched in dimensionless time.

To validate the method (Fig. 3), we proceed as follows. First, we run the numerical decline-curve model from an incorrect initial reservoir pressure of 8,000 psi to a constant BHP of 2,000 psi. In addition, we input the noisy BHP history shown by the brown diamonds in Fig. 5. 

Figure 6 presents the results of the application of the variable pressure drop technique (stage I in Fig. 3) for the synthetic case of Fig. 4. The dashed blue curve is the decline-curve model history-match; the agreement with the gas production history is excellent. The dashed magenta line represents the corrected initial reservoir pressure estimated from the algorithm. This estimate coincides with the true initial reservoir pressure. Finally, the dotted blue curve is the BHP history estimated from the technique, which is in good agreement with the true BHP history. 

Fig. 6. Results of the application of the method for the synthetic case of Fig. 4. The dashed blue lines are the decline-curve model history-match; agreement with the gas production history is excellent. The dashed magenta line represents the corrected initial reservoir pressure estimated from the algorithm. This estimate coincides with the true initial reservoir pressure. Finally, the dotted blue curve is the BHP history, estimated from the technique, which is in good agreement with the true BHP history.

Table 3 compares true and estimated initial reservoir pressure and model’s parameters OGIP_std and τ_i for the synthetic case. The relative errors for these quantities are less than 15%. We attribute these errors to the application of the superposition principle and the use of an incorrect initial reservoir pressure in the 1-D single-phase compressible gas model. The errors are acceptable for the goals of this work. 

Figure 7 shows the probabilistic production history-matches, using the 1-D real gas model, including variable pressure effects (blue curves). The method presents a good coverage of the uncertainty and the variability present in the production rate history (red dots). 

Fig. 7. Probabilistic production history-matches using the 1-D real gas model, including variable pressure effects (blue curves). The method presents a good coverage of the uncertainty and the variability present in the production rate history (red dots).

Figure 8 presents the posterior distribution (top) and trace plots (bottom) of the three parallel Markov chains (C = 3) each with 4,000 simulations (S = 4,000) for the original gas in place (Fig. 8. a.1 and Fig. 8. a.2) and the characteristic time (Figs. 8. b.1 and 8. b.2) parameters. The three chains present similar trace plots; a trace plot shows the observed chain value (x-axis) against the corresponding iteration/simulation number (y-axis), and posterior distributions, illustrating the convergence towards a stationary posterior distribution. A clear, attractive feature of the Bayesian approach over other more traditional approaches, such as Monte Carlo simulation, is that the user does need to specify parameter distribution functions, but the method inherently estimates them. 

Fig. 8. Posterior distribution (top) and trace plots (bottom) of the 3 parallel Markov chains (C = 3), each with 4,000 simulations (S = 4,000) for the original gas in place (Fig. 8.a.1 and Fig 8.a.2) and the characteristic time (Fig. 8.b.1 and Fig. 8.b.2) parameters for well SPE #18. The three chains present similar trace plots and posterior distributions, this fact illustrates the convergence towards a stationary posterior distribution.

SHALE GAS EXAMPLES 

 This section presents the application of the flowchart of Fig. 3 to shale gas wells corresponding to the SPE Data Repository, Data Set #1, Well #17. The first well we analyze is a 9,100-ft lateral-length well from the Haynesville shale. Table 4 presents the reservoir and completion characteristics for Well # 17. 

Figure 9 compares the production history-matches (Fig. 9a) and forecasts (Fig. 9b) of the 1-D real gas model using: pressure-rate-time data (dashed blue curve) and rate-time data (magenta dotted curve). While the pressure-rate-time technique provides an excellent match of the rate history data (red dots), the rate-time analysis does not. In addition, the forecast of the model using pressure-rate-time data is more conservative than the rate-time one (Fig. 9b). 

Fig. 9. Comparison of the production history-matches and forecasts of the 1-D real gas model, using: pressure-rate-time data (dashed blue curve) and rate-time data (magenta dotted curve). While the rate-time-pressure technique provides an excellent history-match to the production data (red dots), rate-time analysis cannot history-match the production of this well (Fig. 9a). In addition, the forecast of the model, using pressure-rate-time data, is more conservative than the rate-time one (Fig. 9b).

Table 5 compares the 1-D real gas model parameters, 20-year EUR, and 20-year recovery factor (RF) values, including variable pressure effects and using only rate-time data of Fig. 9. Including variable pressure effects leads to smaller OGIP_std τ_i, 20-year EUR, and a more realistic 20-year RF for this well. 

Using Equation 23 and Equation 24, we can solve for the fracture half-length x_f and effective gas permeability k_g, respectively. Table 6 compares the 1-D real gas model fracture half-length and gas permeability estimates, including variable pressure effects and using only rate-time data. Using only rate-time data gives unrealistic values for these parameters. 

Figure 10 presents the posterior distribution (top) and trace plots (bottom) of the three parallel Markov chains (C = 3) each with 4,000 simulations (S = 4,000) for the original gas in place (Fig. 10.a.1 and Fig. 10.a.2) and the characteristic time (Fig. 10.b.1 and Fig. 10.b.2) parameters for this well. The three chains present similar trace plots and posterior distributions, illustrating the convergence towards a stationary posterior distribution. 

Fig. 10. Posterior distribution (top) and trace plots (bottom) of the three parallel Markov chains (C = 3) each with 4,000 simulations (S = 4,000) for the original gas in place (Fig. 10.a.1 and Fig. 10.a.2) and the characteristic time (Fig. 10.b.1 and Fig. 10.b.2) parameters for well SPE #17. The three chains present similar trace plots and posterior distributions, this fact illustrates the convergence towards a stationary posterior distribution.

Figure 11 illustrates the probabilistic production history-matches (Fig. 11a) and forecasts (Fig. 11b), using the 1-D real gas model, including variable pressure effects (blue curves). The method presents a good coverage of the uncertainty and the variability present in the production rate history (red dots). 

Fig. 11. Probabilistic production history-matches (Fig. 11a) and forecasts (Fig. 11b), using the 1-D real gas model, including variable pressure effects (blue curves). The method presents a good coverage of the uncertainty and the variability present in the production rate history (red dots).

Figure 12 shows the histograms of the Bayesian posterior distributions of: (a) the 20-year EUR, (b) the fracture half-length, and (c) the effective gas permeability. The histograms also present the deterministic (least-squares value of stage I: L2 fit), the Bayesian P10, P50, and P90 estimates. 

Fig. 12. Histograms of the Bayesian posterior distributions of: (a) the 20-year EUR, (b) the fracture half-length, and (c) the effective gas permeability. The histograms also present the deterministic (least-squares value of stage I: L2 Fit) the Bayesian P10, P50, and P90 estimates.

SPE Data Repository, Data Set 1, Well 18. The second well we analyze is a 9,750-ft lateral-length well from the Haynesville Shale. Table 7 presents the reservoir and completion characteristics for Well 18. 

Figure 13 compares the production history-matches (Fig. 13a) and forecasts (Fig. 13b) of the 1-D real gas model using: pressure-rate-time data (dashed blue curve) and rate-time data (magenta dotted curve). While the pressure-rate-time technique provides an excellent history match to the production data (red dots), rate-time analysis does not. In addition, the forecast of the model using pressure-rate-time data is more conservative, compared to the rate-time one, Fig. 13b. 

Fig. 13. Comparison of the production history-matches and forecasts of the 1-D real gas model, using: pressure-rate-time data (dashed blue curve) and rate-time data (magenta dotted curve). While the rate-time-pressure technique provides an excellent history-match to the production data (red dots), rate-time analysis cannot history-match the production of this well (Fig. 13a). In addition, the forecast of the model, using pressure-rate-time data, is more conservative than the rate-time one (Fig. 13b).

Table 8 compares the 1-D real gas model parameters, 20-year EUR, and 20-year recovery factor (RF) values, including variable pressure effects and using only rate-time data of Fig. 13. Including variable pressure effects leads to smaller OGIP_std, τ_i, 20-year EUR, and a more realistic 20-year RF for this well. 

Using Equation. 23 and Equation 24, we can solve for the fracture half-length x_f and effective gas permeability k_g, respectively. Table 9 compares the 1-D real gas model fracture half-length and gas permeability estimates including variable pressure effects and using only rate-time data. Using only rate-time data gives unrealistic values for these parameters. 

Figure 14 presents the posterior distribution (top) and trace plots (bottom) of the 3 parallel Markov chains (C = 3) each with 4,000 simulations (S = 4,000) for the original gas in place (Fig 14.a.1 and Fig. 14.a.2) and the characteristic time (Fig. 14.b.1 and Fig. 14.b.2) parameters for this well. The three chains present similar trace plots and posterior distributions, illustrating the convergence towards a stationary distribution of the model’s parameters. 

Fig. 14. Posterior distribution (top) and trace plots (bottom) of the 3 parallel Markov chains (C = 3) each with 4,000 simulations (S = 4,000) for the original gas in place (Fig. 14.a.1 and Fig. 14.a.2) and the characteristic time (Fig. 14.b.1 and Fig. 14.b.2) parameters for well SPE #18. The three chains present similar trace plots and posterior distributions, this fact illustrates the convergence towards a stationary posterior distribution.

Figure 15 illustrates the probabilistic production history-matches (Fig. 15a) and forecasts (Fig. 15b) using the 1-D real gas model, including variable pressure effects (blue curves). The method presents a good coverage of the uncertainty and the variability present in the production rate history (red dots). 

Fig. 15. Probabilistic production history-matches (Fig. 15a) and forecasts (Fig. 15b), using the 1-D real gas model, including variable pressure effects (blue curves). The method presents a good coverage of the uncertainty and the variability present in the production rate history (red dots).

Figure 16 shows the histograms of the Bayesian posterior distributions of: (a) the 20-year EUR, (b) the fracture half-length, and (c) the effective gas permeability. The histograms also present the deterministic (least-squares value of stage I: L2 Fit) and the Bayesian P10, P50, and P90 estimates. 

Fig. 16. Histograms of the Bayesian posterior distributions of: (a) the 20-year EUR, (b) the fracture half-length, and (c) the effective gas permeability. The histograms also present the deterministic (least-squares value of stage I: L2 Fit) and the Bayesian P10, P50 and P90 estimates.

CONCLUSIONS

 This article is based on the second paper that extends the deterministic variable-pressure DCA method introduced in the first paper (Ruiz Maraggi et al. 2023) of a two-paper set to: (a) treat compressible gases and (b) include a Bayesian probabilistic calculation. Both papers include a novel method to correct suspected errors in the initial reservoir and bottomhole pressure measurements. The technique history-matches and forecasts production of shale gas wells subject to variable BHP conditions more accurately than traditional DCA. 

The conclusions of this work are the following statements. 

For the synthetic example, the variable-pressure DCA method provides an excellent history-match to the gas production history, despite the presence of errors in both the initial reservoir pressure and the BHP history (Fig. 6). In addition, the method estimates the true model’s parameters and initial reservoir pressure within 15% error (Table 3). 

The proposed method for correcting initial pressure and BHP measurements (Fig. 3) gives a way to infer possible corrections to errors present in the initial reservoir pressure and BHP history for shale gas wells and to estimate the true or actual initial pressure and BHP. 

For the shale gas examples, we show that including variable BHP conditions into the decline-curve model led to more accurate production history-matches and more realistic model estimates ([Equation], [Equation], 20-year EUR, 20-year RF, fracture half-length, and gas permeability), compared to the analysis of decline-curve models using only rate time data (see Fig. 9 and Fig. 13 and Tables 5, 6, 8 and 9). 

The adaptive M-H MCMC algorithm provides a fast convergence of the Markov chain (Figs. 8, 10, and 14). This technique allows us to quantify the uncertainty in the model’s parameters and future production estimates (Fig. 12, Fig. 13 and Fig. 16 and Fig. 17). 

ACKNOWLEDGEMENTS 

This work was partially funded by the State of Texas Advanced Resource Recovery Program (STARR) at the Bureau of Economic Geology. Larry W. Lake holds the Sharon and Shahid Ullah Chair at the Center for Subsurface Energy and the Environment at The University of Texas at Austin. Mark P. Walsh is a Research Fellow at the Hildebrand Department of Petroleum and Geosystems Engineering at The University of Texas at Austin and an oil and gas consultant in Austin, Texas. Frank R. Male is supported through promotion/tenure funds provided internally by The Pennsylvania State University.  

This article is based on URTeC paper 3856826, “Bayesian Variable Pressure Decline-Curve Analysis for Shale Gas Wells,” and is reprinted with permission from the Unconventional Resources Technology Conference, whose permission is required for further use. The paper was presented at the SPE/AAPG/SEG Unconventional Resources Technology Conference, Denver, Colo., June 13–15 2023. The original paper was authorized by the director of the Bureau of Economic Geology. 

About the Authors

Dr. Leopoldo Ruiz Maraggi

Bureau of Economic Geology, University of Texas at Austin

Dr. Leopoldo Ruiz Maraggi is a reservoir engineer within the State of Texas Advanced Resource Recovery program at the Bureau of Economic Geology. His research interests include conventional and unconventional oil and gas resources, subsurface hydrogen and CO2 storage. Dr. Maraggi holds a BS degree in chemical engineering from University of Buenos Aires (Argentina), as well as an MS degree and a PhD in petroleum engineering from The University of Texas at Austin.

Mark P. Walsh

Hildebrant Department of Petroleum and Geosystems Engineering at the University of Texas at Austin

Mark P. Walsh is a consulting reservoir engineer in Austin, Texas, and a research fellow in the Department of Petroleum and Geosystems Engineering at the University of Texas at Austin. He has over 35 years of experience, starting his career with Amoco Production and later serving as a professor of petroleum engineering at Texas A&M University. Dr. Walsh also has been principal advisor at Gaffney, Cline and Associates, chief reservoir engineer at Wapiti Energy and project manager at the Bureau of Economic Geology in Austin. He has a BS degree in chemical engineering from the University of Illinois at Urbana, an MS degree in chemical engineering from the University of Texas at Austin and a PhD in petroleum engineering from the University of Texas at Austin.

Larry W. Lake

Hildebrant Department of Petroleum and Geosystems Engineering at the University of Texas at Austin

Larry W. Lake is a professor in the Hildebrand Department of Petroleum and Geosystems Engineering at The University of Texas at Austin, where he holds the Shahid and Sharon Ullah Chair. He holds a BSE degree and a PhD in chemical engineering from Arizona State University and Rice University, respectively. He is the author or co-author of more than 150 technical papers, four textbooks and the editor of three bound volumes. Dr. Lake has served on the Board of Directors for the Society of Petroleum Engineers, won the 1996 Anthony F. Lucas Gold Medal of the AIME, been awarded the DeGolyer Distinguished Service Award in 2002, and has been a member of the U.S. National Academy of Engineers since 1997. He won the SPE/DOE IOR Pioneer Award in 2000. He was named a Distinguished Graduate of UT in 2022. He has been at the University of Texas since 1978.

Frank R. Male

Pennsylvania State University

Frank R. Male is an assistant research professor of Energy and Mineral Engineering at Penn State University, University Park. His research focuses include interpretable machine learning, reservoir engineering, and risk. He was awarded a PhD in physics from UT Austin in 2015. In 2009, he received a BS degree in physics and a BA degree in political science from Kansas State University.