Ajuste de Curvas de Declinación con Python y R

Modelos de Arps: Exponencial, Hiperbólico y Armónico con regresión no lineal

Ajusta datos de producción reales usando los 3 modelos clásicos de Arps (exponencial, hiperbólico y armónico) con regresión no lineal. Compara el ajuste y el EUR de cada modelo. Código disponible en Python y R.

Curvas de Declinación de Arps

El análisis de curvas de declinación (DCA, Decline Curve Analysis) es una de las herramientas más utilizadas para pronosticar producción y estimar reservas remanentes. Arps (1945) propuso tres modelos empíricos basados en la ecuación general:

\[q(t) = \frac{q_i}{(1 + bD_it)^{1/b}}\]

Donde:

• \(q_i\) = gasto inicial (tasa al inicio de la declinación)
• \(D_i\) = tasa de declinación inicial
• \(b\) = exponente de declinación (\(0 \leq b \leq 1\))

Los tres casos particulares son:

Modelo	Valor de b	Ecuación
Exponencial	\(b = 0\)	\(q(t) = q_i \cdot e^{-D_i t}\)
Hiperbólico	\(0 < b < 1\)	\(q(t) = q_i / (1 + bD_it)^{1/b}\)
Armónico	\(b = 1\)	\(q(t) = q_i / (1 + D_it)\)

Visita el canal @rigopetrodata

Tutoriales de R, Python y Excel para ingeniería petrolera

Ver canal

: Modelos de declinación de Arps. {.striped}

¿Cuál modelo usar?

El exponencial es el más conservador (declina más rápido), el armónico el más optimista (declina más lento), y el hiperbólico es intermedio. En la práctica, se ajustan los tres y se comparan para tomar la mejor decisión.

Los Datos

Trabajamos con datos de producción de gas de un pozo real con 480 días de historial:

Python

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

data = pd.read_csv('REH_18.csv')
DAYS = data['DAYS'].values
RATE = data['RATE'].values

print(f"Registros: {len(data)}")

Registros: 480

print(f"Rango de días: {DAYS.min()} - {DAYS.max()}")

Rango de días: 1 - 480

print(f"Gasto máximo: {RATE.max():.2f}")

Gasto máximo: 6.00

print(f"Gasto final: {RATE[-1]:.2f}")

Gasto final: 0.51

print(data.head(10))

   DAYS  RATE
0     1   2.3
1     2   3.0
2     3   3.2
3     4   3.0
4     5   3.2
5     6   4.2
6     7   4.3
7     8   4.3
8     9   4.3
9    10   4.3

R

data <- read.csv('REH_18.csv')
DAYS <- data$DAYS
RATE <- data$RATE

cat(sprintf("Registros: %d\n", nrow(data)))

Registros: 480

cat(sprintf("Rango de días: %d - %d\n", min(DAYS), max(DAYS)))

Rango de días: 1 - 480

cat(sprintf("Gasto máximo: %.2f\n", max(RATE)))

Gasto máximo: 6.00

cat(sprintf("Gasto final: %.2f\n", tail(RATE, 1)))

Gasto final: 0.51

head(data, 10)

   DAYS RATE
1     1  2.3
2     2  3.0
3     3  3.2
4     4  3.0
5     5  3.2
6     6  4.2
7     7  4.3
8     8  4.3
9     9  4.3
10   10  4.3

El pozo muestra un periodo inicial de incremento (primeros ~17 días, posiblemente limpieza o apertura gradual) seguido de una declinación sostenida. Para el ajuste usaremos todos los datos, aunque en la práctica se podría filtrar el periodo transitorio inicial.

Visualización de los Datos

Python

fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(DAYS, RATE, s=8, alpha=0.6, color='#2d8a4e', label='Datos')
ax.set_xlabel('Tiempo (días)', fontsize=12)
ax.set_ylabel('Qg', fontsize=12)
ax.set_title('Datos de Producción', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Figure 1: Datos de producción de gas. Se observa un periodo inicial de incremento seguido de declinación.

R

library(ggplot2)

ggplot(data, aes(x = DAYS, y = RATE)) +
  geom_point(color = "#2d8a4e", size = 1.5, alpha = 0.6) +
  labs(x = "Tiempo (días)", y = "Qg",
       title = "Datos de Producción") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14))

Figure 2: Datos de producción de gas.

Ajuste del Modelo Hiperbólico

El modelo hiperbólico es el más general de los tres y tiene 3 parámetros a estimar: \(q_i\), \(b\) y \(D_i\). Usamos regresión no lineal para encontrar los valores que minimizan el error cuadrático.

\[q(t) = \frac{q_i}{(1 + bD_it)^{1/b}}\]

Python

En Python usamos curve_fit() del paquete scipy.optimize:

# Definir la función del modelo
def hyperbolic(t, qi, b, Di):
    return qi / ((1 + b * Di * t) ** (1/b))

# Ajustar con curve_fit (regresión no lineal)
coef_hyp, cov_hyp = curve_fit(hyperbolic, DAYS, RATE,
                               p0=[6, 0.5, 0.01],    # valores iniciales
                               maxfev=10000)

qi_h, b_h, Di_h = coef_hyp
print(f"qi = {qi_h:.4f}")

qi = 5.2213

print(f"b  = {b_h:.4f}")

b  = 0.5029

print(f"Di = {Di_h:.6f}")

Di = 0.011202

R

En R usamos nls() (nonlinear least squares):

# Ajustar modelo hiperbólico
fit_hyp <- nls(RATE ~ qi / ((1 + b * Di * DAYS)^(1/b)),
               data = data,
               start = list(qi = 6, b = 0.5, Di = 0.01),
               control = nls.control(maxiter = 1000))

coef_hyp <- coef(fit_hyp)
cat(sprintf("qi = %.4f\n", coef_hyp["qi"]))

qi = 5.2212

cat(sprintf("b  = %.4f\n", coef_hyp["b"]))

b  = 0.5029

cat(sprintf("Di = %.6f\n", coef_hyp["Di"]))

Di = 0.011202

Valores iniciales (p0 / start)

La regresión no lineal requiere valores iniciales razonables para converger. Si el ajuste falla, prueba con diferentes semillas. Una buena estrategia: usar \(q_i\) = gasto máximo observado, \(b\) = 0.5 y \(D_i\) = 0.01.

Ajuste del Modelo Exponencial

El exponencial tiene solo 2 parámetros (\(q_i\) y \(D_i\)) y asume declinación constante:

\[q(t) = q_i \cdot e^{-D_i t}\]

Python

def exponential(t, qi, Di):
    return qi * np.exp(-Di * t)

coef_exp, cov_exp = curve_fit(exponential, DAYS, RATE, p0=[6, 0.005])
qi_e, Di_e = coef_exp
print(f"qi = {qi_e:.4f}")

qi = 4.6442

print(f"Di = {Di_e:.6f}")

Di = 0.006850

R

fit_exp <- nls(RATE ~ qi * exp(-Di * DAYS),
               data = data,
               start = list(qi = 6, Di = 0.005))

coef_exp <- coef(fit_exp)
cat(sprintf("qi = %.4f\n", coef_exp["qi"]))

qi = 4.6442

cat(sprintf("Di = %.6f\n", coef_exp["Di"]))

Di = 0.006850

Ajuste del Modelo Armónico

El armónico también tiene 2 parámetros y es el caso particular donde \(b = 1\):

\[q(t) = \frac{q_i}{1 + D_it}\]

Python

def harmonic(t, qi, Di):
    return qi / (1 + Di * t)

coef_har, cov_har = curve_fit(harmonic, DAYS, RATE, p0=[6, 0.01])
qi_a, Di_a = coef_har
print(f"qi = {qi_a:.4f}")

qi = 5.6238

print(f"Di = {Di_a:.6f}")

Di = 0.017930

R

fit_har <- nls(RATE ~ qi / (1 + Di * DAYS),
               data = data,
               start = list(qi = 6, Di = 0.01))

coef_har <- coef(fit_har)
cat(sprintf("qi = %.4f\n", coef_har["qi"]))

qi = 5.6239

cat(sprintf("Di = %.6f\n", coef_har["Di"]))

Di = 0.017931

Comparación Visual de los 3 Modelos

Python

# Predicciones
hyp_pred = hyperbolic(DAYS, *coef_hyp)
exp_pred = exponential(DAYS, *coef_exp)
har_pred = harmonic(DAYS, *coef_har)

fig, ax = plt.subplots(figsize=(12, 7))
ax.scatter(DAYS, RATE, s=8, alpha=0.4, color='gray', label='Datos', zorder=1)
ax.plot(DAYS, hyp_pred, '-', lw=2.5, color='#2d8a4e', label='Hiperbólico', zorder=3)
ax.plot(DAYS, exp_pred, '--', lw=2, color='#c0392b', label='Exponencial', zorder=2)
ax.plot(DAYS, har_pred, '-.', lw=2, color='#2980b9', label='Armónico', zorder=2)

ax.set_xlabel('Tiempo (días)', fontsize=12)
ax.set_ylabel('Qg', fontsize=12)
ax.set_title('Ajuste de Curvas de Declinación — Modelos de Arps', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Figure 3: Comparación de los 3 modelos de Arps ajustados a los datos de producción.

R

# Predicciones
data$hyp_pred <- predict(fit_hyp)
data$exp_pred <- predict(fit_exp)
data$har_pred <- predict(fit_har)

ggplot(data, aes(x = DAYS)) +
  geom_point(aes(y = RATE), color = "gray", size = 1, alpha = 0.4) +
  geom_line(aes(y = hyp_pred, color = "Hiperbólico"), linewidth = 1.2) +
  geom_line(aes(y = exp_pred, color = "Exponencial"), linewidth = 1, linetype = "dashed") +
  geom_line(aes(y = har_pred, color = "Armónico"), linewidth = 1, linetype = "dotdash") +
  scale_color_manual(values = c("Hiperbólico" = "#2d8a4e",
                                 "Exponencial" = "#c0392b",
                                 "Armónico" = "#2980b9")) +
  labs(x = "Tiempo (días)", y = "Qg", color = "Modelo",
       title = "Ajuste de Curvas de Declinación — Modelos de Arps") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14))

Figure 4: Comparación de los 3 modelos de Arps.

Pronóstico y EUR

La verdadera utilidad del DCA es pronosticar la producción futura y estimar las reservas recuperables (EUR). Extendemos los modelos hasta 1,000 días (520 días más allá de los datos):

Python

# Pronóstico extendido
t_forecast = np.arange(1, 1001)
hyp_fc = hyperbolic(t_forecast, *coef_hyp)
exp_fc = exponential(t_forecast, *coef_exp)
har_fc = harmonic(t_forecast, *coef_har)

# EUR (producción acumulada)
EUR_hyp = np.sum(hyp_fc)
EUR_exp = np.sum(exp_fc)
EUR_har = np.sum(har_fc)

fig, ax = plt.subplots(figsize=(12, 7))

# Zona de pronóstico
ax.axvspan(480, 1000, alpha=0.08, color='orange', label='Pronóstico')

ax.scatter(DAYS, RATE, s=8, alpha=0.4, color='gray', zorder=1)
ax.plot(t_forecast, hyp_fc, '-', lw=2.5, color='#2d8a4e',
        label=f'Hiperbólico (EUR={EUR_hyp:.0f})')
ax.plot(t_forecast, exp_fc, '--', lw=2, color='#c0392b',
        label=f'Exponencial (EUR={EUR_exp:.0f})')
ax.plot(t_forecast, har_fc, '-.', lw=2, color='#2980b9',
        label=f'Armónico (EUR={EUR_har:.0f})')

ax.axvline(x=480, color='orange', ls=':', alpha=0.5)
ax.set_xlabel('Tiempo (días)', fontsize=12)
ax.set_ylabel('Qg', fontsize=12)
ax.set_title('Pronóstico de Producción y EUR', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 1050)

(0.0, 1050.0)

plt.tight_layout()
plt.show()

Figure 5: Pronóstico de producción a 1,000 días con los 3 modelos. La zona sombreada es la extrapolación.

R

t_fc <- 1:1000
qi_h <- coef(fit_hyp)["qi"]; b_h <- coef(fit_hyp)["b"]; Di_h <- coef(fit_hyp)["Di"]
qi_e <- coef(fit_exp)["qi"]; Di_e <- coef(fit_exp)["Di"]
qi_a <- coef(fit_har)["qi"]; Di_a <- coef(fit_har)["Di"]

fc <- data.frame(
  t = t_fc,
  hyp = qi_h / ((1 + b_h * Di_h * t_fc)^(1/b_h)),
  exp = qi_e * exp(-Di_e * t_fc),
  har = qi_a / (1 + Di_a * t_fc)
)

EUR_hyp <- sum(fc$hyp)
EUR_exp <- sum(fc$exp)
EUR_har <- sum(fc$har)

cat(sprintf("EUR Hiperbólico: %.0f\n", EUR_hyp))

EUR Hiperbólico: 791

cat(sprintf("EUR Exponencial: %.0f\n", EUR_exp))

EUR Exponencial: 675

cat(sprintf("EUR Armónico:    %.0f\n", EUR_har))

EUR Armónico:    920

library(tidyr)
fc_long <- pivot_longer(fc, cols = c(hyp, exp, har),
                         names_to = "Modelo", values_to = "Qg")
fc_long$Modelo <- factor(fc_long$Modelo,
                          levels = c("hyp", "exp", "har"),
                          labels = c("Hiperbólico", "Exponencial", "Armónico"))

ggplot() +
  annotate("rect", xmin = 480, xmax = 1000, ymin = -Inf, ymax = Inf,
           alpha = 0.08, fill = "orange") +
  geom_point(data = data, aes(x = DAYS, y = RATE), color = "gray", size = 1, alpha = 0.4) +
  geom_line(data = fc_long, aes(x = t, y = Qg, color = Modelo, linetype = Modelo),
            linewidth = 1) +
  scale_color_manual(values = c("Hiperbólico" = "#2d8a4e",
                                 "Exponencial" = "#c0392b",
                                 "Armónico" = "#2980b9")) +
  geom_vline(xintercept = 480, color = "orange", linetype = "dotted", alpha = 0.5) +
  labs(x = "Tiempo (días)", y = "Qg",
       title = "Pronóstico de Producción y EUR") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14))

Figure 6: Pronóstico de producción a 1,000 días con los 3 modelos.

Tabla Resumen de Resultados

Python

# R-cuadrado
SS_res_hyp = np.sum((RATE - hyp_pred)**2)
SS_res_exp = np.sum((RATE - exp_pred)**2)
SS_res_har = np.sum((RATE - har_pred)**2)
SS_tot = np.sum((RATE - np.mean(RATE))**2)

R2_hyp = 1 - SS_res_hyp / SS_tot
R2_exp = 1 - SS_res_exp / SS_tot
R2_har = 1 - SS_res_har / SS_tot

print(f"{'Parámetro':<15} {'Exponencial':>14} {'Hiperbólico':>14} {'Armónico':>14}")

Parámetro          Exponencial    Hiperbólico       Armónico

print("-" * 59)

-----------------------------------------------------------

print(f"{'qi':<15} {qi_e:>14.4f} {qi_h:>14.4f} {qi_a:>14.4f}")

qi                      4.6442         5.2213         5.6238

print(f"{'b':<15} {'0 (fijo)':>14} {b_h:>14.4f} {'1 (fijo)':>14}")

b                     0 (fijo)         0.5029       1 (fijo)

print(f"{'Di':<15} {Di_e:>14.6f} {Di_h:>14.6f} {Di_a:>14.6f}")

Di                    0.006850       0.011202       0.017930

print(f"{'R²':<15} {R2_exp:>14.4f} {R2_hyp:>14.4f} {R2_har:>14.4f}")

R²                      0.8350         0.8600         0.8397

print(f"{'EUR (1000 d)':<15} {EUR_exp:>14.0f} {EUR_hyp:>14.0f} {EUR_har:>14.0f}")

EUR (1000 d)               675            791            920

R

# R-cuadrado
SS_tot <- sum((RATE - mean(RATE))^2)
R2_hyp <- 1 - sum(residuals(fit_hyp)^2) / SS_tot
R2_exp <- 1 - sum(residuals(fit_exp)^2) / SS_tot
R2_har <- 1 - sum(residuals(fit_har)^2) / SS_tot

results <- data.frame(
  Parámetro = c("qi", "b", "Di", "R²", "EUR (1000 d)"),
  Exponencial = c(sprintf("%.4f", qi_e), "0 (fijo)",
                  sprintf("%.6f", Di_e), sprintf("%.4f", R2_exp),
                  sprintf("%.0f", EUR_exp)),
  Hiperbólico = c(sprintf("%.4f", qi_h), sprintf("%.4f", b_h),
                  sprintf("%.6f", Di_h), sprintf("%.4f", R2_hyp),
                  sprintf("%.0f", EUR_hyp)),
  Armónico = c(sprintf("%.4f", qi_a), "1 (fijo)",
               sprintf("%.6f", Di_a), sprintf("%.4f", R2_har),
               sprintf("%.0f", EUR_har))
)

knitr::kable(results, caption = "Comparación de modelos de Arps", align = 'lrrr')

Comparación de modelos de Arps

Parámetro	Exponencial	Hiperbólico	Armónico
qi	4.6442	5.2212	5.6239
b	0 (fijo)	0.5029	1 (fijo)
Di	0.006850	0.011202	0.017931
R²	0.8350	0.8600	0.8397
EUR (1000 d)	675	791	920

Observaciones:

• El modelo hiperbólico tiene el mejor R² (0.86) ya que tiene un parámetro más de libertad (\(b\))
• El valor ajustado de \(b = 0.50\) es consistente con yacimientos de gas convencional (típicamente \(b\) entre 0.3 y 0.7)
• El exponencial es el más conservador: predice el menor EUR (675) y declina más rápido
• El armónico es el más optimista: predice el mayor EUR (920) y declina más lento
• La diferencia entre EUR exponencial y armónico es de ~36% — esta incertidumbre es importante para la evaluación económica

Precaución con el pronóstico

Los modelos de Arps son empíricos y la extrapolación más allá de los datos siempre tiene incertidumbre. Mientras más corto sea el historial, mayor la incertidumbre del pronóstico. Para yacimientos no convencionales (shale), los modelos de Arps con \(b > 1\) pueden sobreestimar las reservas — en esos casos se recomienda usar modelos específicos como SEDM o Duong.

Descarga

• 📥 REH_18.csv — Datos de producción (480 días)

Referencias

• Arps, J.J. (1945). Analysis of Decline Curves. Transactions of the AIME, 160(1), 228-247.
• Ahmed, T. (2019). Reservoir Engineering Handbook, 5th ed. Gulf Professional Publishing.
• Sun, H. (2015). Advanced Production Decline Analysis and Application. Gulf Professional Publishing.

---

🛒

Propiedades De Los Fluidos Del Yacimiento

Disponible en Mercado Libre

Ver oferta →

---