Balance de Materia: Ajuste de Np y Verificación de Índices de Empuje

Parte 3: Simulación de la producción acumulada mediante regresión no lineal y validación del mecanismo de empuje

Yacimientos
MBE
Python
R

Tercera parte del balance de materia para yacimientos de aceite. Simulamos Np a partir de un valor de N obtenido por Havlena-Odeh o Campbell, lo ajustamos con regresión no lineal (Levenberg-Marquardt) para mejorar el ajuste al historial de producción, y verificamos el mecanismo de empuje dominante con los índices de empuje.

Author

Rigoberto Chandomi

Published

August 2, 2026

Introducción

En la Parte 2 usamos el método de Havlena-Odeh para estimar el OOIP (N) a partir de la pendiente de F vs Et. Ese valor de N es un buen punto de partida, pero no siempre reproduce con precisión la producción acumulada histórica (\(N_p\)).

En esta tercera parte damos un paso adicional:

  1. Simular \(N_p\) con el N obtenido previamente (ej. de Havlena-Odeh o Campbell)
  2. Comparar contra el \(N_p\) real y detectar desviaciones
  3. Ajustar N con regresión no lineal (Levenberg-Marquardt) para minimizar el error
  4. Recalcular los índices de empuje con el N ajustado y verificar el mecanismo dominante
¿Por qué ajustar Np?

El gráfico de Havlena-Odeh es sensible a la dispersión de los primeros puntos y a errores de medición. Simular Np con el N estimado y compararlo contra el histórico real es una prueba de consistencia adicional — si Np calculado se aleja del observado, especialmente en etapas avanzadas de depleción, el N asumido probablemente necesita ajuste.

La Ecuación de Np

Para un yacimiento sin acuífero (\(W_e = 0\)) y sin casquete de gas (\(m = 0\)), despejando \(N_p\) de la MBE:

\[N_p = \frac{N B_{oi}}{B_o + (R_p - R_s)B_g} \left[ \frac{(B_o - B_{oi}) + (R_{si} - R_s)B_g}{B_{oi}} + \frac{c_w S_{wc} + c_f}{1 - S_{wc}} \Delta p \right]\]

Esta ecuación nos permite simular la producción acumulada esperada para cualquier valor de N, y compararla contra el historial real.

Datos del Ejemplo

Usamos un historial de producción con presión promedio, Np, Gp y propiedades PVT por fecha:

Ver código 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv("MB_ex.csv", parse_dates=["Date"], dayfirst=True)

# === PROPIEDADES DE REFERENCIA (primer registro) ===
Boi = df["Bo"].iloc[0]
Rsi = df["Rs"].iloc[0]
Pi  = df["Pavg"].iloc[0]

# === PARÁMETROS DE ROCA Y AGUA ===
Swc = 0.15
cf  = 3.5e-6
cw  = 7.0e-6

print(f"Registros: {len(df)}")
Registros: 31
Ver código 
print(f"Boi = {Boi}, Rsi = {Rsi}, Pi = {Pi} psia")
Boi = 1.376, Rsi = 0.8, Pi = 5200.0 psia
Ver código 
print(f"Rango de fechas: {df['Date'].min().date()} a {df['Date'].max().date()}")
Rango de fechas: 2003-01-01 a 2010-06-23
Ver código 
library(dplyr)
Warning: package 'dplyr' was built under R version 4.3.3

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
Ver código 
BM_data <- read.csv("MB_ex.csv")
BM_data$Date <- as.Date(BM_data$Date, format = "%d/%m/%Y")

Swc <- 0.15
cf  <- 3.5e-6
cw  <- 7.0e-6

Boi <- BM_data$Bo[1]
Rsi <- BM_data$Rs[1]
Pi  <- BM_data$Pavg[1]

cat(sprintf("Registros: %d\n", nrow(BM_data)))
Registros: 31
Ver código 
cat(sprintf("Boi = %.4f, Rsi = %.1f, Pi = %.2f psia\n", Boi, Rsi, Pi))
Boi = 1.3760, Rsi = 0.8, Pi = 5200.00 psia

Paso 1 — Simular Np con el N inicial

Partimos de un valor de N obtenido previamente (Havlena-Odeh o Campbell). En este ejemplo usamos N = 250 MMSTB como estimación inicial:

Ver código 
def np_calc(N, d):
    Et_over_Boi = ((d["Bo"] - Boi) + (Rsi - d["Rs"]) * d["Bg"]) / Boi
    Efw_term = (Pi - d["Pavg"]) * (cw * Swc + cf) / (1 - Swc)
    return (N * Boi / (d["Bo"] + d["Bg"] * (d["Rp"] - d["Rs"]))) * (Et_over_Boi + Efw_term)

N_inicial = 250
df["Np_cal"] = np_calc(N_inicial, df)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(df["Date"], df["Np"], 'o', ms=8, color='#2d8a4e', label='Np observado', zorder=5)
ax.plot(df["Date"], df["Np_cal"], '-', color='#c0392b', lw=2, label=f'Np simulado (N = {N_inicial} MMSTB)')
ax.set_xlabel('Fecha', fontsize=12)
ax.set_ylabel('Np (MMSTB)', fontsize=12)
ax.set_title('Ajuste inicial de Np', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.xticks(rotation=45)
(array([11688., 12053., 12418., 12784., 13149., 13514., 13879., 14245.,
       14610.]), [Text(11688.0, 0, '2002'), Text(12053.0, 0, '2003'), Text(12418.0, 0, '2004'), Text(12784.0, 0, '2005'), Text(13149.0, 0, '2006'), Text(13514.0, 0, '2007'), Text(13879.0, 0, '2008'), Text(14245.0, 0, '2009'), Text(14610.0, 0, '2010')])
Ver código 
plt.tight_layout()
plt.show()
Figure 1: Np observado vs Np simulado con N = 250 MMSTB. Las desviaciones al final del historial sugieren que N puede ajustarse.
Ver código 
library(ggplot2)

np_calc <- function(N, d) {
  Et_over_Boi <- ((d$Bo - Boi) + (Rsi - d$Rs) * d$Bg) / Boi
  Efw_term <- (Pi - d$Pavg) * (cw * Swc + cf) / (1 - Swc)
  (N * Boi / (d$Bo + d$Bg * (d$Rp - d$Rs))) * (Et_over_Boi + Efw_term)
}

N_inicial <- 250
BM_data$Np_cal <- np_calc(N_inicial, BM_data)

ggplot(BM_data) +
  geom_point(aes(x = Date, y = Np), color = "#2d8a4e", size = 3) +
  geom_line(aes(x = Date, y = Np_cal), color = "#c0392b", linewidth = 1) +
  labs(x = "Fecha", y = "Np (MMSTB)",
       title = "Ajuste inicial de Np",
       subtitle = sprintf("N = %d MMSTB", N_inicial)) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14),
        axis.text.x = element_text(angle = 45, hjust = 1))
Figure 2: Np observado vs Np simulado con N = 250 MMSTB.
Interpretación

Si el ajuste inicial se separa del histórico real — típicamente hacia el final, cuando el yacimiento ya ha depletado más presión — el N usado necesita corrección. Un ajuste visual no es suficiente para cuantificar el error; el siguiente paso usa regresión no lineal.

Paso 2 — Ajuste no lineal de N (Levenberg-Marquardt)

Usamos el algoritmo de Levenberg-Marquardt para encontrar el valor de N que minimiza la diferencia entre \(N_p\) observado y calculado:

Ver código 
from scipy.optimize import least_squares

def residuos(params, d):
    N = params[0]
    return d["Np"].values - np_calc(N, d).values

resultado = least_squares(residuos, x0=[N_inicial], args=(df,))
N_fit = resultado.x[0]

df["Np_cal"] = np_calc(N_fit, df)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(df["Date"], df["Np"], 'o', ms=8, color='#2d8a4e', label='Np observado', zorder=5)
ax.plot(df["Date"], df["Np_cal"], '-', color='#c0392b', lw=2, label=f'Np simulado (N = {N_fit:.1f} MMSTB)')
ax.set_xlabel('Fecha', fontsize=12)
ax.set_ylabel('Np (MMSTB)', fontsize=12)
ax.set_title('Np ajustado con regresión no lineal', fontsize=14, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.xticks(rotation=45)
(array([11688., 12053., 12418., 12784., 13149., 13514., 13879., 14245.,
       14610.]), [Text(11688.0, 0, '2002'), Text(12053.0, 0, '2003'), Text(12418.0, 0, '2004'), Text(12784.0, 0, '2005'), Text(13149.0, 0, '2006'), Text(13514.0, 0, '2007'), Text(13879.0, 0, '2008'), Text(14245.0, 0, '2009'), Text(14610.0, 0, '2010')])
Ver código 
plt.tight_layout()
plt.show()
Figure 3: Np observado vs Np simulado con el N ajustado por mínimos cuadrados no lineales.
Ver código 
print(f"\n✅ N ajustado = {N_fit:.2f} MMSTB (inicial: {N_inicial} MMSTB)")

✅ N ajustado = 253.26 MMSTB (inicial: 250 MMSTB)
Ver código 
library(minpack.lm)

residFun <- function(parS, observed, d) {
  observed - np_calc(parS$N, d)
}

parStart <- list(N = N_inicial)

fit_model <- nls.lm(par = parStart, fn = residFun, observed = BM_data$Np,
                     d = BM_data, control = nls.lm.control(nprint = 0))

N_fit <- coef(fit_model)[["N"]]
BM_data$Np_cal <- np_calc(N_fit, BM_data)

ggplot(BM_data) +
  geom_point(aes(x = Date, y = Np), color = "#2d8a4e", size = 3) +
  geom_line(aes(x = Date, y = Np_cal), color = "#c0392b", linewidth = 1) +
  labs(x = "Fecha", y = "Np (MMSTB)",
       title = "Np ajustado con regresión no lineal",
       subtitle = sprintf("N = %.1f MMSTB", N_fit)) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14),
        axis.text.x = element_text(angle = 45, hjust = 1))

cat(sprintf("N ajustado = %.2f MMSTB (inicial: %d MMSTB)\n", N_fit, N_inicial))
N ajustado = 253.26 MMSTB (inicial: 250 MMSTB)
Figure 4: Np observado vs Np simulado con el N ajustado por mínimos cuadrados no lineales.

Paso 3 — Índices de Empuje con el N ajustado

Con el N ajustado, recalculamos los índices de empuje:

\[DDI = \frac{N\left[(B_o - B_{oi}) + (R_{si} - R_s)B_g\right]}{N_p\left[B_o + (R_p - R_s)B_g\right]}\] \[CDI = \frac{NB_{oi}(1+m) \left( \frac{c_wS_{wc}+c_f}{1-S_{wc}} \right) \Delta p}{N_p(B_o+(R_{si}-R_s)B_g)}\]

\[CDI = 1 - DDI\]

El punto inicial (\(N_p = 0\)) se excluye por indeterminación matemática (división entre cero).

Ver código 
mask = df["Np"] > 0
df_di = df[mask].copy()

df_di["DDI"] = N_fit * ((df_di["Bo"] - Boi) + (Rsi - df_di["Rs"]) * df_di["Bg"]) / \
               (df_di["Np"] * (df_di["Bo"] + df_di["Bg"] * (df_di["Rp"] - df_di["Rs"])))
df_di["CDI"] = 1 - df_di["DDI"]

fig, ax = plt.subplots(figsize=(10, 6))
ax.stackplot(df_di["Date"], df_di["DDI"], df_di["CDI"],
             colors=['#2d8a4e', '#2980b9'],
             labels=['DDI (Expansión aceite + gas)', 'CDI (Compactación + agua connata)'],
             alpha=0.85)

ax.axhline(y=1.0, color='#c0392b', ls='--', lw=1, alpha=0.7)
ax.set_xlabel('Fecha', fontsize=12)
ax.set_ylabel('Índice de empuje', fontsize=12)
ax.set_title('Índices de Empuje (Drive Indices)', fontsize=14, fontweight='bold')
ax.legend(fontsize=10, loc='lower left')
ax.set_ylim(0, 1.0)
(0.0, 1.0)
Ver código 
ax.set_xlim(df_di["Date"].min(), df_di["Date"].max())
(np.float64(12144.0), np.float64(14783.0))
Ver código 
ax.grid(True, alpha=0.2, axis='y')
plt.xticks(rotation=45)
(array([12053., 12418., 12784., 13149., 13514., 13879., 14245., 14610.]), [Text(12053.0, 0, '2003'), Text(12418.0, 0, '2004'), Text(12784.0, 0, '2005'), Text(13149.0, 0, '2006'), Text(13514.0, 0, '2007'), Text(13879.0, 0, '2008'), Text(14245.0, 0, '2009'), Text(14610.0, 0, '2010')])
Ver código 
plt.tight_layout()
plt.show()
Figure 5: Índices de empuje vs tiempo, graficados como área apilada entre 0 y 1.
Ver código 
library(tidyr)

mask <- BM_data$Np > 0
df_di <- BM_data[mask, ]

df_di$DDI <- N_fit * ((df_di$Bo - Boi) + (Rsi - df_di$Rs) * df_di$Bg) /
             (df_di$Np * (df_di$Bo + df_di$Bg * (df_di$Rp - df_di$Rs)))
df_di$CDI <- 1 - df_di$DDI

df_long <- pivot_longer(df_di[, c("Date", "DDI", "CDI")], cols = c(DDI, CDI),
                         names_to = "Index", values_to = "Value")
df_long$Index <- factor(df_long$Index, levels = c("DDI", "CDI"))

ggplot(df_long, aes(x = Date, y = Value, fill = Index)) +
  geom_area(position = "stack", alpha = 0.85) +
  scale_fill_manual(values = c("DDI" = "#2d8a4e", "CDI" = "#2980b9"),
                     labels = c("DDI (Expansión aceite + gas)",
                                "CDI (Compactación + agua connata)")) +
  geom_hline(yintercept = 1.0, color = "#c0392b", linetype = "dashed") +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
  labs(x = "Fecha", y = "Índice de empuje",
       title = "Índices de Empuje (Drive Indices)", fill = "Mecanismo") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 14),
        axis.text.x = element_text(angle = 45, hjust = 1))
Figure 6: Índices de empuje vs tiempo, graficados como área apilada entre 0 y 1.

Interpretación de los Resultados

¿Qué significan los índices?

Interpretación de los drive indices.
Índice Mecanismo Típico cuando…
DDI alto (>0.8) Expansión de aceite y gas liberado Yacimiento volumétrico, sin soporte de presión
CDI significativo (>0.1) Compactación de roca y expansión de agua connata Yacimientos profundos, alta compresibilidad
WDI alto (>0.5) Invasión de agua del acuífero Acuífero activo
SDI alto (>0.3) Expansión del casquete de gas Casquete de gas significativo

En este ejemplo, DDI domina durante todo el historial de producción (~0.90-0.93), lo que confirma que el yacimiento produce principalmente por expansión del aceite y del gas liberado (depletion drive). El CDI contribuye entre un 7-10%, algo mayor que en el caso volumétrico puro de la Parte 2, consistente con la mayor caída de presión observada en este ejemplo.

Verificación

La suma DDI + CDI debe ser ≈ 1.0 en cada punto del historial. El área apilada facilita esta verificación visualmente: si la línea superior se aleja de 1.0, hay un error en el N ajustado o en el modelo asumido (sin acuífero, sin casquete de gas).

Resumen

Resultado Valor
N inicial (Havlena-Odeh / Campbell) 250 MMSTB
N ajustado por regresión no lineal Valor obtenido por Levenberg-Marquardt
Mecanismo dominante Depletion drive (DDI ≈ 0.90-0.93)
Contribución de compactación CDI ≈ 0.07-0.10
Modelo validado Ajuste consistente de Np y DDI+CDI ≈ 1.0

Referencias

  • Havlena, D. & Odeh, A.S. (1963). The Material Balance as an Equation of a Straight Line. JPT, 15(8).
  • Ahmed, T. (2019). Reservoir Engineering Handbook, 5th ed. Gulf Professional Publishing. Cap. 12-13.
  • Sanni, M. (2019). Petroleum Engineering: Principles, Calculations, and Workflows. Wiley/AGU.
  • Dake, L.P. (1978). Fundamentals of Reservoir Engineering. Elsevier. Cap. 3.
🛒
Propiedades De Los Fluidos Del Yacimiento

Disponible en Mercado Libre

Ver oferta →