Estudio regional para la evaluación de la formación Vaca Muerta: generación de transecta geofísica
Integración con estudio geomecánico e interpretación secuencial semiautomática. Exposición mural de gigantografía Neuquén, cuenca Neuquina
Por Juan Tavella (JT Reservorios), Javier Carrero (consultor), David Epelboim (Geoinfo), Sigfrido Nielsen (Geoinfo),
Federico Späth (YPF - coordinador), Hugo Aguirre (Shell Argentina), José Luis Capuano (Capex), Gabriel Chao (TotalEnergies)
Tomás D’Biassi (YPF), Victoria Lazzari (Phoenix Global Resources), Diego Lenge (Oilstone),
Enzo Luna Carabajal (Pampa Energía), Juan Quiroga (Vista Energy), Neptalí Requena (Tecpetrol),
María Dolores Vallejo (Chevron Argentina), Luis Vernengo (PanAmerican Energy) y Emilse Zunino (Pluspetrol)
Este estudio complementa el presentado en Conexplo 2014 y se apoya en la información aportada por 12 compañías operadoras, con más de 5000 km2 de información sísmica migrada preapilada en tiempo provenientes de 19 relevamientos 3D y 2D y la información de 17 pozos. La sísmica es seleccionada a lo largo de un callejón de 300 m de ancho y más de 360 km de longitud y tratada adecuadamente para convivir en una grilla 3D de subsuperficie unificada de más de 22.000 km2.
El yacimiento no convencional asociado a la formación Vaca Muerta del Jurásico Superior de la cuenca Neuquina abarca una extensión de 36.000 km2 en las provincias de Neuquén, Mendoza, Río Negro y La Pampa. Este estudio se basó en un perfil regional de orientación preferencial norte-sur, que se extiende desde el flanco norte del alto estructural conocido como Dorsal de Huincul, dirigiéndose hacia el norte para atravesar buena parte de la zona de máximo potencial de la formación. Luego de atravesar otro alto: el Dorso de los Chihuidos, remata cerca del límite entre las provincias de Neuquén y Mendoza, donde la formación pierde desarrollo e interés económico (Figura 1).
En el desarrollo de los no convencionales a nivel mundial ha sido de mucha importancia la definición de nuevos flujos de trabajo para la caracterización y la definición de la estrategia de desarrollo de los campos. Factores como la geoquímica y geomecánica han tomado mayor protagonismo en la estimación del comportamiento de producción esperado de los pozos. Ya no solo es importante encontrar zonas ricas en hidrocarburo, sino también identificar los niveles con mayor eficiencia en el proceso de fracturación hidráulica, necesaria dada la baja permeabilidad primaria de la formación. Es decir, el reservorio como tal no producirá sin una correcta operación de estimulación por más que se tengan las condiciones geológicas y petrofísicas de mayor prospectividad en la zona. Aunque este trabajo no tiene por objetivo directo la evaluación de un área específica para definición de locaciones, igualmente se beneficia de la experiencia adquirida en los últimos 15 años en el desarrollo de los flujos de trabajo mencionados.
El estudio, que se complementaría con el realizado en 2014 y presentado en CONEXPLO 2014 (González et al., 2016), se apoya en el aporte de la información de 12 compañías operadoras, actores principales en la historia reciente del desarrollo de la formación: se recibieron más de 5000 km2 de información sísmica migrada preapilado en tiempo provenientes de 19 relevamientos 3D y 2D y la información de 17 pozos. La sísmica será seleccionada a lo largo de un callejón de 300 m de ancho y más de 360 km de longitud y tratada adecuadamente para convivir en una grilla 3D de subsuperficie unificada de más de 22.000 km2.
La estrategia de procesamiento adoptada es audaz: hacer un merge preapilado, pero postmigración. Es decir, se intentará contar con buena parte de los beneficios de un merge desde disparos: tener un conjunto consistente capaz de oficiar de entrada a procesos de caracterización preapilado. Pero se evitará un trabajo que estaría totalmente fuera de escala con relación a los objetivos regionales planteados. La imposibilidad de hacer deconvolución y compensación de amplitudes preapilado consistentes en superficie obliga a poner mucho énfasis en el tratamiento de las sumas parciales de ángulo: conformación espectral, ecualización de amplitudes y un innovador tratamiento de amplitud preapilado “AVO compliant” con máximo control de pozo, se constituyeron en actores principales.
Ya en las etapas de caracterización e integración, la información de los pozos, previamente editada y normalizada, se usará para generar un juego de electrofacies, clave para dar sentido a las sismofacies y estimación del carbono orgánico total .Se aprovecharán los resultados del trabajo de procesamiento para hacer una inversión preapilado, que se complementará con las bondades de las redes neurales, todo ayudado por el abundante y bien distribuido set de puntos de control aportado por los pozos. Se estará en condiciones de predecir la impedancia S para aumentar el alcance de la caracterización. Así es como se podrán estimar parámetros clave para entender el comportamiento geomecánico de la formación e integrar los resultados al estudio de la transecta geomecánica elaborada en paralelo con este trabajo (Sosa Massaro et al., 2022).
Se contará con el auxilio de herramientas de interpretación estratigráfica semiautomática (Qayyum et al. 2012, 2013), que permitirá definir superficies de alcance regional, techo y base del sistema Quintuco-Vaca Muerta y siete límites de secuencias. Estás últimas son fundamentales para delimitar unidades y clasificar en sismofacies.
En el terreno de la interpretación cuantitativa se confiará en las sismofacies, que controladas por las electrofacies, se espera que permitan identificar, propagar y mapear las condiciones del paleoambiente que, como es aceptado, controlan buena parte de las características que hoy exhibe la formación. Los resultados serán expuestos en la exhibición mural de gigantografía y posters en CONEXPLO 2022.
En la figura 2 se observa el flujo de trabajo diseñado para cumplir con los objetivos mencionados.
Predicción del contenido orgánico total (COT)
Desde los estudios de las rocas madres en la exploración de las diversas cuencas en el mundo, se han desarrollado metodologías empíricas para la predicción del contenido de materia orgánica a partir de los registros de pozo. Este es el caso de Schmoker (1981) y Fertl y Rieke (1980) que emplean el registro de GR (rayos gamma) y GR espectral para establecer relaciones entre los valores de radiactividad y la materia orgánica. Dallenbach et al. (1983) combinan el uso del registro de GR con el tiempo de tránsito compresional para obtener un parámetro derivado que se relaciona linealmente con la riqueza en materia orgánica.
Por su parte, Schmoker (1979), Schmoker y Hester (1983) y Abdel-Rahman y Wali (1984) utilizan el registro de densidad para determinar la presencia de arcillas bituminosas.
Para la determinación de presencia de carbón y horizontes ricos en materia orgánica Rieke et al. (1980), Lawrence et al. (1984) y Herron (1986) emplearon herramientas del tipo carbono/oxígeno.
Nixon (1973), Meissner (1978), Mendelson (1985), Schmoker y Hester (1989) y Morel (1999) estudiaron el impacto que genera la presencia de materia orgánica sobre las lecturas de los perfiles resistivos, así como las distintas concentraciones producen variaciones en los registros.
Veiga y Orchuela (1989) a partir del uso combinado de los registros de GR, sónico, radiactivos y eléctricos identificaron niveles generadores de hidrocarburos en la formación Vaca Muerta.
Carpentier et al. (1989) combinan el uso de las propiedades resistivas y de propagación de ondas en la estimación de la materia orgánica.
Passey et al. (1990) desarrollaron una metodología denominada “D log R” sobre la base del análisis de la respuesta que los registros de tiempo de tránsito compresional y resistividad exhiben frente a la presencia de materia orgánica.
Huang y Williamson (1996) utilizaron técnicas basadas en la aplicación de redes neuronales para la cuantificación de la materia orgánica a partir del uso combinado de registros de GR, tiempo de tránsito compresional y resistivos.
Cuervo et al. (2016) consigue una regresión múltiple entre los registros de densidad, sónico y resistivo para la predicción de la totalidad del carbono orgánico (incluyendo hidrocarburo in situ) en la formación Vaca Muerta.
Para el análisis regional de la transecta se utilizaron los datos de laboratorio de materia orgánica de los pozos MBE.x-2, LCa.xp-2001, LCav.x-8 y ET.xp-2002, los cuales se encuentran distribuidos a lo largo de la sección en estudio, y se entrenó una red neural a partir de los registros Gamma ray, densidad, sónico compresional y resistividad. De esta manera para cada uno de los pozos se pudo obtener una curva de COT (Figura 3).
Para la predicción a lo largo de toda la transecta del contenido de materia orgánica, se buscó una relación cuadrática empírica general de tal forma que se pudiese realizar el cálculo a partir del atributo generado de la inversión sísmica. En la figura 4, en su parte derecha, se muestra el gráfico cruzado entre la Impedancia P y la predicción de COT de todos los pozos a lo largo de la transecta donde se obtuvo a relación empírica para luego extenderlo al dominio sísmico. En la parte izquierda, en la pista para cada pozo se ve, de izquierda a derecha: la impedancia y la comparación entre el COT de laboratorio (puntos rojos), la predicción del COT por redes neurales (curva azul) y el COT calculado a partir de la impedancia P según la relación empírica conseguida.
Clasificación de electrofacies
Para la clasificación de electrofacies se realizó un análisis de clúster (Figura 5) para el entrenamiento de una red neural a partir de los registros de densidad, la predicción de TOC, sónico compresional, gamma ray y resistividad. De acuerdo con las descripciones litológicas de la Formación Vaca Muerta y de los rangos de valores de los registros de la clasificación, se le asignó los nombres a cada facies. En la figura 6 se observa un corte con todos los pozos que muestra la correlación de las electrofacies con el COT.
La transecta se construyó a partir del aporte de 12 compañías operadoras que cedieron el uso de información sísmica preapilado y posapilado de 19 procedencias diferentes, según el detalle de la tabla 1. Se completó con la información de los pozos, cuya ubicación se puede apreciar en el mapa de la figura 1, junto con los relevamientos sísmicos y la traza adoptada para la transecta.
Carga y selección de datos
La traza adoptada presenta tendencia general norte-sur, aunque con varios quiebres para adaptarla al dato sísmico y de pozos con que se contaba. La traza se aprecia en el mapa de la figura 1 y su definición analítica por medio de sus vértices en la tabla 2.
Se cargó en forma individual e independiente cada conjunto de dato con sus propios parámetros geométricos. Luego, a partir de lo cargado, se procedió a selecciona para todos los datos de una faja de 300 m de ancho alrededor de la traza adoptada, con el objetivo de permitir correr con eficiencia procesos multitraza.
Adecuación de la información
A cada conjunto se le aplicó una compensación estática de tiempo constante (gual valor para todas las trazas) calculado a partir de los valores nominales de plano de referencia y velocidad de reemplazo, con el objeto de referenciarlos a una superficie de 700 m de altitud respecto del nivel del mar. Los valores en tiempo aplicados a cada conjunto, así como la velocidad de reemplazo y plano original usados para el cálculo, se detallan en la tabla 1.
Se aplicaron correcciones de fase, salto en tiempo y amplitud, que se calcularon con base en los apilados de cada conjunto. Posteriormente, se separará el dato en tres sumas parciales de ángulo de incidencia y se refinarán las compensaciones a las sumas de ángulos medios, como se explica en el apartado “Procesamiento unificado”. A continuación, se presenta una breve explicación de cada compensación.
Fase
El objetivo es llevar toda la sísmica a una misma relación de fase, condición clave para los datos que tienen como destino fusionarse para, finalmente, constituirse en la entrada a la inversión de traza. Se llevó a cabo en dos etapas: la primera, basada en el modelo convolucional de la traza sísmica del apilado y, otra posterior, con los datos ya en la grilla unificada de carácter empírico, de ajuste fino, entre aquellos datos que mostraron todavía diferencias con sus vecinos. Esta última corrección se hizo sobre la base de las sumas de ángulo medio.
En este primer abordaje de la compensación por fase se confió en que el modelo convolucional, a través del amarre pozo-sísmica, proporcionara el valor medio de rotación fase a aplicar al dato para optimizar la correlación con el pozo. Si bien la estimación de la fase a partir de lo observado en un pozo tiene su debilidad, se consideró el procedimiento más sólido para abordar el problema. La segunda etapa optimizará la fase de algunos datos, a la vez que oficiará de control y revisión final. En la tabla 3 se pueden apreciar los valores de compensación de fase aplicados.
Tiempo
Al igual que para la fase, es deseable que los datos que deberán fusionarse en el procesamiento unificado exhiban concordancia en la posición temporal de los marcadores sísmicos en la zona de overlap o de aproximación. En el paso anterior se aplicaron compensaciones temporales en base a parámetros nominales para llevar a una referenciación común (700 m). Sin embargo, dado que cada conjunto fue procesado con distintos procedimientos para el cálculo de correcciones estáticas, es normal que persistan diferencias. Aunque de segundo orden, estas diferencias deben ser compensadas para garantizar los resultados en zonas de unión. Estas correcciones estáticas (un valor por conjunto) se determinaron por correlación cruzada de las trazas de los apilados y elección de valores de compromiso. En una etapa posterior, ya en la grilla unificada, se aplicarán correcciones estáticas de borde y basculaciones para atacar las diferencias remanentes. Los valores de estática aplicadas a cada dato se observan en la tabla 3.
Amplitud
La amplitud global de un dato sísmico es totalmente relativa dependiendo de los procedimientos de escalado que la modifican arbitrariamente. Esto explica los niveles tan dispares de amplitud de los datos recibidos, que debieron ser compensados en una primera etapa aplicando un escalar a cada conjunto. Para estimarlos se calcularon factores inversamente proporcionales a la amplitud media de los apilados en la ventana Quintuco-Tordillo, representativa de la zona de interés. Para este fin se interpretaron en forma expeditiva versiones preliminares de los marcadores mencionados. La aplicación de estos factores multiplicativos llevó las amplitudes a valores compatibles en todos los datos. Posteriormente, una vez aplicada la conformación espectral, como se verá en “Procesamiento unificado”, se calcularán ya con base en las sumas de ángulo medio, escalares más enfocados. Adicionalmente, se aplicarán correcciones de amplitud de bordes en los datos que lo requirieron.
Al igual que la fase, el espectro de amplitud debe guardar consistencia en los datos que se integrarán. Para compensar las diferencias naturales en contenido espectral, la estrategia adoptada consistió, primeramente, en una restricción a la banda pasante 12-45 hz por medio de la aplicación de un filtro de esquinas 8-12-45-55 hz. Posteriormente, para compensar bandas de frecuencia decaídas no se usó balanceo espectral formal, sino que se hicieron refuerzos puntuales para lograr un espectro razonablemente plano dentro de la banda pasante.
En la parte izquierda de la figura 7 están los espectros de los 19 datos recibidos después de aplicado el filtro pasabanda. En color se destacan solo aquellos siete que exhibieron apartamientos importantes de la condición plana dentro de la banda pasante. En la parte derecha se muestra como luce el gráfico luego de las compensaciones aplicadas a estos siete datos: el diseño de filtros pasabanda en las zonas de debilidad y posterior suma produce resultados eficientes. La suma modifica el nivel de amplitud global, que ya puede considerarse definitivo, así los datos están en condiciones de entrar en la etapa de ajuste fino de amplitudes, como se explica a continuación.
Procesamiento unificado
La estrategia adoptada para encarar el procesamiento subsiguiente se basó en considerar al conjunto de datos como un único relevamiento de sísmica tridimensional. Para ello fue necesario definir una grilla unificada de subsuperficie bidimensional que los contuviera a todos. La orientación elegida fue este-oeste para las inlines, justificado en que los datos recibidos exhibían orientaciones, aunque variadas, con predominio de aquella orientación. En cuanto al tamaño de celda, se optó por 60 x 60 m, mayor a cualquiera de los tamaños de celda de origen (Tabla1), con lo que se evitaron las celdas vacías y la necesidad de interpolación. La degradación de resolución espacial asociada a la adopción de una celda amplia se consideró irrelevante en virtud del carácter regional del estudio. El despliegue de esta grilla puede apreciarse en los ejes en color rojo en el mapa de la figura 1. En resumen, los parámetros de esta grilla, de más de 22.000 km2, son los siguientes:
Orientación de inlines (xlines crecientes) Este
Tamaño de celda 60 x 60 m
Coordenadas del origen
(IL 1, XL 1 extremo sudoeste) X = 2444490 Y = 5698290
Rango de Inlines 1 – 3200
Rango de Xlines 1 – 1930
Si se pretende fundir en un único conjunto los diversos datos, además de aplicar los procesos de adecuación descriptos en la sección anterior, es indispensable redefinir los aspectos geométricos vinculados a la grilla. Es decir que hay que proveer a cada traza la relación entre las coordenadas geográficas y el par Inline/Xline que le corresponde en la grilla.
Con los datos en la condición descripta en el apartado anterior se procedió a cargar los 19 conjuntos. En la figura 8 se muestra un perfil de la transecta resultante para un apilado de offset medio. Son casi 400 km
lineales de perfil en el dominio del tiempo sísmico doble. Se puede apreciar que el ensamble es satisfactorio tanto en tiempo (no hay saltos apreciables) como en amplitudes globales, y no muestran fronteras con cambios abruptos. Otro aspecto para destacar es la presencia de zonas con carencia de información entre distintos datos. Estos intervalos serán completados en etapas posteriores a la inversión, con el fin ofrecer una sección continua para hacer más eficientes procesos tipo propagación e interpretación semiautomática y también para mejorar la calidad de despliegue y producir una imagen más elegante en la exhibición de la gigantografía.
A partir de este dato unificado se implementará un flujo de procesamiento tendiente a optimizar el ensamble, mejorar las relaciones de amplitud, tanto posapilado como preapilado, con vistas a la inversión.
El ángulo con que una propagación elástica incide sobre una interfaz es la verdadera variable de la que depende la reflectividad. Por esa razón fue necesario convertir los datos preapilado del dominio del offset al del ángulo de incidencia. En este punto conviene repasar las columnas asociadas a la “configuración recibida” en la tabla 1. Se advierten dos cosas: que algunos datos ya vinieron en el dominio del ángulo, por lo que no necesitarán ser convertidos y que hay datos que llegaron en forma de tres sumas parciales, también en el dominio del ángulo. En la siguiente sección se analizará la estrategia para encarar la inversión haciendo frente a la diversidad de formatos presentes a la entrada.
La conversión a ángulo requiere de un campo velocidad de ondas compresionales sobre el cual trazar rayos y evaluar el ángulo de incidencia para cada nivel. No fue posible contar con velocidades interválicas (PSDM) confiables para todos los datos recibidos, por lo que se decidió generar un modelo de velocidad muy suavizado construido mediante la propagación de los sónicos compresionales con la asistencia de cuatro horizontes.
El hecho de que algunos datos llegaron en forma de tres sumas parciales en ángulo forzó a que el trabajo de inversión se realice precisamente usando tres sumas parciales como entrada. Por cuestiones de consistencia y para evitar discontinuidades no conviene tener multiplicidades distintas a la entrada. Por esta razón es que se decidió degradar aquellos datos que vinieron en forma de gather PSTM, de multiplicidades mucho mayores a 3.
Conviene en este punto recordar que el dato recibido para la línea 2D HV89-202 es un apilado, por no estar disponible ningún tipo de información preapilado. Esta circunstancia impide calcular sumas parciales y se supondrá que el apilado representa la suma parcial de ángulo medio. En la sección “Compensación AVO Compliant” se retoma este tema.
En función de los valores de ángulos de los datos recibidos en forma de sumas parciales (Tabla1), se adoptaron los siguientes valores para degradar los datos recibidos en formato gather (ver tabla a continuación).
De esta forma se considera que el error que se comete al reasignar los valores es insignificante para los objetivos del proyecto, al evitar procesos de reconversión.
Con estos valores se realizaron sumas parciales a partir de los datos recibidos en forma de gathers en offset, que habían sido previamente convertidos al dominio del ángulo de incidencia.
A partir de este punto, habiéndose generado las sumas parciales, se proseguirá la secuencia de procesamiento sobre la suma de ángulos medios. Una vez definida esta secuencia se la hará extensiva a las sumas cercanas y lejanas.
La elección de los ángulos medios como patrón no es caprichosa, sino que obedece al hecho de que es el apilado más robusto, en donde la mayoría de los relevamientos tienen mejor población de trazas y la consiguiente superioridad en calidad y relación señal/ruido.
Las estimaciones de las compensaciones de fase se hicieron por correlación cruzada sobre los datos confluyentes en cada unión, aplicando distintos valores de fase a uno de ellos. Se aplicó la rotación de fase que produjo la correlación óptima. Con esta metodología se revisó cada unión, pero ahora con base en la suma media en vez del apilado; solo tres datos merecieron compensación adicional, como se observa en tabla 3.
En relación con la amplitud, recordemos que, a diferencia del primer paso de compensación, ahora se evalúan sobre los apilados medios que, además, ya tienen la conformación espectral aplicada. Los escalares resultantes se observan en la tabla 3.
Adicionalmente, se aplicaron compensaciones de amplitud asociadas a deficiencias en bordes, en general decaimientos, hecho muy frecuente en virtud de la baja multiplicidad propia de los bordes y su efecto negativo sobre la migración. Estas correcciones tomaron la forma de un mapa de escalares multiplicativos. Como se observa en la tabla 3, solo 7 de los 19 datos necesitaron de estas correcciones. Estas compensaciones, junto con las de estáticas de borde, que se explicarán en la siguiente sección, serán las responsables de la eficiencia de las “transiciones suaves” que se practicarán.
Las zonas de overlap o superposición entre dos datos fueron optimizadas por transición suave, lo que requiere cierto grado de similitud entre los datos duplicados. Las compensaciones de amplitud de borde explicadas contribuyen a esta similitud. Pero no solo las amplitudes son deficientes en las zonas de borde, también la forma estructural de los marcadores sísmicos se ve comprometida. Esta deformación, que atenta contra la similitud, puede ser compensada por una corrección estática. Para calcularla se hizo una doble interpretación en la zona de superposición —una sobre cada dato— y, en base a la diferencia se diseñó la corrección. En la figura 9 se ilustra este procedimiento tomando como ejemplo la zona de superposición entre Aguada Baguales y Loma La Lata. A la izquierda se ven secciones de ambos datos alrededor de la zona de superposición; las líneas roja y azul son las interpretaciones de un marcador importante para Loma La Lata y Aguada Baguales, respectivamente. Se representan además en línea punteada, por comparación, para apreciar la distorsión estructural relativa entre ambos conjuntos en los bordes, atribuida a deficiencias en la migración. Como se mencionó, es muy importante compensar estas diferencias antes de proceder a suavizar la unión. Para este fin, usando las interpretaciones que se observan en el gráfico central, se calcularon correcciones estáticas para ambos conjuntos que minimizaron el efecto.
Al aplicar las correcciones, se obtiene las secciones de la derecha en donde se observa cómo han cambiado la forma de los marcadores para ambas, de manera de tender a igualarlas. Para poner de manifiesto el cambio ocurrido se transportaron los horizontes. Finalmente, en el extremo derecho se hace una representación conjunta de los datos corregidos, con Aguada Baguales en celeste transparente. Es muy claro el efecto positivo del procedimiento, que provocó coincidencia en la zona de superposición, que dejaron preparado los datos para sacar óptimos beneficios del proceso de Transición Suave. Por fuera de la zona de overlap, la corrección es nula. En la tabla 3 se muestran los datos a los que se necesitó aplicar este procedimiento.
Se implementaron transiciones graduales de un conjunto al otro en las zonas de superposición, aquellos segmentos de la transecta en donde hay dos versiones del dato sísmico. Esto se logró mediante la técnica de sumas pesadas, para lo que se genera una función de peso, que toma valor 1 en uno de los extremos de la zona de superposición y decrece linealmente hasta tomar valor cero en coincidencia con el extremo opuesto. En la figura 10 se esquematiza lo explicado, en donde se observa la función de peso y su inversa. La solución suavizada se obtiene haciendo
DATO 1 * (Función de peso) + DATO 2 * (Función de inversa)
El procedimiento se implementó en 13 de las 20 transiciones, y su efecto puede apreciarse en la comparación de la figura 11. Aquí volvemos a la zona de transición entre Loma La Lata y Aguada Baguales. A la izquierda se ve la unión antes de todas las correcciones y compensaciones descriptas, mientras que a la derecha se aprecia como luce después; la mejora experimentada no necesita explicación.
En la figura 12, a través de un detalle, se aprecia los beneficios de la secuencia de procesamiento, que incluye la interpolación en zonas de carencia y las mejoras en resolución (se explica más adelante), detalle y carácter.
La misma secuencia de procesamiento unificado descripta hasta este punto se aplicó a las sumas parciales cercana y lejana. Se puede decir que en este punto ya se contaría con un set de información apta para entrar a la inversión preapilado. Sin embargo, nada garantiza que las amplitudes de apilados cercanos, medios y lejanos guarden una relación consistente con curvas nominales de AVO acorde a la ecuación de Zoeppritz (Simm y Bacon, 2014). Asimismo, es esperable que, dada la procedencia variada de la información, además de la eventual distorsión de amplitudes, esta sea variable de dato a dato. Si pretendemos que la inversión preapilado produzca resultados que predigan el modelo terrestre con algún grado de fidelidad, debemos tomar recaudos para verificar y, eventualmente, corregir las relaciones de amplitud.
La estrategia de control y compensación de amplitudes adoptada hace máximo aprovechamiento de la abundante y bien distribuida información de pozos con información elástica completa (velocidades compresionales y de corte y densidad) y se puede resumir en los siguientes pasos, que en todos los casos usan como ventana temporal de análisis la asociada al sistema Quintuco-Vaca Muerta:
El punto 4 se beneficia por el hecho de contar con al menos un pozo por dato sísmico.
En la tabla 4 se puede ver para cada dato, de izquierda a derecha, las amplitudes nominales estimadas en los gathers sintéticos para los tres rangos de ángulo, las relaciones Medio/Cercano y Medio/Lejano estimadas con base en la sísmica real y en el extremo derecho los escalares calculados para que al aplicarlos a las sumas cercanas y lejanas llevan las relaciones de amplitud a las nominales observadas en los sintéticos.
En este punto es oportuno advertir que los datos de El Trapial, Curamched y Los Toldos (en rojo en “configuración recibida” de la tabla 1 y en celeste en la tabla 4) no tenían contenido suficiente de offset largos para producir una suma parcial de ángulo lejano robusta, dada la profundidad de Vaca Muerta. Se optó por aplicar un escalar sobre la base de los de las vecindades para obtenerlo por multiplicación desde el ángulo medio. Desde ya que esta estimación produce un dato incapaz de predecir anomalías en la relación Vp/Vs, pero cumple con el objetivo de completar el dato preapilado, y de dejar como legado un dato que pudiera oficiar de entrada a futuros estudios basados en la amplitud.
Otro tanto sucedió con la línea 2D HV89-202 (en azul en la tabla 1), que por no contar con información preapilado, se consideró al apilado recibido como suma parcial en ángulo medio, para luego generar Cercano y Lejano con base en las relaciones de amplitud de las vecindades.
Proceso e integración de datos del SSIS (Seismic Sequence Interpretation System)
En el ambiente asociado a la formación Vaca Muerta es de capital importancia, para la caracterización e interpretación, la definición de los límites de secuencias. Con el objetivo de robustecer y dar apoyo a esta interpretación se implementó el proceso llamado SSIS, que consiste en una interpretación semiautomática de un set de horizontes (o eventos) de muy alta densidad, en el caso de la transecta su número total supera los 5200. Como entrada se utilizó la amplitud de la transecta, tal como quedó luego la secuencia descripta en procesamiento sísmico: unificado en un solo survey y con todas las correcciones y compensaciones aplicadas. La interpolación que rellenó las zonas de carencia y le confirió continuidad que es clave para la performance del algoritmo de interpretación.
Cada evento u horizonte interpretado representa una línea (o superficie en caso de tratarse de un 3D) de tiempo geológico constante. No necesariamente son reflectores sísmicos, aunque en muchos casos pueden coincidir con estos; ejemplo de esto es una discordancia que se interpretará como una línea de tiempo, aunque no sea un reflector sísmico y a partir del cual se irán desprendiendo eventos (u horizontes), cada uno asociado, como se mencionó, a una línea de tiempo geológico constante. Del mismo modo, un límite de secuencia seguramente tendrá cambios de facies en su extensión lateral, por lo que la impedancia sísmica cambiará y, por ende, la reflexión sísmica también. De este modo, a partir de este denso set de eventos, se genera una sección como el gráfico obtenido de la Memoria 26 (Figura 13). En nuestro caso el eje vertical será el tiempo (sísmico) de reflexión convertido a profundidad, ya que la transecta se transformó a profundidad. Con lo cual, si bien los eventos u horizontes representan un tiempo geológico constante, cada uno se halla graficado a un tiempo relativo, ya que no se ha realizado la correspondencia necesaria para invertirla al tiempo geológico absoluto. Para realizar este paso sería necesario generar el diagrama de Wheeler y asignar tiempos geológicos absolutos a las secuencias analizadas.
En la figura 13, las líneas de tiempo constante de la parte superior son los eventos sísmicos interpretados por el SSIS, en la parte inferior se ha transformado al diagrama de Wheeler como se menciona en el informe.
En la figura 14 se muestra un detalle de la zona de Parva Negra con representación compuesta: los colores de fondo representan el módulo de Young, mientras que los horizontes en color verde y blanco son los calculados de modo semiautomático siguiendo los límites de secuencia de diverso orden. Son, por lo tanto, líneas tiempo geológico constante. La relación entre ellos revela diversas geometrías que ayudan a inferir los ambientes paleo-sedimentarios y, como consecuencia, las facies geológicas. En la figura 15 se aprecia la amplitud sísmica en alta resolución para el mismo sector de la transecta, pero en colores blanco y negro. Los horizontes están representados en color verde.
Inversión sísmica simultánea preapilado
Se construyó por interpolación de la información de 17 pozos según la inversa del cuadrado de la distancia y con la guía de 10 superficies: Quintuco, Tordillo y 8 secuencias interpretadas dentro del intervalo. De esta manera se obtuvieron tres volúmenes para impedancias P y S y densidad. La frecuencia se restringió a la banda de carencia sísmica estimada en 0-0-10-17 hz. Cabe destacar que, dado el carácter tortuoso de la traza de la transecta, la construcción del modelo se benefició por la geometría 3D. Implementar un pseudo 2D a lo largo del perfil habría sido más rápido y sencillo, pero hubiera introducido distorsión –severa en ciertas zonas–, debido a la mala evaluación de las distancias.
Se invirtió usando tres ondículas –una para cada suma parcial– que, como se ve en la figura 16, resultaron casi idénticas, tanto en espectro como en fase. Este hecho ofrece un diagnóstico positivo sobre el trabajo de procesamiento realizado. También son muy buenos síntomas la altísima correlación obtenida y la homogeneidad de predicción en todos los pozos, máxime teniendo en cuenta que la entrada es la unión de 19 datos.
Redes neurales
La sección de la figura 17 es la impedancia P producida mediante predicción por redes neurales. Esta revela resultados excelentes, de mayor detalle que la inversión y también muy homogéneos. Pone de manifiesto la capacidad de discriminación de las variaciones dentro del sistema Quintuco-Vaca Muerta. Debidamente interpretadas estas variaciones resultarán clave para entender el paleoambiente y su influencia en la roca como reservorio.
La figura 16 en el extremo derecho muestra el gráfico de correlación de la solución de la impedancia por redes neurales, que se obtuvo entrenando con los siguientes parámetros:
La correlación efectiva de entrenamiento resultó casi óptima (99%), pero más impactante es la correlación de validación superior a 96%. Recordemos que la validación se obtiene prediciendo todos los pozos, pero escondiendo uno a uno el que se está prediciendo, por lo que constituye una medida confiable de la probabilidad de éxito de una nueva locación. Por último, es de destacar la eficiencia del entrenamiento, teniendo en cuenta la variedad de procedencia del dato sísmico y la extensión de casi 400 km de norte a sur, ambos aspectos que tienden a comprometer la aplicabilidad del método.
Se volverá a recurrir a las redes neurales para refinar los atributos para la integración con geomecánica.
Interpolación en zonas de carencia de información sísmica
Si comparamos la figura 17 con el apilado de la figura 8, notaremos que ya no están presente los huecos sin dato, producto de la carencia de información. El método de interpolación implementado consistió en interpretar 10 niveles con gran detalle dentro de la ventana de interés, que luego asistirán a la propagación que completará las zonas sin información. La interpretación de detalle no se hizo guiada por la amplitud, sino por la impedancia de alta resolución proveniente de las redes neurales. Debido a que a nivel de dato invertido –no convolucionado– son mucho más evidentes los niveles estratigráficos y mucho mayor la certidumbre al conectar eventos a ambos lados de la zona de carencia. Una vez completada la interpolación, se procedió a propagar la impedancia a todo el hueco, confiriéndole gran continuidad a toda la sección de la transecta.
El procedimiento y sus resultados colaboran en dos aspectos:
Es muy difícil apreciar la calidad del resultado obtenido a la escala de la sección de la figura 17, por eso se suma la figura 18, que muestra en detalle el resultado de la interpolación entre Los Toldos y La Escalonada al norte, y rellena el hueco vinculado al flanco norte del Dorso de los Chihuidos. De izquierda a derecha, impedancia por redes neurales con interpretación auxiliar, impedancia interpolada y, en el extremo derecho, la amplitud obtenida por convolución con ondícula 8-12-90-100 hz, previo cálculo de la reflectividad a partir de la impedancia. Este último volumen, muy mejorado en contenido espectral y resolución, intervino en la interpretación estructural final de detalle y contribuyó como entrada a la clasificación de sismofacies, como se verá más adelante.
Conversión a profundidad
Para la exhibición de los resultados la conversión se hizo mediante Kriging con drift externo. Las velocidades de conversión calculadas por la relación marcadores estratigráficos/marcadores sísmicos en los 15 pozos oficiaron de dato duro, mientras que la velocidad compresional estimada por la inversión lo hizo como dato blando o drift externo.
Cálculo de propiedades de reservorio e integración al estudio geomecánico
A partir de las propiedades predichas en forma directa por la inversión –Impedancias P y S y densidad–, se estuvo en condiciones de calcular cualquier otra propiedad elástica de la roca mediante operación aritmética basada en relaciones proporcionadas por la teoría de la elasticidad.
Contar con los volúmenes de las propiedades permitió, a su vez, estimar propiedades que describen la roca como reservorio. Por medio de relación provenientes del estudio de los pozos y del estudio geomecánico, se derivaron el carbono orgánico total (COT) y el módulo de Young estático vertical.
A partir del estudio en los pozos se obtuvo una relación entre la impedancia P y el COT:
COT[w/w] = 0.338 - 5.3 x 10-5 x Zp + 2.07 x 10-9 x Zp2
El volumen de impedancia P se obtuvo por inversión y refinamiento por redes neurales y el COT queda en valores de fracción de peso. En la figura 19 se observa la sección de la transecta con la distribución de contenido orgánica obtenida por la relación mencionada y cuyo gráfico cruzado se muestra en la parte inferior.
La situación en esta predicción es algo más complicada, debido a la debilidad de la predicción sísmica de la densidad y el carácter vertical del módulo de Young.
La información sísmica con que contamos y sus derivados por inversión son de carácter isotrópico, por lo que no hay posibilidades de calcularlo por relaciones geomecánicas que, precisamente, requieren de propiedades anisotrópicas.
Entonces se diseñó una estrategia alternativa que combina el análisis elástico y nuevamente las redes neurales:
En la figura 20 se muestra un perfil en profundidad de un despliegue combinado de predicciones geofísicas con parámetros geomecánicos, como expresión de la integración lograda entra ambas comisiones. El background en color es la estimación, mediante la secuencia explicada, del Módulo de Young Vertical. Sobreimpuesto se despliega la estimación en todos los pozos, que exhibe una excelente correlación (98,2%). Adicionalmente, en los siete pozos comunes a los dos estudios, se incorporan las curvas estimadas por el estudio geomecánico (Sosa Massaro et al., 2022): esfuerzos in situ, presión de poro, más los diagramas de roseta con la distribución de fracturas.
Clasificación de sismofacies
Las sismofacies son volúmenes de categorías discretas, producto de la clusterización basada en la variación de la morfología de los atributos sísmico. Adecuadamente interpretadas y supervisadas por las electrofacies, constituyen la mejor herramienta para entender el paleoambiente y las tendencias sedimentarias que, junto con la diagénesis y demás fenómenos posdepositacionales, determinan las características que hoy exhibe la roca reservorio. Esto se explicará en más detalle en la sección “Integración”.
Para generarlas se usaron como volúmenes de entrada:
En los tres casos se clasificó en siete categorías y, previo al merge, se hizo una reasignación del código numérico para que el volumen resultante se vinculara en forma más natural con las categorías de las electrofacies. Esto facilitará la interpretabilidad y supervisión por parte de las electrofacies, verdaderas portadoras del significado sedimentológico y paleoambiental, ya desde su concepción.
En la figura 21 se muestra una sección en profundidad del volumen de sismofacies resultante, con el despliegue sobreimpuesto de las electrofacies en los 15 pozos.
Reconstrucción estratigráfica paleoambiental
Por último, una breve mención a lo que tal vez sea el producto de mayor valor agregado del proyecto: en la figura 22 se muestra una secuencia desde arriba hacia abajo, desde lo más antiguo a lo más reciente, que visibiliza ciclo por ciclo la evolución de la formación hasta llegar, en el panel inferior, a la distribución de facies sedimentarias tal como la interpretamos a partir de los resultados del trabajo regional. Este ejercicio ayuda a entender y relacionar aspectos, como la distribución de contenido orgánico y su vínculo con el ambiente de depósito, que permitió su preservación y posterior maduración en las zonas más profundas de la cuenca. También echa luz sobre la aparición de sedimentos de origen continental en el sur, hacia la Dorsal de Huincul. Este tipo de análisis es importante para sustentar y entender la zonificación observada en materia de potencial económico de la formación. El análisis paleogeográfico se inició a partir de la interpretación de un nivel del mar imaginario, ya que no se podría horizontalizar al tope de la formación Quintuco, debido a las pérdidas de material sedimentario con la erosión representada en la Discordancia Intravalanginiana. Para ello se utilizó como guía la predicción de la materia orgánica, cuya hipótesis es que las mayores concentraciones estaban relacionadas a las áreas más profundas (anóxicas) de la cuenca. La relación entre la interpretación de los ciclos, su morfología y el contenido de materia orgánica permitió establecer un horizonte de referencia como posible nivel del mar para el momento de la depositación de la formación Vaca Muerta. El objetivo de ello es poder entender la heterogeneidad de las facies sedimentarias de acuerdo con su posición en la cuenca y como fueron evolucionando en el tiempo en cada uno de los ciclos estratigráficos.
De acuerdo con el punto anterior se concluye que la secuencia muestra ser una alternativa válida cuando no se cuenta con datos de disparos o no se tiene presupuesto para una unión desde tiros de campo.
Además de lo descripto, de este trabajo se desprende un juego de datos preapilado “AVO compliant” que demostró su calidad y deja la puerta abierta a futuros trabajos de caracterización regional.
Los resultados se integraron muy bien a aquellos provenientes del estudio de Geomecánico. El carácter preapilado de la inversión posibilitó la predicción confiable de la relación Vp/Vs. Finalmente, las redes neurales aportaron para llegar a la predicción del Módulo de Young Vertical, clave para el estudio geomecánico. La sección lograda será usada como base en una de las gigantografías de la Exposición Mural CONEXPLO 2022.
Las sismofacies, alimentadas por toda la gama de productos de procedencia sísmica, mostraron claramente sentido estratigráfico y sedimentológico, claves para entender el ambiente depositacional y explicar los cambios observados en las propiedades de la formación.
La interpretación cuantitativa esbozada en este documento puede ser enriquecida e impulsa a revisitar los productos obtenidos, en la certeza de que se podrá llegar a conclusiones valiosas de índole regional.
La interpretación estratigráfica semiautomática se aplicó con confiabilidad y sus resultados colaboraron en la definición de los límites de secuencias.
Agradecemos a las siguientes empresas por la cesión de derechos de uso y exhibición de datos sísmicos y de pozo: YPF, Shell Argentina, Capex, TotalEnergies, Pan American Energy, Phoenix Global Resources, Oilstone, Geopark, Pampa Energía, Vista Wintershall, Tecpetrol, Chevron Argentina y Pluspetrol. Y también a las empresas que cedieron tiempo de su personal: JT Reservorios y Geoinfo.
González, G., Vallejo, D., Kietzmann, D., Marchal, D., Desjardins, P., González Tomassini, F., Gómez Rivarola, L. y Domínguez, R. F. (2016). Transecta regional de la Formación Vaca Muerta: Integración de sísmica, registros de pozos, coronas y afloramientos, IAPG-AGA, 244.
Hampson, D., Russell, B. y Bankhead, B. (2005). Simultaneous inversion of pre-stack seismic data, SEG Technical Program Expanded Abstracts: 1633-1637. December 07.
Herron, M., Matteson, A., (1993). Elemental Composition and Nuclear Parameters of Some Common Sedimentary Minerals. Nuclear Geophysics, 7(3), 383-406.
Hester, T., Schmoker, J., (1987). Determination of Organic Content from Formation-Density Logs, Devonian-Mississippian Woodford Shale, Anadarko Basin, Oklahoma. USGS Open-Fole Report 87-20.
Howell, J., Veiga, G., Spalletti, L. y Schwarz, E. (2015). The Neuquén Basin, Argentina - A Case Study in Sequence Stratigraphy and Basin Dynamics.
McDonough, K. J., Bouanga, E., Pierard, C., Sterne, E. J., Granath, J. W., Danforth, A. & Gross, J. S. (2012). Submarine Fan Chronostratigraphy From Wheeler-Transformed ION Basin SPAN Seismic Data, Late Cretaceous – Tertiary, Offshore Tanzania. PESA News Resources.
Mallet, J. L. (2004). Space-time mathematical framework for sedimentary geology: Mathematical Geology, 36, 1-32. 2008, Numerical earth models: European Association of Geoscientists and Engineers. Education tour series.
Mayer, C. y Sibbit, A. (1980). Global, a New Approach to Computer-Processed Log Interpretation, paper SPE 9341 presented at the SPE Annual Technical Conference and Exhibition, Dallas, September 21-24.
Qayyum, F., De Groot, P., Yasin, J. & Akhter, G. (2012). Building a Sequence Stratigraphic Framework from HorizonCube and Well Data. 74th EAGE Conference & Exhibition incorporating SPE EUROPEC, Copenhagen, Denmark.
Qayyum, F., Hemstra, N. & De Groot, P. (2013). Dense sets of seismic horizons - A New Approach to Stratigraphic Interpretation. Geophysical Corner, AAPG Explorer, 68-69.
Qayyum, F., Hemstra, N. & Singh, R. (2013). A Modern Approach to Build 3D Sequence Stratigraphic Framework. Oil and Gas Journal, 111(10).
Sánchez, M., Asurmendi, E. y Armas, P. Subgrupo Río Colorado Grupo Neuquén: Registros de paleosismicidad en la cuenca de antepaís andina, Cuenca Neuquina, Provincias de Neuquén y Río Negro.
Slatt, R., Abousleiman, Y. (2011). Multi-scale, Brittle-Ductile Couplets in Unconventional Gas Shales: Merging Sequence Stratigraphy and Geomechanics. Search and Discovery Article #80181. AAPG.
Simm, R., Bacon, M. (2014). Seismic Amplitude: An Interpreter’s Handbook, Cambridge University Press. Sosa Massaro, A., Frydman, M., Ezequiel Lombardo, E., Canatelli, A., Arguello, J., De Barrio, T., Marchal, D., Teran, C., Hryb, D., Díaz, E., Cruz, L., Paris, M. y Nawratill, A. Transecta geomecánica de la formación Vaca Muerta. Un nuevo proyecto de coopetición, 2022. A presentar en CONEXPLO 2022.
Wolak, J., Hemstra, N., Ochoa, J. & Pelissier, M. (2013). Reconstruction of depocenter evolution through time using relative stratigraphic thickness. The Leading Edge, 32(2), 172-177.
Varela, R., Hasbani, J., (2017). A rock mechanics laboratory characterization of Vaca Muerta formation. American Rock Mechanics Association, 16(2), 222-231.
Figura 1. Mapa de ubicación del proyecto (Modificado de Sánchez et al., 2004). Mapa base (Sísmica y Pozo).
Figura 2. Flujo de trabajo.
Figura. 3 Predicción y calibración del COT en los pozos que contaban con datos de laboratorio.
Figura 4. Relación entre carbono orgánico total (COT) y la impedancia P.
Figura 5. Análisis de clúster y clasificación de electrofacies a partir de la información de 15 pozos a lo largo de la transecta.
Figura 6. Correlación de las electrofacies y predicción de COT en los pozos a lo largo de la transecta.
Tabla 1. Información Sísmica Recibida. Parámetros.
Tabla. 2 Definición de la traza adoptada.
Tabla 3. Parámetros, valores y opciones adoptados en el tratamiento de los datos.
Figura 7. Conformación espectral.
Figura 8. Perfil en tiempo de la suma parcial de offset medio.
Figura 9. Ejemplo de corrección estática de overlap para la unión Loma La Lata-Aguada Baguales.
Figura10. Funciones de peso para “transición suave” en zonas de superposición.
Figura11. Mejoras experimentadas debido a la aplicación de la secuencia completa de correcciones y compensaciones en la zona de overlap entre Loma La Lata y Aguada Baguales.
Figura 12. Comparación de datos antes y después de la secuencia de procesamiento, incluidos la interpolación y el mejoramiento en frecuencia: izquierda 40-50 hz, derecha 90-100 hz.
Tabla 4. Compensación AVO compliant. Generación de escalares para compensar las amplitudes preapilado.
Figura 13. Esquema de interpretación (Memoria26).
Figura 14. Interpretación semiautomática sobre la sección de impedancia P.
Figura 15. Interpretación semiautomática sobre la amplitud en alta resolución por reflectividad y convolución (0-0-90-100 hz) en gama de grises.
Figura 16. Correlación para impedancias S y P por inversión e impedancia P mejorada por redes neurales.
Figura 17. Sección por los pozos de la impedancia P en su versión refinada por redes neurales.
Figura 18. Ejemplo de interpolación en la zona de carencia Los Toldos-La Escalonada.
Figura 19. Distribución de la materia orgánica a lo largo de la transecta.
Figura 20. Integración geofísica-geomecánica.
Figura 21. Sección de sismofacies con electrofacies y COT en los pozos.
Figura 22. Secuencia de reconstrucción paleoambiental por intervalos.
> SUMARIO DE NOTAS
“Un trabajo que mejoró el grado de entendimiento regional”
Estudio regional para la evaluación de la formación Vaca Muerta: generación de transecta geofísica :Integración con estudio geomecánico e interpretación secuencial semiautomática. Exposición mural de gigantografía Neuquén, cuenca Neuquina
Transecta geomecánica de la formación Vaca Muerta. Un nuevo proyecto de coopetición
Monitoreo microsísmico de superficie durante la fractura hidráulica de dos pozos horizontales en la formación Vaca Muerta
Desarrollo multilanding de la formación Vaca Muerta en el yacimiento La Amarga Chica, cuenca Neuquina. Experiencias, lecciones aprendidas y desafíos
Caracterización integral a escala de cuenca de la formación Vaca Muerta como reservorio no convencional
Caracterización morfológica y estructural de una megaestructura de tipo MTD (mass-transport deposits) de escala sísmica en el sistema Vaca Muerta - Quintuco, cuenca Neuquina
> Ver todas las notas