Prueba de Inyectividad en Python: Movilidad, Permeabilidad y Skin

Ejemplo 9.1 de Pressure Transient Testing con la tabla completa de presiones

Analiza una prueba de inyección a caudal constante con los datos tabulados del ejemplo 9.1: diagnóstico log–log, ajuste semilogarítmico, movilidad, skin y comparación con el libro.

¿Qué podemos obtener de una prueba de inyectividad?

Al inyectar a caudal constante, la presión de fondo aumenta. Su evolución permite estimar la movilidad de la formación y la resistencia adicional cerca del pozo, representada por el skin.

Desarrollaremos el ejemplo 9.1 de Pressure Transient Testing, de Lee, Rollins y Spivey. Utilizaremos los 15 pares de tiempo y presión de la tabla 9.1, transcritos del libro, junto con sus propiedades de roca y fluido. No se generan presiones sintéticas.

Esta es una prueba de inyectividad: el tiempo se cuenta desde el comienzo de la inyección y la presión aumenta. En un falloff, el registro comienza después del cierre y la presión disminuye; las ecuaciones de interpretación no deben intercambiarse.

Visita el canal @rigopetrodata

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

Ver canal

Datos del ejercicio

Antes de la prueba, los demás pozos del yacimiento permanecieron cerrados durante varias semanas para estabilizar la presión. El ejemplo considera un yacimiento con inyección de agua y una interpretación de movilidad aproximadamente uniforme.

Parámetro	Valor
Magnitud del caudal de inyección	100 STB/d
Presión inicial, pi	449 psia
Radio del pozo, rw	0.25 ft
Espesor, h	16 ft
Porosidad, φ	0.15
Compresibilidad total, ct	7.7 × 10⁻⁶ psi⁻¹
Viscosidad, μ	1.0 cP
Factor volumétrico, B	1.0 RB/STB
Densidad del fluido en el pozo	62.4 lbm/ft³
Última medición	6.70 h

El libro usa caudal negativo para inyección. Aquí definimos qinj = 100 como magnitud positiva y escribimos las ecuaciones con esa convención. Las presiones son de fondo y absolutas, por lo que no necesitamos una conversión desde cabeza.

Disponemos de la tabla de presión, los parámetros para calcular k y skin y los resultados resueltos para comprobarlos. El libro también proporciona un punto de ajuste de curvas tipo; lo usaremos después como dato gráfico del control, sin presentarlo como un ajuste automático realizado en Python.

1. Leer la tabla de tiempo y presión

Descargar los datos del ejemplo 9.1. Guarda el archivo junto al notebook. Necesitamos pandas, NumPy y Matplotlib: pip install pandas numpy matplotlib.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

datos = pd.read_csv("inyectividad_ejemplo_9_1.csv")
qinj, B, mu = 100, 1.0, 1.0
pi, h, phi, ct, rw = 449, 16, 0.15, 7.7e-6, 0.25

datos["dp_psi"] = datos.p_fondo_psia - pi
datos

	t_h	p_fondo_psia	dp_psi
0	0.0100	546.3	97.3
1	0.0159	587.3	138.3
2	0.0253	642.3	193.3
3	0.0403	711.0	262.0
4	0.0642	785.6	336.6
5	0.1020	855.2	406.2
6	0.1630	910.2	461.2
7	0.2590	949.1	500.1
8	0.4120	976.4	527.4
9	0.6560	997.8	548.8
10	1.0400	1016.5	567.5
11	1.6600	1034.2	585.2
12	2.6400	1051.2	602.2
13	4.2100	1068.0	619.0
14	6.7000	1084.5	635.5

En la tabla 9.2, las diferencias de presión corresponden aparentemente a pi = 448.9 psi, mientras el enunciado indica 449 psi. Conservamos el valor del enunciado; por eso nuestras Δp son 0.1 psi menores. Las presiones transcritas de la tabla 9.1 permanecen intactas.

2. Revisar presión y derivada

La derivada respecto al logaritmo natural del tiempo ayuda a reconocer cambios de régimen:

\[p'=\frac{d\Delta p}{d\ln t}\]

Usaremos una diferencia numérica sencilla con np.gradient. En el gráfico omitimos los extremos de la derivada, donde la estimación es unilateral y menos robusta. Este cálculo no pretende reproducir exactamente el procedimiento de derivación de la tabla 9.2.

datos["derivada_psi"] = np.gradient(datos.dp_psi, np.log(datos.t_h))

plt.loglog(datos.t_h, datos.dp_psi, "o-", label="Incremento de presión")
plt.loglog(datos.t_h.iloc[1:-1], datos.derivada_psi.iloc[1:-1],
           "s-", label="Derivada respecto a ln(t)")
plt.xlabel("Tiempo de inyección (h)")
plt.ylabel("Incremento de presión y derivada (psi)")
plt.legend()
plt.show()

Figura 1: Presión y derivada calculadas con los datos tabulados. El tramo tardío se aproxima a una meseta de derivada.

La respuesta temprana está afectada por almacenamiento y transición. Hacia el final, la derivada se aproxima a un valor constante, compatible con flujo radial. En una recta semilogarítmica de pendiente m, la meseta esperada es m/ln(10), no m.

El libro sitúa el comienzo del tramo semilogarítmico alrededor de 0.53 h. Para reproducir la pendiente de su solución, que usa las mediciones de 1.66 y 6.70 h, ajustaremos los cuatro puntos desde 1.66 h. Más adelante compararemos otras ventanas.

3. Ajustar la recta semilogarítmica

Para una prueba de inyección a caudal constante:

\[p_{wf}=p_{1h}+m\log_{10}\left(\frac{t}{1\,h}\right)\]

La pendiente m es positiva. El intercepto corresponde a la presión sobre la recta a una hora.

radial = datos.loc[datos.t_h >= 1.66]
m, p1h = np.polyfit(np.log10(radial.t_h), radial.p_fondo_psia, 1)

print(f"m = {m:.3f} psi/ciclo")
print(f"p1h = {p1h:.2f} psia")
print(f"Meseta de derivada esperada = {m/np.log(10):.2f} psi")

m = 82.996 psi/ciclo
p1h = 1016.07 psia
Meseta de derivada esperada = 36.04 psi

plt.semilogx(datos.t_h, datos.p_fondo_psia, "o", label="Tabla 9.1")
plt.semilogx(radial.t_h, radial.p_fondo_psia, "s", label="Puntos del ajuste")
t_linea = np.geomspace(1, 6.7, 50)
plt.semilogx(t_linea, p1h+m*np.log10(t_linea), label="Recta ajustada")
plt.xlabel("Tiempo de inyección (h)")
plt.ylabel("Presión de fondo (psia)")
plt.legend()
plt.show()

Figura 2: Ajuste semilogarítmico de los cuatro puntos tardíos. La línea se extiende hasta una hora para leer p1h.

Obtenemos m ≈ 82.996 psi/ciclo y p1h ≈ 1016.07 psia. La solución gráfica del libro reporta 83 psi/ciclo y 1015 psia. La diferencia de aproximadamente 1 psi en el intercepto proviene del ajuste numérico frente a la lectura gráfica.

🛒

Propiedades De Los Fluidos Del Yacimiento

Disponible en Mercado Libre

Ver oferta →

4. Calcular movilidad, permeabilidad y skin

Con qinj positivo, las ecuaciones de campo son:

\[\frac{k}{\mu}=\frac{162.6q_{inj}B}{mh},\qquad k=\frac{162.6q_{inj}B\mu}{mh}\]

\[s=1.151\left[\frac{p_{1h}-p_i}{m} -\log_{10}\left(\frac{k}{\phi\mu c_t r_w^2}\right)+3.23\right]\]

Se utilizan mD, cP, ft, horas, psi y STB/d. Los coeficientes de la ecuación de skin incorporan esas unidades.

movilidad = 162.6*qinj*B/(m*h)
k = movilidad*mu
skin = 1.151*((p1h-pi)/m
             - np.log10(k/(phi*mu*ct*rw**2)) + 3.23)
dp_skin = 141.2*qinj*B*mu*skin/(k*h)

print(f"Movilidad = {movilidad:.3f} mD/cP")
print(f"Permeabilidad efectiva = {k:.3f} mD")
print(f"Skin = {skin:.3f}")
print(f"Presión adicional por skin = {dp_skin:.1f} psi")

Movilidad = 12.245 mD/cP
Permeabilidad efectiva = 12.245 mD
Skin = 2.110
Presión adicional por skin = 152.1 psi

Los resultados son k ≈ 12.245 mD y s ≈ 2.110, consistentes con 12.2 mD y 2.1 del libro. El skin positivo indica resistencia adicional a la inyección. Puede representar daño de formación y otros efectos de terminación; no identifica por sí solo su causa.

La movilidad y la permeabilidad tienen el mismo valor numérico porque μ = 1 cP, pero son magnitudes diferentes. En condiciones multifásicas, la permeabilidad interpretada puede ser efectiva al agua, no necesariamente la absoluta de la roca.

5. Comprobar la selección de la ventana

Comparemos tres inicios posibles, manteniendo la última medición. Esto permite ver cuánto afectan al resultado los puntos que todavía se aproximan al régimen radial.

ventanas = []
for inicio in [0.53, 1.04, 1.66]:
    tramo = datos.loc[datos.t_h >= inicio]
    pendiente, intercepto = np.polyfit(np.log10(tramo.t_h), tramo.p_fondo_psia, 1)
    permeabilidad = 162.6*qinj*B*mu/(pendiente*h)
    s = 1.151*((intercepto-pi)/pendiente
              - np.log10(permeabilidad/(phi*mu*ct*rw**2)) + 3.23)
    ventanas.append([inicio, len(tramo), pendiente, permeabilidad, s])

pd.DataFrame(ventanas, columns=["Inicio_h", "Puntos", "m", "k_mD", "Skin"])

	Inicio_h	Puntos	m	k_mD	Skin
0	0.53	6	85.581479	11.874649	1.867266
1	1.04	5	83.966282	12.103073	2.016392
2	1.66	4	82.995814	12.244593	2.109774

El umbral de 0.53 h incluye desde la primera medición disponible posterior, 0.656 h. La elección del tramo no debe basarse únicamente en conseguir el valor esperado: también debe ser consistente con la derivada, el almacenamiento y la historia de caudal.

6. Revisar el radio investigado y el banco de agua

El libro utiliza la siguiente estimación de radio de investigación:

\[r_i=\sqrt{\frac{kt}{948\phi\mu c_t}}\]

Además, considera dos años de inyección previa a 100 STB/d y un incremento de saturación de agua de 0.4. El radio del banco se estima mediante un balance volumétrico:

\[W_i=q_{inj}B(2\times365),\qquad r_{wb}=\sqrt{\frac{5.615W_i}{\pi h\phi\Delta S_w}}\]

tiempos = np.array([0.53, 6.70])
radios = np.sqrt(k*tiempos/(948*phi*mu*ct))
Wi = qinj*B*2*365
delta_Sw = 0.4
r_banco = np.sqrt(5.615*Wi/(np.pi*h*phi*delta_Sw))

print(f"Radio a 0.53 h: {radios[0]:.1f} ft")
print(f"Radio a 6.70 h: {radios[1]:.1f} ft")
print(f"Radio del banco de agua: {r_banco:.1f} ft")

Radio a 0.53 h: 77.0 ft
Radio a 6.70 h: 273.7 ft
Radio del banco de agua: 368.7 ft

t = np.linspace(0.01, 6.70, 100)
plt.plot(t, np.sqrt(k*t/(948*phi*mu*ct)), label="Radio investigado")
plt.axhline(r_banco, linestyle="--", label="Banco de agua")
plt.xlabel("Tiempo de inyección (h)")
plt.ylabel("Radio (ft)")
plt.legend()
plt.show()

Figura 3: Radio de investigación estimado frente al radio del banco de agua. Ambos son aproximaciones del modelo.

Los radios son aproximadamente 77 ft al inicio radial, 274 ft al final y 369 ft para el banco de agua. El libro reporta 77, 273 y 369 ft; el pequeño cambio del radio final se explica por usar k sin redondear.

El radio investigado permanece dentro del banco de agua estimado. Esto respalda una interpretación de movilidad uniforme en la zona investigada, bajo los supuestos de banco radial y propiedades homogéneas. No demuestra que todo el yacimiento tenga la misma movilidad.

7. Comprobación adicional con el punto de curvas tipo

La solución también obtiene \(C_D\approx355\) y \(s\approx2.0\) mediante curvas tipo. Podemos reproducir sus operaciones usando los valores que el libro lee de la figura 9.5:

• Tiempo de ajuste: 0.8 h.
• Coordenada adimensional \(t_D/C_D\): 100.
• Parámetro de la curva: \(C_De^{2s}=2\times10^4\).

CD = 0.0002637*k/(phi*mu*ct*rw**2) * (0.8/100)
C = CD*phi*ct*h*rw**2/0.8936
skin_tipo = 0.5*np.log(2e4/CD)

print(f"CD = {CD:.1f}")
print(f"C = {C:.6f} bbl/psi")
print(f"Skin por punto de curvas tipo = {skin_tipo:.2f}")

CD = 357.8
C = 0.000463 bbl/psi
Skin por punto de curvas tipo = 2.01

Se obtienen aproximadamente CD = 358 y s = 2.01. La pequeña diferencia con CD = 355 del libro procede de los valores redondeados del ajuste gráfico. Este bloque utiliza el punto de ajuste suministrado por el libro; no realiza una nueva búsqueda ni un ajuste completo de curvas tipo.

8. Contrastar con la solución del libro

control = pd.DataFrame({
    "Magnitud": ["m (psi/ciclo)", "p1h (psia)", "k (mD)", "Skin semilog",
                 "Radio final (ft)", "Banco de agua (ft)", "CD", "Skin curvas tipo"],
    "Python": [m, p1h, k, skin, radios[-1], r_banco, CD, skin_tipo],
    "Libro": [83, 1015, 12.2, 2.1, 273, 369, 355, 2.0]
})
control.round(3)

	Magnitud	Python	Libro
0	m (psi/ciclo)	82.996	83.0
1	p1h (psia)	1016.067	1015.0
2	k (mD)	12.245	12.2
3	Skin semilog	2.110	2.1
4	Radio final (ft)	273.725	273.0
5	Banco de agua (ft)	368.660	369.0
6	CD	357.835	355.0
7	Skin curvas tipo	2.012	2.0

Los resultados numéricos concuerdan con la interpretación publicada dentro del redondeo y la lectura gráfica declarados. La presión inicial se conserva en 449 psia; el intercepto y la pendiente se calculan a partir de las presiones tabuladas, sin modificar los datos para forzar coincidencias.

El ejercicio muestra por qué conviene interpretar primero el régimen de flujo, después estimar movilidad y skin, y finalmente revisar si el radio investigado es compatible con el modelo de fluido adoptado.

Referencia

Lee, J., Rollins, J. B. y Spivey, J. P. (2003). Pressure Transient Testing. Society of Petroleum Engineers. Capítulo 9, Injection-Well Testing, sección 9.2, ejemplo 9.1, tablas 9.1 y 9.2, figuras 9.3–9.5 y ecuaciones 9.3–9.4. Páginas impresas 168–171 (páginas 182–185 del archivo PDF consultado).