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

Parte 2: Simulación de la producción acumulada mediante mínimos cuadrados por el origen y validación del mecanismo de empuje

Yacimientos

MBE

Python

Segunda 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 mínimos cuadrados lineales para comparar contra el historial de producción, y examinamos los índices y el cierre del balance dominante con los índices de empuje.

Author

Rigoberto Chandomi

Published

August 2, 2026

Introducción

En la Parte 1 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 segunda parte damos un paso adicional:

Simular \(N_p\) con el N obtenido previamente (ej. de Havlena-Odeh o Campbell)
Comparar contra el \(N_p\) real y detectar desviaciones
Ajustar N con mínimos cuadrados lineales por el origen, porque Np es lineal en N si PVT, presión y Rp están fijados
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, se deben revisar N y los supuestos. La desviación también puede proceder de PVT, presión promedio, producción o un mecanismo omitido.

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 reconstruye \(N_p\) condicionado a la presión, PVT y \(R_p\) históricos. No es un pronóstico independiente: \(R_p=G_p/N_p\) procede del historial y no es el GOR instantáneo. Este ejemplo no incluye producción de agua; si existiera, habría que incorporar \(W_p B_w\).

Datos del Ejemplo

Usamos el mismo CSV de la Parte 1. Sus unidades son: Np en MMSTB, Gp en \(10^9\) scf, Rs y Rp en Mscf/STB y Bg en RB/Mscf; M=\(10^3\), MM=\(10^6\). Por tanto, \((R_p-R_s)B_g\) es RB/STB. Se conserva el redondeo de las columnas publicado en el CSV. La presión es la promedio del yacimiento, no pwf.

Python
R

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:

Python
R

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 la reconstrucción se separa del histórico, revisa los datos y el modelo antes de atribuir toda la diferencia a N. El siguiente paso cuantifica el ajuste. Minimizar errores de Np no equivale a minimizar errores de F: son objetivos con ponderaciones distintas.

Paso 2 — Ajuste lineal de N

Con las propiedades e historia fijadas, la ecuación tiene la forma \(N_{p,i}=N a_i\). No hace falta Levenberg–Marquardt para esta única incógnita:

\[N_{ajustado} = \frac{\sum_i a_i N_{p,i}}{\sum_i a_i^2}, \qquad a_i=N_{p,calculado}(N=1)\]

Python
R

Ver código

a = np_calc(1.0, df)
validos = np.isfinite(a) & np.isfinite(df["Np"]) & (a > 0)
N_fit = np.dot(a[validos], df.loc[validos, "Np"]) / np.dot(a[validos], a[validos])
df["Np_cal"] = np_calc(N_fit, df)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(df["Date"], df["Np"], 'o', color='#2d8a4e', label='Np observado')
ax.plot(df["Date"], df["Np_cal"], '-', color='#c0392b', label=f'Np reconstruido (N={N_fit:.2f} MMSTB)')
ax.set(xlabel='Fecha', ylabel='Np (MMSTB)', title='Ajuste de N por mínimos cuadrados lineales')
ax.legend()
ax.grid(True, alpha=.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 y reconstruido con mínimos cuadrados lineales por el origen.

Ver código

print(f"N ajustado = {N_fit:.4f} MMSTB")

N ajustado = 253.2646 MMSTB

Ver código

a <- np_calc(1, BM_data)
validos <- is.finite(a) & is.finite(BM_data$Np) & a > 0
N_fit <- sum(a[validos]*BM_data$Np[validos])/sum(a[validos]^2)
BM_data$Np_cal <- np_calc(N_fit, BM_data)

ggplot(BM_data) +
  geom_point(aes(Date, Np), color="#2d8a4e", size=3) +
  geom_line(aes(Date, Np_cal), color="#c0392b", linewidth=1) +
  labs(x="Fecha", y="Np (MMSTB)", title="Ajuste de N por mínimos cuadrados lineales",
       subtitle=sprintf("N = %.2f MMSTB", N_fit)) + theme_minimal()
cat(sprintf("N ajustado = %.4f MMSTB\n", N_fit))

N ajustado = 253.2646 MMSTB

Figure 4: Np observado y reconstruido con mínimos cuadrados lineales por el origen.

Con el CSV publicado, N≈253.26 MMSTB. El ajuste de F frente a Et por el origen de la Parte 1 da aproximadamente 253.76 MMSTB. La diferencia responde al objetivo minimizado; no implica una contradicción entre las ecuaciones.

Paso 3 — Índices de empuje y cierre independiente

Para este caso, sin acuífero, casquete ni producción de agua:

\[F=N_p[B_o+(R_p-R_s)B_g],\quad E_o=(B_o-B_{oi})+(R_{si}-R_s)B_g\] \[E_{fw}=B_{oi}\frac{c_w S_{wc}+c_f}{1-S_{wc}}(P_i-p)\] \[DDI=\frac{N E_o}{F},\qquad CDI=\frac{N E_{fw}}{F}\] \[\varepsilon=\frac{F-N(E_o+E_{fw})}{F},\qquad DDI+CDI=1-\varepsilon\]

No se define CDI como 1−DDI: eso obligaría a cerrar el balance e impediría verificarlo. El estado inicial, F=0, se excluye de estos cocientes. No se recorta el gráfico a 1, para no ocultar una suma mayor que la unidad.

Python
R

Ver código

F = df["Np"]*(df["Bo"]+(df["Rp"]-df["Rs"])*df["Bg"])
Eo = df["Bo"]-Boi+(Rsi-df["Rs"])*df["Bg"]
Efw = Boi*(cw*Swc+cf)/(1-Swc)*(Pi-df["Pavg"])
mask = F > 0
df_di = df.loc[mask, ["Date"]].copy()
df_di["DDI"] = N_fit*Eo[mask]/F[mask]
df_di["CDI"] = N_fit*Efw[mask]/F[mask]
df_di["error_pct"] = 100*(F[mask]-N_fit*(Eo[mask]+Efw[mask]))/F[mask]

fig, (ax, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
ax.stackplot(df_di["Date"], df_di["DDI"], df_di["CDI"],
             colors=['#2d8a4e','#2980b9'], labels=['Aceite + gas disuelto','Roca + agua connata'])
ax.axhline(1, color='#c0392b', ls='--')
ax.set_ylabel('Fracción de F')
ax.legend()
ax2.plot(df_di["Date"], df_di["error_pct"], color='#c0392b')
ax2.axhline(0, color='gray', ls='--')
ax2.set(xlabel='Fecha', ylabel='Error de cierre (%)')
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 5: Aportes independientes y error de cierre del balance.

Ver código

print(df_di[["DDI","CDI","error_pct"]].agg(['min','max']))

          DDI       CDI  error_pct
min  0.609349  0.094443  -4.405631
max  0.903622  0.434708   0.813412

Ver código

library(tidyr)
F <- BM_data$Np*(BM_data$Bo+(BM_data$Rp-BM_data$Rs)*BM_data$Bg)
Eo <- BM_data$Bo-Boi+(Rsi-BM_data$Rs)*BM_data$Bg
Efw <- Boi*(cw*Swc+cf)/(1-Swc)*(Pi-BM_data$Pavg)
mask <- F > 0
df_di <- BM_data[mask, "Date", drop=FALSE]
df_di$DDI <- N_fit*Eo[mask]/F[mask]
df_di$CDI <- N_fit*Efw[mask]/F[mask]
df_di$error_pct <- 100*(F[mask]-N_fit*(Eo[mask]+Efw[mask]))/F[mask]
df_long <- pivot_longer(df_di[,c("Date","DDI","CDI")], cols=c(DDI,CDI),
                        names_to="Index", values_to="Value")
ggplot(df_long, aes(Date,Value,fill=Index)) + geom_area(position="stack",alpha=.85) +
  geom_hline(yintercept=1,linetype="dashed",color="#c0392b") +
  labs(x="Fecha",y="Fracción de F",title="Aportes calculados por separado") + theme_minimal()
ggplot(df_di,aes(Date,error_pct)) + geom_line(color="#c0392b") +
  geom_hline(yintercept=0,linetype="dashed") +
  labs(x="Fecha",y="Error de cierre (%)") + theme_minimal()
print(sapply(df_di[,c("DDI","CDI","error_pct")],range))

           DDI       CDI  error_pct
[1,] 0.6093488 0.0944432 -4.4056307
[2,] 0.9036215 0.4347075  0.8134117

Figure 6: Aportes independientes, sin forzar su suma a uno.

Figure 7: Aportes independientes, sin forzar su suma a uno.

Interpretación de los resultados

Índice	Aporte representado
DDI	Expansión del aceite y del gas originalmente disuelto
CDI	Reducción de volumen poroso y expansión de agua connata
WDI	Entrada de agua, si el modelo incluye un acuífero
SDI	Expansión de casquete inicial, si el modelo lo incluye

No se usan umbrales universales para clasificar un yacimiento. En este historial y con N ajustado, DDI varía aproximadamente de 0.61 a 0.90 y CDI de 0.09 a 0.43. Roca y agua son especialmente relevantes en los primeros registros; no es correcto asignar un CDI de 7–10% a todo el historial.

La suma DDI+CDI varía aproximadamente entre 0.992 y 1.044. El error de cierre, con el signo definido arriba, va de −4.41% a +0.81%. El ajuste es aproximado, no una identidad exacta impuesta al gráfico.

Un cierre cercano a cero respalda consistencia con los datos y supuestos, pero no valida un modelo único. Si falla, deben revisarse presión promedio, PVT, producción, N y posibles mecanismos omitidos. N es petróleo original en sitio y no equivale a reservas recuperables.

Resumen

Resultado	Valor o criterio
N inicial de comparación	250 MMSTB
N ajustado minimizando error de Np	253.26 MMSTB
N ajustado minimizando error de F (Parte 1)	253.76 MMSTB
DDI / CDI	Calculados independientemente, sin normalización
Verificación	Comparación histórica y error de cierre por registro

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. 11.
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 →