Simulación de Pruebas de Presión con el Algoritmo de Stehfest

Inversión numérica de Laplace para verificar k, C y skin con datos reales

Yacimientos
Pruebas de Presión
Python
R

Usa el algoritmo de Stehfest para simular la respuesta de presión de un pozo con almacenamiento y daño, y compara con datos reales de una prueba de decremento. Código en Python y R.

Author

Rigoberto Chandomi

Published

August 1, 2026

Introducción

Cuando analizas una prueba de presión — ya sea por método semilog, curvas tipo o TDS — obtienes valores de permeabilidad (\(k\)), coeficiente de almacenamiento (\(C\)) y factor de daño (\(s\)). Pero, ¿cómo verificas que esos parámetros son correctos?

La respuesta es simular la respuesta de presión con esos parámetros y compararla con los datos reales. Si la simulación ajusta los datos, los parámetros son confiables.

El problema es que la solución analítica del yacimiento con almacenamiento y daño está en el espacio de Laplace — no se puede invertir analíticamente al espacio del tiempo. Ahí es donde entra el algoritmo de Stehfest (1970): una técnica numérica que invierte la transformada de Laplace usando solo 8 coeficientes precalculados.

¿Por qué Stehfest?

Es el método estándar en la industria para inversión numérica de Laplace en pruebas de presión. Software como Saphir, PanSystem y Ecrin lo usan internamente. Con este algoritmo, puedes hacer lo mismo en 20 líneas de código.

Video

El Algoritmo de Stehfest

La inversión numérica se realiza con la fórmula:

\[f(t) = \frac{\ln(2)}{t} \sum_{i=1}^{N} V(i) \cdot \bar{f}\left(i \cdot \frac{\ln(2)}{t}\right)\]

Donde \(\bar{f}(s)\) es la función en el espacio de Laplace evaluada en \(s = i \cdot \ln(2)/t\), y \(V(i)\) son los coeficientes de Stehfest precalculados.

Para \(N = 8\):

Los 8 coeficientes de Stehfest con N=8. Estos números son todo lo que necesitas.
\(V(1)\) \(V(2)\) \(V(3)\) \(V(4)\) \(V(5)\) \(V(6)\) \(V(7)\) \(V(8)\)
-0.3333 48.3333 -906 5464.6667 -14376.6667 18730 -11946.6667 2986.6667

Solución en Laplace: Yacimiento con WBS y Skin

La presión adimensional en el espacio de Laplace para un yacimiento homogéneo infinito con almacenamiento y daño es:

\[\bar{P}_{wD} = \frac{1}{u} \left[ \frac{K_0(\sqrt{u}) + s\sqrt{u}K_1(\sqrt{u})}{\sqrt{u}K_1(\sqrt{u}) + C_D u \left[K_0(\sqrt{u}) + s\sqrt{u}K_1(\sqrt{u})\right]} \right]\]

Donde \(K_0\) y \(K_1\) son las funciones de Bessel modificadas de segunda clase, y las variables adimensionales son:

\[t_D = \frac{0.0002637 \cdot k \cdot t}{\phi \mu c_t r_w^2} \qquad C_D = \frac{0.8936 \cdot C}{\phi c_t h r_w^2} \qquad p_D = \frac{k h \cdot \Delta p}{141.2 q B \mu}\]

Datos del Ejemplo

Datos del pozo y yacimiento.
Parámetro Valor Unidades
Qo 125 STB/D
h 32 ft
φ 0.22 fracción
Bo 1.125 RB/STB
Pi 2,750 psia
ct 1.09×10⁻⁵ psi⁻¹
rw 0.25 ft
μ 2.122 cp

Parámetros estimados (del análisis previo con curvas tipo):

Parámetro Valor
k 22.92 md
s 5.78
C 0.005 bbl/psi

Paso 1 — Cargar Datos y Gráfico Diagnóstico

Ver código 
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import kn  # Funciones de Bessel K0 y K1
import pandas as pd

# === DATOS ===
Qo = 125; h = 32; phi = 0.22; Bo = 1.125
Pi = 2750; ct = 1.09e-5; rw = 0.25; vis = 2.122

# Resultados del análisis previo
k = 22.92; s = 5.78; C = 0.005

# Cargar datos de la prueba
data = pd.read_csv('Datos_DD.csv')
t = data['t'].values
pwf = data['pwf'].values

# Calcular Δp y derivada
dp = Pi - pwf
dpdt = np.zeros(len(t))
for i in range(1, len(t)):
    dpdt[i] = (dp[i] - dp[i-1]) / (np.log(t[i]) - np.log(t[i-1]))

# Gráfico diagnóstico
fig, ax = plt.subplots(figsize=(10, 7))
ax.loglog(t[1:], dp[1:], 'o', ms=6, color='#2980b9', alpha=0.7, label='Δp (datos)')
ax.loglog(t[1:], dpdt[1:], 'o', ms=6, color='#2d8a4e', alpha=0.7, label="tΔp' (datos)")
ax.set_xlabel('Tiempo (hr)', fontsize=12)
ax.set_ylabel('Δp, tΔp\' (psi)', fontsize=12)
ax.set_title('Gráfico Diagnóstico Log-Log', fontsize=14, fontweight='bold')
ax.legend(fontsize=11); ax.grid(True, which='both', alpha=0.3)
plt.tight_layout(); plt.show()
Figure 1: Gráfico diagnóstico log-log: Δp y derivada de Bourdet vs tiempo.
Ver código 
library(ggplot2)

Qo <- 125; h <- 32; phi <- 0.22; Bo <- 1.125
Pi <- 2750; ct <- 1.09e-5; rw <- 0.25; vis <- 2.122
k <- 22.92; s <- 5.78; C <- 0.005

data_DD <- read.csv('Datos_DD.csv')
data_DD$dp <- Pi - data_DD$pwf
data_DD$dpdt <- c(0, diff(data_DD$dp) / diff(log(data_DD$t)))

ggplot(data_DD[-1,]) +
  geom_point(aes(x=t, y=dp, color="Δp"), size=2, alpha=0.7) +
  geom_point(aes(x=t, y=dpdt, color="tΔp'"), size=2, alpha=0.7) +
  scale_x_log10() + scale_y_log10() +
  scale_color_manual(values=c("Δp"="#2980b9", "tΔp'"="#2d8a4e")) +
  labs(x="Tiempo (hr)", y="Δp, tΔp' (psi)",
       title="Gráfico Diagnóstico Log-Log", color="") +
  theme_minimal() + theme(plot.title=element_text(face="bold", size=14))
Figure 2: Gráfico diagnóstico log-log.

Paso 2 — Implementar el Algoritmo de Stehfest

Ver código 
def stehfest_inversion(tD, cD, skin):
    """
    Inversión numérica de Laplace con el algoritmo de Stehfest (N=8).
    Calcula PwD en el espacio del tiempo para un yacimiento homogéneo
    infinito con almacenamiento (cD) y daño (skin).
    """
    V = np.array([-0.3333, 48.3333, -906, 5464.6667,
                  -14376.6667, 18730, -11946.6667, 2986.6667])
    
    m = len(tD)
    PwD = np.zeros(m)
    
    for j in range(m):
        a = np.log(2) / tD[j]
        i_vec = np.arange(1, len(V) + 1)
        u = i_vec * a
        ru = np.sqrt(u)
        
        # Funciones de Bessel modificadas K0 y K1
        K0 = kn(0, ru)
        K1 = kn(1, ru)
        
        # Solución en Laplace: yacimiento infinito con WBS y skin
        numerator = K0 + skin * ru * K1
        denominator = ru * K1 + cD * u * (K0 + skin * ru * K1)
        PwD_laplace = (1/u) * (numerator / denominator)
        
        # Inversión de Stehfest
        PwD[j] = a * np.sum(V * PwD_laplace)
    
    return PwD

print("✅ Función stehfest_inversion() definida")
✅ Función stehfest_inversion() definida
Ver código 
Stehfest_inversion <- function(tD, cD, skin) {
  
  V <- c(-0.3333, 48.3333, -906, 5464.6667,
         -14376.6667, 18730, -11946.6667, 2986.6667)
  
  m <- length(tD)
  PwD <- numeric(m)
  
  for (j in 1:m) {
    a <- log(2) / tD[j]
    i <- 1:length(V)
    u <- i * a
    ru <- sqrt(u)
    
    K0 <- besselK(ru, 0)
    K1 <- besselK(ru, 1)
    
    numerator <- K0 + skin * ru * K1
    denominator <- ru * K1 + cD * u * (K0 + skin * ru * K1)
    PwD_laplace <- (1/u) * (numerator / denominator)
    
    PwD[j] <- a * sum(V * PwD_laplace)
  }
  
  return(PwD)
}

cat("✅ Función Stehfest_inversion() definida\n")
✅ Función Stehfest_inversion() definida

Paso 3 — Simular y Comparar con Datos Reales

Ver código 
# Variables adimensionales
tD = (0.0002637 * k * t) / (phi * vis * ct * rw**2)
cD = (0.8936 * C) / (phi * ct * h * rw**2)

# Simular presión con Stehfest
PwD = stehfest_inversion(tD[1:], cD, s)  # Excluir t=0

# Convertir a dimensionales
dp_sim = (PwD * 141.2 * Bo * vis * Qo) / (h * k)
pwf_sim = Pi - dp_sim

# Derivada de la simulación
dpdt_sim = np.zeros(len(dp_sim))
for i in range(1, len(dp_sim)):
    dpdt_sim[i] = (dp_sim[i] - dp_sim[i-1]) / (np.log(t[i+1]) - np.log(t[i]))

# === GRÁFICO LOG-LOG: DATOS vs SIMULACIÓN ===
fig, ax = plt.subplots(figsize=(11, 7))

# Datos reales (puntos)
ax.loglog(t[1:], dp[1:], 'o', ms=7, color='#2980b9', alpha=0.5, label='Δp (datos)')
ax.loglog(t[1:], dpdt[1:], 'o', ms=7, color='#2d8a4e', alpha=0.5, label="tΔp' (datos)")

# Simulación (líneas)
ax.loglog(t[1:], dp_sim, '-', lw=2.5, color='#2980b9', label='Δp (Stehfest)')
ax.loglog(t[2:], dpdt_sim[1:], '-', lw=2.5, color='#2d8a4e', label="tΔp' (Stehfest)")

ax.set_xlabel('Tiempo (hr)', fontsize=12)
ax.set_ylabel('Δp, tΔp\' (psi)', fontsize=12)
ax.set_title('Datos Reales vs Simulación con Stehfest', fontsize=14, fontweight='bold')
ax.legend(fontsize=10); ax.grid(True, which='both', alpha=0.3)

# Anotar parámetros
ax.text(0.02, 0.05, f'k = {k} md\nC = {C} bbl/psi\ns = {s}',
        transform=ax.transAxes, fontsize=11, verticalalignment='bottom',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

plt.tight_layout(); plt.show()
Figure 3: Comparación datos reales vs simulación con Stehfest. Las líneas continuas (simulación) ajustan los puntos (datos reales), validando los parámetros k, C y s.
Ver código 
print(f"\n✅ Simulación completada")

✅ Simulación completada
Ver código 
print(f"   Δp simulado final: {dp_sim[-1]:.1f} psi")
   Δp simulado final: 760.6 psi
Ver código 
print(f"   Δp medido final:   {dp[-1]:.1f} psi")
   Δp medido final:   808.9 psi
Ver código 
print(f"   Error: {abs(dp_sim[-1]-dp[-1])/dp[-1]*100:.1f}%")
   Error: 6.0%
Ver código 
tD <- (0.0002637 * k * data_DD$t) / (phi * vis * ct * rw^2)
cD <- (0.8936 * C) / (phi * ct * h * rw^2)

PwD <- Stehfest_inversion(tD[-1], cD, s)

data_DD$dp_sim <- c(0, (PwD * 141.2 * Bo * vis * Qo) / (h * k))
data_DD$dpdt_sim <- c(0, diff(data_DD$dp_sim) / diff(log(data_DD$t)))

ggplot(data_DD[-1,]) +
  geom_point(aes(x=t, y=dp, color="Δp datos"), size=2, alpha=0.5) +
  geom_point(aes(x=t, y=dpdt, color="tΔp' datos"), size=2, alpha=0.5) +
  geom_line(aes(x=t, y=dp_sim, color="Δp Stehfest"), linewidth=1) +
  geom_line(aes(x=t, y=dpdt_sim, color="tΔp' Stehfest"), linewidth=1) +
  scale_x_log10() + scale_y_log10() +
  scale_color_manual(values=c("Δp datos"="#2980b9", "tΔp' datos"="#2d8a4e",
                               "Δp Stehfest"="#1a5276", "tΔp' Stehfest"="#1a472a")) +
  labs(x="Tiempo (hr)", y="Δp, tΔp' (psi)",
       title="Datos Reales vs Simulación con Stehfest", color="") +
  theme_minimal() + theme(plot.title=element_text(face="bold", size=14))
Figure 4: Datos reales vs simulación con Stehfest.

Paso 4 — Gráfico de Presión de Fondo

Ver código 
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(t, pwf, 'o', ms=5, color='#2980b9', alpha=0.5, label='Pwf (datos)')
ax.plot(t[1:], pwf_sim, '-', lw=2.5, color='#c0392b', label='Pwf (Stehfest)')
ax.set_xlabel('Tiempo (hr)', fontsize=12)
ax.set_ylabel('Pwf (psia)', fontsize=12)
ax.set_title('Presión de Fondo: Datos vs Simulación', fontsize=14, fontweight='bold')
ax.legend(fontsize=11); ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()
Figure 5: Presión de fondo fluyente: datos medidos vs simulación.
Ver código 
data_DD$pwf_sim <- Pi - data_DD$dp_sim

ggplot(data_DD) +
  geom_point(aes(x=t, y=pwf, color="Datos"), size=2, alpha=0.5) +
  geom_line(aes(x=t, y=pwf_sim, color="Stehfest"), linewidth=1) +
  scale_color_manual(values=c("Datos"="#2980b9", "Stehfest"="#c0392b")) +
  labs(x="Tiempo (hr)", y="Pwf (psia)",
       title="Presión de Fondo: Datos vs Simulación", color="") +
  theme_minimal() + theme(plot.title=element_text(face="bold", size=14))
Figure 6: Pwf medido vs simulado.

¿Qué Aprendimos?

El algoritmo de Stehfest permite cerrar el ciclo del análisis de pruebas de presión:

Datos → Análisis (semilog, curvas tipo, TDS) → k, C, s → Stehfest → Simulación → Verificación

Si la simulación ajusta los datos, los parámetros son confiables. Si no ajusta, hay que revisar el modelo o los parámetros. Este paso de verificación es el que separa un análisis profesional de uno académico.

Mejora potencial

El siguiente paso sería usar regresión no lineal para optimizar k, C y s automáticamente — minimizar la diferencia entre los datos y la simulación de Stehfest. Eso lo veremos en un post futuro.

Referencias

  • Stehfest, H. (1970). Algorithm 368: Numerical Inversion of Laplace Transforms. Communications of the ACM, 13(1).
  • Spivey, J. & Lee, J. (2013). Applied Well Test Interpretation. SPE.
  • Lee, J., Rollins, J. & Spivey, J. (2003). Pressure Transient Testing. SPE.
  • Sun, H. (2015). Advanced Production Decline Analysis. Gulf Professional Publishing.
🛒
Propiedades De Los Fluidos Del Yacimiento

Disponible en Mercado Libre

Ver oferta →