Modelo estocástico automático de fracción de arcilla en 3D a partir de datos de sondeo y sondeo tTEM

Scientific Reports volumen 12, número de artículo: 17112 (2022) Citar este artículo

752 Accesos

1 Citas

Detalles de métricas

En la mayoría de las zonas urbanizadas y agrícolas de Europa central, el subsuelo poco profundo está constituido por depósitos cuaternarios que suelen ser las capas más utilizadas (bombeo de agua, geotermia poco profunda, excavación de materiales). Todos estos depósitos a menudo están entrelazados de manera compleja, lo que conduce a una alta variabilidad espacial y una alta complejidad. Los datos geofísicos pueden ser una fuente rápida y fiable de información sobre el subsuelo. Aún así, la integración de estos datos puede llevar mucho tiempo, carece de una interpolación realista en un espacio 3D completo y la incertidumbre final a menudo no está representada. En este estudio, proponemos una nueva metodología para combinar perforaciones y datos geofísicos con incertidumbre en un marco automático. Una función traductora que varía espacialmente y que predice la fracción de arcilla a partir de la resistividad se invierte utilizando la descripción de los pozos como puntos de control. Se combina con un marco de interpolación estocástica 3D basado en un algoritmo de estadísticas de puntos múltiples y una función aleatoria gaussiana. Este novedoso flujo de trabajo permite incorporar de forma robusta los datos y su incertidumbre y requiere menos intervención del usuario que los flujos de trabajo ya existentes. La metodología se ilustra para datos electromagnéticos transitorios remolcados (tTEM) desde tierra y datos de pozos del alto valle del Aare, Suiza. En esta ubicación, se obtuvo un modelo 3D realista de alta resolución espacial de la fracción arcillosa sobre todo el valle. El conjunto de datos muy denso permitió demostrar la calidad de los valores predichos y sus correspondientes incertidumbres mediante validación cruzada.

Los acuíferos cuaternarios aluviales poco profundos se utilizan con frecuencia para el suministro de agua o la explotación de energía geotérmica a poca profundidad. En este contexto, a menudo es necesario abordar una amplia gama de cuestiones asociadas, como evaluar los recursos de aguas subterráneas, estudiar la migración de contaminantes, evaluar las interacciones con las aguas superficiales o evitar una superposición de zonas de influencia entre pozos geotérmicos vecinos. Todas estas preguntas sólo pueden responderse adecuadamente después de modelar la estructura y las heterogeneidades internas de esos acuíferos.

Estos modelos suelen construirse en varios pasos1,2. Para los modelos a pequeña escala, es común el uso de pozos como única fuente de datos. Sin embargo, este enfoque a menudo ignora la mayor parte de la heterogeneidad espacial y, por lo tanto, puede conducir a modelos inadecuados y conclusiones erróneas. Los sondeos son sólo una fuente de información para inferir la distribución local y vertical de las facies. A menudo son de ayuda limitada para estimar estructuras 3D complejas. Cuando el área de interés es amplia, aumentar el número de perforaciones para reducir la incertidumbre a un nivel aceptable suele ser difícil, lento y costoso. Una solución es combinar datos geofísicos y datos de pozos menos costosos. Los datos geofísicos a menudo se interpretan y combinan manualmente en un modelo estructural, que luego se completa con simulaciones litológicas o de facies y, finalmente, con propiedades físicas. Un flujo de trabajo de este tipo ha demostrado ser eficaz a la hora de generar modelos geológicos a partir de datos sísmicos o electromagnéticos aéreos, por ejemplo3,4. Pero estos pasos de modelado son complejos y a menudo necesitan un software diferente. Además, incluso si a menudo se utilizan algunos métodos estocásticos4,5, las incertidumbres no siempre se propagan a través de todo el flujo de trabajo. A menudo, algunos pasos se consideran deterministas y el modelo geológico resultante es el que se ajusta a la mayor parte del conocimiento integral disponible6,7,8,9. Por último, cuando se trabaja con un gran conjunto de datos, la construcción manual del modelo estructural utilizando tanto datos de sondeos como geofísicos puede llevar mucho tiempo.

Por lo tanto, existe la necesidad de un enfoque automático que sea capaz de integrar múltiples tipos de datos y producir modelos estructurales o paramétricos. Por ejemplo, la generación rápida de modelos de arcilla 3D con un algoritmo automático podría ser de gran beneficio para las autoridades locales que a menudo no tienen la capacidad de ejecutar la integración manual de los datos.

Sin embargo, en la mayoría de los casos, los principales datos disponibles a lo largo de los pozos son descripciones litológicas y no valores de resistividad o densidad. Pero estos parámetros físicos son los que infunden los datos geofísicos. Vincular directamente resistividades y litologías es difícil debido a la amplia variedad de factores que afectan la resistividad10 y a la descripción a menudo incompleta de las facies litológicas. La relación más fundamental es la Ley de Archie11, que vincula la resistividad con la saturación, la conductividad del agua, la tortuosidad y la porosidad. Esta relación empírica se basa en el supuesto de que la matriz no es conductora, supuesto que no es válido tan pronto como tenemos la presencia de minerales arcillosos que conducen corriente en su superficie. Esta conductividad de la superficie de los poros dependerá del área de la superficie, el tamaño del grano, el tipo de arcilla y el contenido de arcilla. La estimación de todos estos parámetros que pueden tener una alta variabilidad espacial es complicada. Una revisión reciente12 señala también problemas de escala. La mayoría de las leyes empíricas de laboratorio se miden a escala central, donde la muestra está en el rango de una docena de centímetros cúbicos. Y la ampliación de dichas leyes a la escala de campo no es sencilla. Finalmente, la mayoría de las descripciones litológicas asociadas a los sondeos son cualitativas y no cuantitativas12. Además, es demasiado simplista aplicar una función que vincule la descripción litológica a una única resistividad. Cada descripción litológica puede asociarse con una amplia gama de resistividad con cierta superposición entre diferentes litologías13.

Una respuesta a esa cuestión es definir la probabilidad de tener una facies litológica determinada condicionada al valor de resistividad. Estas distribuciones de probabilidad pueden estimarse a partir de un muestreo de pozos y modelos de resistividad coincidentes13. Sin embargo, este enfoque no tiene en cuenta la distribución espacial de los pozos e ignora la posible dependencia espacial entre el tipo de sedimento y su resistividad. La función de distribución de probabilidad (PDF) se calcula a partir de pozos y modelos de resistividad en todo el dominio. Si una facies litológica determinada es siempre más resistiva en una subárea del dominio, la PDF calculada para todo el dominio no reflejará esto.

Para resolver ese problema, en lugar de intentar estimar una única FDP, Foged et al.14 propusieron un método basado en la inversión de una función traductora que varía espacialmente entre la resistividad y la CF. La función proporciona el mejor ajuste entre el CF observado en los pozos y el CF calculado a partir de los modelos de resistividad. Este método tiene la ventaja de no depender de ninguna estimación previa de parámetros y sólo se infiere a partir de datos observados. También tiene la ventaja de transformar la descripción litológica en una variable continua, haciendo posible el escalado, teniendo en cuenta la coubicación de la función. Sin embargo, incluso si este método puede estimar el CF en la posición donde existen modelos de resistividad geofísica, aún necesita ser interpolado para obtener un modelo continuo 3D completo.

En este sentido, tras predecir el FC, Vilhelmsen et al.15 utilizaron un enfoque geoestadístico. Existen varios métodos geoestadísticos16,17,18 que son capaces de interpolar datos y producir modelos y simulaciones realistas. Su uso está muy extendido en disciplinas como la evaluación de riesgos, la gestión de recursos, la minería, la ingeniería petrolera o la modelización geológica. La estadística de puntos múltiples19,20 (MPS) es una de estas técnicas geoestadísticas. Es un método no paramétrico que se basa en el uso de una Imagen de Entrenamiento (TI) para inferir la variabilidad espacial de una o múltiples variables. MPS ha demostrado ser capaz de generar una variabilidad espacial compleja y realista en una amplia variedad de situaciones21,22,23,24,25. Vilhelmsen et al.15 agruparon el CF en unidades y lo utilizaron como TI en un procedimiento MPS. Los Datos Duros se definieron como las zonas donde el cluster es más seguro, principalmente en resistividades muy bajas y altas. Las otras áreas sólo están limitadas por datos blandos. Esto permite reflejar un tipo de incertidumbre en los datos de CF y mostrar variaciones en zonas inciertas del modelo. Pero este método requiere elegir el valor umbral discreto en el que la pertenencia a un grupo se vuelve segura. Además, no permite la propagación de la incertidumbre a partir de los datos. La incertidumbre sobre la resistividad o la CF se pierde después de la agrupación. Por lo tanto, un valor que estaría cerca de los límites del grupo pero asociado con una gran incertidumbre se considerará tan seguro como un valor que entraría en un grupo más allá de toda duda. Además, debido a la forma de la función traductora, pequeños cambios en el valor de resistividad en la zona de transición entre arcilla y no arcilla pueden tener enormes impactos en el CF estimado. Finalmente, dado que la incertidumbre de los datos se utiliza como normalización en la función objetivo de la inversión CF, sostenemos que también debe considerarse al aplicar la función e interpolar los resultados. Esto muestra la limitación de utilizar un TI determinista cuando se debe tener en cuenta la incertidumbre.

En este artículo, proponemos un flujo de trabajo extendido para generar automáticamente un modelo de fracción de arcilla 3D, con una sólida propagación de la incertidumbre desde los datos hasta el modelo final. El documento está estructurado de la siguiente manera. Primero presentamos los tres pasos principales de la metodología: la inversión determinista CF, el marco de interpolación estocástica y, finalmente, la implementación de la validación cruzada. Luego, presentamos la aplicación de la metodología para un conjunto de datos geofísicos electromagnéticos transitorios (tTEM) remolcados desde tierra26 adquiridos en el valle superior del Aare, Suiza.

La metodología propuesta permite generar automáticamente un modelo CF 3D de un valle cuaternario completo. Un aspecto clave de este enfoque es confiar en los datos mismos para inferir la mayoría de los parámetros automáticamente. Las siguientes subsecciones de este documento describirán en detalle todos los pasos. Pero primero proporcionemos una descripción general rápida del flujo de trabajo completo como se describe en la Fig. 1.

Resumen de los principales pasos de la metodología propuesta. El área sombreada resalta la parte estocástica del flujo de trabajo.

Las entradas son, por un lado, un denso conjunto de datos de registros de resistividad obtenidos invirtiendo las mediciones geofísicas (aquí datos tTEM) y, por otro lado, un conjunto de datos mucho más escaso de registros geológicos que pueden usarse para estimar la referencia. CF a lo largo de los pozos. Para producir un modelo CF 3D en todo el dominio, primero invertimos una función traductora CF (Paso 1 en la Fig. 1) utilizando el método de Foged et al.14 y ambos conjuntos de datos. Luego, la función de traducción se utiliza para estimar el CF a partir de la resistividad de todos los sondeos geofísicos disponibles (incluidos los lugares donde no hay datos de pozo disponibles). La principal novedad del flujo de trabajo propuesto es la metodología empleada para interpolar los valores de CF resultantes en 3D (Paso 2 en la Fig. 1). Respecto a trabajos anteriores, hay dos mejoras importantes. Una es evitar tener que clasificar los datos; los valores de CF se utilizan directamente. La segunda, y quizás la más importante, es que consideramos las diversas fuentes de incertidumbre y las propagamos en todo el flujo de trabajo. Para ello, la interpolación se realiza con un algoritmo de Estadística de Puntos Múltiples, utilizando los propios datos de CF en 3D como Imagen de Entrenamiento (TI). La densa cobertura geofísica permite utilizar directamente los datos sin tener que añadir información externa. Para tener en cuenta las incertidumbres tanto en los pozos como en los datos tTEM, utilizamos funciones aleatorias gaussianas para simular mapas de error, que se incluyen en los datos de acondicionamiento y las imágenes TI antes de las simulaciones MPS. Finalmente (Paso 3 en la Fig. 1), se realiza una validación cruzada para comprobar la calidad de los resultados. Una parte de los datos de acondicionamiento no se utilizan durante la interpolación y se comparan con los resultados. Luego evaluamos cómo funciona el método tanto en términos de valor predicho como de incertidumbre predicha.

Esta parte de nuestra metodología (Paso 1 en la Fig. 1) sigue de cerca el trabajo de Christiansen et al.27 ampliado por Foged et al.14. Aplicamos su idea de invertir una función que llena el vacío entre resistividad y CF. Dado que el método se ha descrito en detalle en estas publicaciones anteriores, aquí sólo se presenta un breve resumen. Este enfoque se desarrolló por primera vez para evaluar el riesgo de contaminación por nitratos27 y se denominó Espesor de arcilla acumulada (ACT). Posteriormente se propuso una extensión para estimar el FC con la misma metodología14. El supuesto principal es que los principales cambios en la resistividad en depósitos saturados y no consolidados son causados ​​por variaciones en la cantidad de arcilla. Teniendo esto en cuenta, se puede construir una función traductora que vincule la resistividad y el CF.

Sin embargo, es necesario discutir algunos desafíos. Primero, los modelos de resistividad derivados de los datos geofísicos son siempre una versión suavizada de la realidad, debido a la huella del instrumento, la capacidad de resolución y los procedimientos de inversión. En segundo lugar, la composición del propio sedimento arcilloso puede cambiar y provocar diferentes valores de resistividad en diferentes lugares con descripciones litológicas similares. Esto implica que una función de traducción basada en el uso únicamente de resistividad como entrada no puede representar adecuadamente el vínculo entre resistividad y conductividad hidráulica.

(a) Ilustración del ajuste automático de la función de traducción en un punto de datos. (b) Ejemplo de un registro litológico del componente principal del USCS. (c) CF estimado asociado en una cuadrícula de capas regulares (\(CF_{log}\)), derivado de la descripción litológica del pozo. La línea es la CF predicha utilizando los modelos de resistividad y las funciones de traducción invertida. (\(CF_{rho}\)) (d) Modelo de resistividad a partir de la inversión suave de los datos tTEM en la posición del pozo. Para cuatro capas, se muestra la función de traducción invertida utilizada.

Por lo tanto, Foged et al.14 propusieron un procedimiento basado en la inversión de dos parámetros de una función traductora que predice un CF a partir de un valor de resistividad y permite que estos parámetros varíen en el espacio. En otras palabras, dos valores de resistividad idénticos y alejados entre sí no necesariamente corresponderán al mismo valor de CF. Los datos de entrada básicos, en este procedimiento, son los registros geológicos de los pozos, como se ilustra en la Fig. 2b, y los perfiles de resistividad geofísica en las mismas ubicaciones, como el que se muestra como una curva roja en la Fig. 2d. Los dos conjuntos de datos se adquirieron utilizando métodos independientes y están ubicados juntos.

Con base en la descripción geológica, la fracción de arcilla \(CF_{log}\) a lo largo del pozo se estima dividiendo el registro geológico en intervalos de profundidad regulares y calculando la proporción de arcilla litológica para cada intervalo. Debido a que la descripción geológica es a menudo cualitativa, el CF resultante es incierto y el resultado es un valor medio y un rango de CF para cada intervalo. Este resultado se ilustra como cuadros grises en la Fig. 2c.

El siguiente paso es definir una función traductora parametrizada:

donde erfc es la función de error complementaria, \(\rho\) es la resistividad, y \(m_{up}\) y \(m_{low}\) son los valores de resistividad en los que la función devuelve un peso de 0,975 y 0,025 respectivamente. Los parámetros \(m_{up}\) y \(m_{low}\) se identifican durante la inversión y pueden considerarse como los límites de resistividad solo para arcilla/solo arena. La Figura 2a ilustra la inversión de \(m_{up}\) y \(m_{low}\) para un punto de datos. Antes de la inversión, el \(CF_{rho}\) predicho no está de acuerdo con la cantidad observada de arcilla \(CF_{log}\) descrita en el pozo colocado, lo que resulta en un desajuste. La inversión ajustará \(m_{up}\) y \(m_{low}\) para reducir el desajuste.

El ejemplo anterior muestra el principio para un único punto de datos, pero el problema está correlacionado en 3D y apunta a minimizar un desajuste global. Dos funciones de traducción vecinas no pueden tener parámetros drásticamente diferentes. La Figura 2d muestra, por ejemplo, cómo la función del traductor variará con la profundidad de un solo pozo. La función se describe para cuatro profundidades diferentes en la Fig. 2d, pero está definida para todas las capas.

Para obtener este tipo de resultado, el procedimiento de inversión sigue los siguientes pasos. Primero, el residuo de datos global se define como el error cuadrático promedio entre el CF predicho (\(CF_{rho}\)) y el CF medido en los pozos (\(CF_{log}\)):

Los datos residuales están normalizados por \(\sigma ^2\), la varianza combinada de \(CF_{rho}\) y \(CF_{log}\), y N el número de datos condicionantes. Luego se agregan parámetros de regularización para garantizar que los puntos vecinos no muestren variaciones bruscas y poco realistas. La regularización limitará la variación espacial de \(m_{up}\) y \(m_{low}\). La regularización se expresa como:

donde \(N_{con}\) es el número de pares de parámetros, A es la diferencia en el espacio logarítmico entre los dos pares de parámetros investigados y \(e_i\) es un factor dependiente de la distancia. Cuanto más alejados estén dos puntos, mayor puede ser su diferencia. Finalmente, la función objetivo completa es:

Los parámetros óptimos (\(m_{up}\) y \(m_{low}\)) se obtienen minimizando Q usando un esquema iterativo de Gauss-Newton con una modificación de Marquardt27. El ajuste entre el \(CF_{rho}\) y el \(CF_{log}\) predicho no es perfecto, ya que los parámetros de la función traductora se ven afectados no sólo por los datos coubicados (modelos de resistividad del pozo) sino también por todos los vecinos (a través de la regularización). La inversión busca un mínimo global.

En el ejemplo mostrado en 2, podemos indicar cómo la función traductora se adapta para identificar estructuras que no están bien representadas en el modelo de resistividad, como la capa de arcilla poco profunda.

Al final del paso de estimación de la fracción de arcilla, cada registro de resistividad se asocia con un registro \(CF_{rho}\) estimado, obtenido aplicando la función traductora con los parámetros óptimos \(m_{up}\) y \(m_ {bajo}\).

El último paso es propagar la incertidumbre de los valores de resistividad al CF. De hecho, las resistividades se obtuvieron mediante una inversión geofísica que es capaz de estimar la incertidumbre sobre las resistividades. Luego, las funciones de traducción se pueden aplicar no solo a los datos de resistividad sino también a sus incertidumbres, lo que da como resultado un \(CF_{rho}\) en los puntos de adquisición de tTEM con incertidumbre. Para el siguiente paso de simulación, el \(CF_{log}\) se utilizará como datos condicionantes para las celdas en las que tenemos información del pozo. \(CF_{rho}\) se utilizará en otros lugares.

Para generar un modelo 3D completo de CF y resistividad, los valores obtenidos en el paso anterior deben interpolarse para cubrir el espacio donde no se han adquirido datos geofísicos o de pozo. Esta parte del flujo de trabajo corresponde al Paso 2 en la Fig. 1 y es la principal parte novedosa del flujo de trabajo propuesto. Utilizamos el método MPS de muestreo directo28 en este proceso. La principal ventaja de utilizar un enfoque MPS es que puede aprender automáticamente los patrones espaciales de las estructuras regionales a partir del conjunto de datos muy denso construido en el paso anterior y utilizarlo para representar la variabilidad espacial y la incertidumbre en las áreas interpoladas. Sin embargo, la aplicación estándar de técnicas MPS supone que el TI y los datos duros (HD) son deterministas. Aquí proponemos un método para ir un paso más allá y tener en cuenta la incertidumbre en estos datos de entrada. Por lo tanto, incluimos datos duros (HD) e imágenes de entrenamiento (TI) no deterministas en el algoritmo MPS. El método general, como se muestra en el Paso 2 en la Fig. 1, consistirá en aplicar el algoritmo MPS muchas veces con diferentes TI y diferentes HD para generar un conjunto de modelos 3D interpolados de FQ. A partir de este conjunto de simulaciones, derivaremos distribuciones de probabilidad para la CF y la resistividad en cualquier ubicación en el dominio 3D.

El principio general de los algoritmos MPS es llenar una cuadrícula de simulación de forma iterativa mientras se reproducen los patrones del TI. El algoritmo de muestreo directo28 es un algoritmo MPS versátil basado en un remuestreo aleatorio y condicional del TI. En este artículo, utilizamos la implementación de Deesse29 del muestreo directo. Para generar una simulación, primero se llena la cuadrícula 3D con los datos físicos disponibles. Luego, el algoritmo visita aleatoriamente todas las ubicaciones restantes de la cuadrícula. Para simular un valor en una ubicación determinada, se extraen los n HD más cercanos y los valores ya simulados para definir un patrón de datos. Luego, el patrón de datos se utiliza para buscar ubicaciones que tengan patrones similares en el TI. Durante esa búsqueda, se calcula la distancia entre el patrón de datos y los patrones encontrados en el TI. Si la distancia está por debajo de un valor umbral (t), entonces los dos patrones se consideran similares y el valor de píxel en la ubicación faltante se toma del TI y se copia en la cuadrícula de simulación. Para acelerar el algoritmo y evitar copiar y pegar directamente el TI, solo se escanea una fracción (f) del TI. Los tres parámetros n, t y f los elige de antemano el usuario. Meerschman et al.30 ofrecen algunas recomendaciones prácticas para la selección de esos parámetros. Uno de los puntos fuertes del algoritmo de muestreo directo es que puede tratar variables continuas y múltiples simultáneamente. Significa que se puede proporcionar un TI que contenga varias variables y el algoritmo puede simular una o varias variables condicionadas a los demás datos disponibles. Esto se describe en detalle en varios artículos28,31. Se puede utilizar para describir la presencia de tendencias en la imagen de entrenamiento y en la cuadrícula de interpolación utilizando variables auxiliares22,24.

Para interpolar el CF y la resistividad, utilizamos DeeSse y diseñamos un problema de simulación de cuatro variables: dos variables auxiliares y dos variables principales. Cada una de ellas tiene un número específico de vecinos n y un umbral t, mientras que la fracción escaneada f es común a todas las variables. Las dos variables secundarias son (1) la profundidad de la celda y (2) el norte de la celda. Su propósito es describir nuestro conocimiento previo de las tendencias espaciales. Esta elección refleja nuestra expectativa de encontrar ciertos patrones más preferiblemente en una profundidad o área determinada en nuestro dominio interpolado. Por ejemplo, las estructuras más profundas del TI tenderán a reproducirse más profundamente en la simulación. Sin embargo, los valores de umbral que seleccionamos son altos y n es pequeño, lo que limita solo parcialmente el algoritmo. También elegimos incluir solo el norte, ya que esperamos tener la mayoría de las variaciones del patrón a lo largo del valle (orientado aproximadamente de N a S). Esta suposición se basa en una inspección visual de los modelos de resistividad, el conocimiento geológico y los pozos. Las principales variables de simulación son (3) la transformada logarítmica de la resistividad y (4) la fracción de arcilla. Para estas variables, elegimos un mayor número de vecinos n y un umbral menor t ya que estas son las variables de interés. Se simulan las cuatro variables. También activamos la opción Pirámide Gaussiana del código DeeSse29. Al hacerlo, los patrones espaciales de los datos se analizan y modelan en múltiples escalas. El algoritmo utiliza filtros gaussianos para crear una pirámide de imágenes de escala más gruesa coubicadas que se utilizan conjuntamente para entrenamiento y simulación32. Este método mejora la calidad de la simulación entre áreas densamente cubiertas y áreas más dispersas, al tiempo que puede ser sensible a diferentes escalas de variación.

Para ejecutar el algoritmo MPS e interpolar el CF, necesitamos proporcionar un TI y HD 3D. Los HD son simplemente los datos puntuales 3D derivados del Paso 1 en los sondeos y sondeos geofísicos. Debido a que la densidad espacial de los datos tTEM es muy alta (consulte la aplicación de ejemplo a continuación), este conjunto de datos suele ser muy denso y puede usarse directamente como TI. Esta situación corresponde al llamado problema de llenado de huecos, en el que necesitamos interpolar sólo algunas partes de un conjunto de datos que ya es muy denso. El MPS y, en particular, el muestreo directo han demostrado ser muy eficientes para estos problemas31,33. En estos casos se utilizan los mismos datos que HD y TI. El supuesto detrás de esta decisión de modelado es que la cobertura de los datos es suficiente para representar adecuadamente las estadísticas espaciales de las variables que deben interpolarse.

En la práctica, los puntos 3D HD se colocan en la cuadrícula de simulación 3D por adelantado. Estos puntos son estimaciones de los valores de CF y resistividad provenientes de los pozos o de la aplicación de la función de traducción en los datos tTEM. Las variables auxiliares (profundidad y norte) también se calculan y almacenan para cada ubicación en la cuadrícula. En este punto de la metodología, podríamos simplemente aplicar el algoritmo MPS y obtener un conjunto de simulaciones estocásticas que representen la incertidumbre debida a la interpolación. Esto es lo que normalmente se hace cuando se utiliza MPS. El TI y el HD son deterministas.

Pero los valores de CF y resistividad tienen una incertidumbre asociada que ya se estimó en el paso anterior del flujo de trabajo. Para tener en cuenta también esta fuente de incertidumbre, una posibilidad sería considerar solo los puntos más ciertos como HD, por ejemplo, las resistividades extremadamente bajas o altas, que casi con certeza están asociadas con puntos enteramente arcillosos o no arcillosos, respectivamente. Pero esta idea tiene dos desventajas principales. Primero, se debe tomar una decisión para determinar los límites de resistividad superior e inferior. Esto iría en contra de la idea de implementar el procedimiento lo más automatizado posible. En segundo lugar, la distribución de probabilidad de los valores del TI debe ser similar a la distribución del HD. Si son diferentes, las simulaciones pueden tender a sobrerrepresentar los puntos arcillosos y no arcillosos en la simulación. En términos más generales, no estaríamos utilizando todos los HD disponibles y podríamos asignar una resistividad o un CF a una celda completamente fuera del rango de incertidumbre derivado de las mediciones de campo. Por lo tanto, necesitamos una mejor manera de dar cuenta de la incertidumbre sobre estos datos.

Superamos estos desafíos mediante el uso combinado de MPS y la función aleatoria gaussiana (GRF), que nos permite realizar MPS con una imagen de entrenamiento (TI) y un conjunto de datos duros (HD) que no son deterministas. Los modelos GRF son bien conocidos18,34. La variabilidad espacial se modela utilizando distribuciones paramétricas multigaussianas. Estos modelos se definen con un modelo de covarianza o variograma que representa la variabilidad espacial. Se pueden generar múltiples realizaciones, con o sin datos condicionantes. En la metodología propuesta, el modelo GRF se utiliza para representar la incertidumbre de las mediciones de los datos de CF y resistividad. En la práctica, para cada simulación MPS de CF y resistividad, se generan un TI y un HD diferentes.

Principio para la generación de un conjunto de datos duros e imágenes de entrenamiento para tener en cuenta la incertidumbre en HD.

La Figura 3 esboza el principio general de generación de los TI. El HD (y el TI) originales se ven ligeramente perturbados al agregar algo de ruido dentro del rango de incertidumbre estimado para estos puntos de datos. La simulación del ruido se realiza mediante simulaciones 3D GRF incondicionales. Las simulaciones de GRF tienen una media de 0 y una varianza de 1. Para tener en cuenta la posible correlación espacial del ruido, se ajusta automáticamente un variograma anisotrópico 3D a los datos de TI. Utilizamos un algoritmo reflectante de región de confianza35 para optimizar el umbral y los rangos multidireccionales de dos contribuciones (una gaussiana y una exponencial). Luego, el alféizar se normaliza para obtener una varianza de 1. Los valores de CF y resistividad para 3D TI (y 3D HD) en cada iteración se definirán sumando el HD con un ruido reescalado:

siendo \(\text {TI}_{i}(x,y,z)\) el valor de TI para la simulación i en la posición (x, y, z) en 3D, \(\text {data}(x siendo ,y,z)\) el valor estimado determinista de la variable de interés en esa ubicación, \(\sigma (x,y,z)\) la desviación estándar de la incertidumbre estimada y \(GRF_{i}( siendo x,y,z)\) el valor GRF simulado incondicional en la posición (x, y, z) para la iteración i. El valor final simulado (CF por ejemplo) varía dentro del rango de su incertidumbre. En el ejemplo presentado en la Fig. 3, el TI inicial derivado del CF estimado se modifica ligeramente agregando un ruido aleatorio correlacionado que depende de las incertidumbres locales. Se realiza la misma operación para los conjuntos de datos geofísicos y de pozo utilizando el mismo modelo GRF.

En resumen, se genera un TI 3D único y su correspondiente conjunto 3D HD para cada simulación y se proporciona como entrada al algoritmo de muestreo directo que simulará las partes faltantes en la cuadrícula 3D. Esto dará como resultado un conjunto de realizaciones en 3D. Luego podemos acumular todas las simulaciones y calcular la media y la desviación estándar para cualquier ubicación. Las celdas en las que tenemos datos condicionantes tendrán la misma media y desviación estándar que los datos condicionantes originales.

Para verificar el desempeño de la metodología propuesta, implementamos un paso de validación cruzada. La validación cruzada permite cuantificar los errores asociados a un modelo. El principio es simple: se crea un subconjunto de datos; contiene una parte aleatoria del conjunto de datos original. Luego se aplica la interpolación estocástica utilizando solo el subconjunto de datos como HD y la simulación resultante se compara con los datos excluidos. Se pueden utilizar varios indicadores de error36. En nuestro caso, un muestreo aleatorio de los densos datos de CF no crea brechas suficientemente grandes para producir una estimación del error representativa. Los puntos que faltan están demasiado restringidos por los vecinos. Para crear una interrupción mayor en los datos, asignamos a cada punto un grupo, basado en una agrupación 3D de k-medias37. La agrupación se realiza utilizando las coordenadas espaciales (x,y,z) del modelo. El propósito de este paso es crear múltiples grupos espaciales de tamaños similares que serán excluidos aleatoriamente de las simulaciones. En nuestro caso, definimos 28 grupos. Luego se compara el valor simulado y la desviación estándar de estas zonas con el valor real. El error y el error normalizado para cada punto del modelo se definen como

siendo \(\varepsilon\) y \(\varepsilon _{norm}\) respectivamente el error y el error normalizado, n el número total de simulaciones, verdadero es el punto de datos no incluido en el conjunto de datos interpolados, \(sim_i\ ) es el valor simulado y sigma es la desviación estándar de \(sim_i\) sobre las n simulaciones que no incluyeron estos puntos de datos. Estos dos indicadores se calculan puntualmente. El error normalizado es un indicador importante ya que muestra qué tan bien estimamos la incertidumbre, siendo el objetivo una relación cercana a 1.

(a) Distribución de la calidad y profundidad de las capas después de la normalización. (b) Descripción general de la cobertura de los datos tTEM y mapa de situación general (creado con QGIS V3.22.9, qgis.org. El mapa base está disponible gratuitamente en la Oficina Federal Suiza de Topografía swisstopo).

El área de estudio está situada en el centro de Suiza y cubre una sección de aproximadamente 20 km del valle del Alto Aare entre las ciudades de Thun y Berna. El valle del Alto Aare presenta una geología cuaternaria típica de los valles alpinos: el basamento suele tener unos cientos de metros de profundidad y está cubierto por un complejo entrelazamiento de depósitos glaciares, fluvioglaciares y fluviales38,39. Se han identificado múltiples ciclos de avances y retrocesos glaciales en la cuenca suiza, provocando múltiples cambios en los procesos de deposición40. En el caso particular del Alto Valle del Aar, se han identificado al menos cuatro ciclos glaciales en las perforaciones41. Sin embargo, una descripción completa de la litoestratigrafía basada únicamente en descripciones de pozos es casi imposible porque depósitos similares de diferentes edades pueden estar superpuestos y entrelazados. Un acuífero superficial constituido por el depósito aluvial más reciente está presente en todo el valle. El nivel freático se encuentra entre unas pocas docenas de centímetros y hasta 2 metros por debajo de la superficie. El conocimiento general del acuífero sugiere que el espesor del acuífero superior oscila entre 4 y 20 m dependiendo de la zona. Además, se ha identificado un segundo acuífero más profundo en algunas perforaciones profundas separadas del superficial por una capa de arcilla. La extensión, conectividad y espesores exactos de ambos acuíferos no se conocen ampliamente. Sin embargo, en la zona se utilizan 6.350 pozos de bombeo (geotérmicos poco profundos o de agua potable) y 5.300 pozos de inyección (reinyección de agua después de una bomba de calor geotérmica). En este contexto, la realización de un modelo 3D sería de gran beneficio para las autoridades locales para evaluar la vulnerabilidad del acuífero superior. Debido a la altura del nivel freático, se pueden asumir condiciones de saturación en casi toda la altura del modelo.

En la zona de interés se encuentran descritos litológicamente 1542 sondeos. El alto valle del Aare también fue uno de los lugares de prueba elegidos para el proyecto Geoquat del Instituto Topográfico Nacional Suizo39. En este contexto, realizaron la digitalización y estandarización de los datos del pozo. Todas las descripciones de los pozos se convirtieron a descripciones estándar del USCS42. Además, se añadió un valor de control de calidad para cada capa que evalúa la confiabilidad de la información geológica: a cada capa se le asignó una calificación del uno al cinco, dependiendo de la precisión de la descripción. uno corresponde a una capa con sólo una descripción básica cuando cinco corresponde a una capa donde se han realizado mediciones de laboratorio y se encuentra disponible una descripción litológica multifase completa. Se realizó una estimación del contenido de arcilla en los pozos utilizando las pautas del USCS y la incertidumbre se escaló de acuerdo con la ley. Una capa de arcilla mal descrita (Calidad 1) tendrá, por ejemplo, un valor de CF de \(1 \pm 0,5\). Consideramos que hasta el 50% de la capa arcillosa podría ser material no arcilloso. Por otro lado, uno bien descrito (Calidad 5) tendrá un valor de CF de \(1 \pm 0.08\). No se tienen en cuenta las capas faltantes, artificiales o no descritas.

La Figura 4a muestra la distribución de la calidad de las capas y su profundidad asociada. Dado que la mayor parte de la exploración de pozos se realiza con fines geotérmicos poco profundos o geotécnicos, observamos que alrededor del 60% de los datos se encuentran por encima de los 10 m de profundidad y que la proporción sigue aproximadamente una distribución de ley de potencia decreciente con la profundidad.

En enero de 2020, se llevó a cabo un estudio electromagnético transitorio remolcado (tTEM) desde tierra en el valle superior del Aare. Los datos se invirtieron utilizando diferentes regularizaciones (nítidas y suaves) y los modelos de resistividad resultantes se utilizan en este estudio. Todo el flujo de trabajo de procesamiento, los detalles sobre la regularización y los datos se describen en Neven et al.26. El conjunto de datos cubre alrededor de 1500 hectáreas, con un espacio entre líneas de 25 metros. La frecuencia de muestreo después del procesamiento es de aproximadamente 1 modelo de resistividad cada 10 metros a lo largo de las líneas. La inversión de los datos se realizó utilizando el código de inversión AarhusInv (https://hgg.au.dk/software/aarhusinv)43. El residuo promedio de la inversión es 0,52 para la regularización aguda, lo que significa que tendemos a tener un ajuste excelente entre los datos predichos de nuestros modelos de resistividad y las mediciones de campo. Para estimar la profundidad de la investigación, utilizamos la matriz de sensibilidad jacobiana de la última iteración44. Al hacerlo, podemos identificar la profundidad exacta a la que cada modelo de resistividad está mal representado en los datos registrados. Por debajo de esta profundidad, los modelos de resistividad están cegados. Nuevamente, para obtener una descripción completa de la inversión del conjunto de datos y el control de calidad, se remite a los lectores a Neven et al.26.

La inversión de la función de traducción se realizó en el conjunto de datos de inversión brusca tTEM20AAR26. Se tuvieron en cuenta 57.862 modelos de resistividad de 30 capas y se cegaron a la profundidad de investigación estándar. La Figura 4b describe la cobertura de datos geofísicos en el valle después del procesamiento. La mayoría de los campos fuera de las ciudades están mapeados. Las ciudades, así como los alrededores de las carreteras y las vías del tren, son claramente visibles, ya que el acoplamiento electromagnético en tales entornos prohíbe el uso de un método inductivo y se deja al descubierto. La incertidumbre sobre la resistividad se estimó utilizando la matriz de covarianza de la última iteración45. Incluso si este método no reemplaza completamente la incertidumbre que se puede derivar de una inversión estocástica, tiene la ventaja de poder reflejar la incertidumbre relativa en las celdas del modelo y al mismo tiempo es relativamente rápido de calcular.

Resultados de la función de traductor invertido en 50.000 puntos TEM. Además, la incertidumbre promediada se calculó con una media móvil de 0,05 CF de ancho.

El desajuste de datos promedio [ver Ec. (2)], o residual, de la aplicación de la función traductora a los modelos geofísicos es 0,38. Significa que la diferencia entre \(CF_{rho}\) y \(CF_{log}\) predichos está más de dos veces dentro del rango de incertidumbre. La principal fuente de desajuste es la cuestión de la resolución. Debido a la huella del equipo geofísico y a la inversión de mínimos cuadrados amortiguados realizada en los datos TEM, tendemos a tener una transición más suave en los datos geofísicos que en la realidad. Dado que estamos comparando los resultados con datos de pozo, que tienden a reflejar estas transiciones bruscas, pueden ser una fuente de desajuste entre el \(CF_{rho}\) predicho y el \(CF_{log}\) real. Algunas capas delgadas de materiales resistivos en una capa conductora grande no quedarán atrapadas por la resolución del equipo pero aumentarán el valor residual. Sin embargo, la principal ventaja de la metodología14 es que la función traductora se adapta para contrarrestar estos efectos. La Figura 5 muestra los resultados de la función de traducción invertida y la incertidumbre promediada en 50.000 puntos TEM. La mayor parte de la incertidumbre se concentra en los puntos centrales que presentan una resistividad intermedia, ya que una pequeña incertidumbre sobre la resistividad en esta área afectará fuertemente la CF predicha. El rango de valores en el que se produce la transición entre arcilla y no arcilla con mayor probabilidad es entre 90 y 40 ohmios. Dichos valores son coherentes con la interpretación manual esperada ya que las Arcillas Cuaternarias de Aare tienen una resistividad entre 10 y 50 Ohmios. Los modelos se ciegan a la profundidad de investigación estándar obtenida de los modelos de resistividad y luego se pasan al paso de simulación geoestadística.

Los resultados presentados aquí se basan en 200 bucles de simulación de 10 realizaciones cada uno. Para cada bucle, se genera un nuevo conjunto TI y HD con una simulación GRF y se extrae un subconjunto aleatorio diferente para la validación cruzada. Los grupos se constituyeron utilizando un algoritmo K-Means de las coordenadas espaciales de los datos de CF. Como se mencionó anteriormente, el algoritmo de Deesse está controlado por tres parámetros principales: la fracción escaneada, el umbral y el número de vecinos. La fracción de exploración TI es común a todas las variables y se fijó en 5%. El patrón se comparará con el TI hasta que se alcance la fracción escaneada o hasta que se respete el umbral en todas las variables. El umbral se establece en 10% en tres nodos vecinos para la variable auxiliar y en 1% en 24 nodos para las variables principales. Si nunca se alcanza el valor umbral, se selecciona el mejor candidato y se marca el punto para volver a simularlo al final de la simulación. La ruta de escaneo en el TI es aleatoria. El resultado es un conjunto de 2000 simulaciones, basadas en 200 TI diferentes. El modelo 3D resultante tiene una dimensión de celda de 50 m por 50 my una resolución vertical de 2 m. La superficie recorrida es de 35 km2, para un tiempo de cálculo de 10 h en un cluster de CPU. Las secciones transversales verticales del modelo promediado se muestran en la Fig. 6. Hay una tendencia clara en el modelo entre el lado norte y sur del modelo. Tal variación era esperada debido a las variaciones en la forma del valle y está corroborada por los sondeos. La incertidumbre de los datos es mayor en las celdas más profundas del modelo, donde no hay perforaciones ni datos geofísicos que limiten las simulaciones. También denotan las zonas de transición entre áreas arcillosas y no arcillosas y reflejan la incertidumbre sobre la profundidad exacta de la transición. La forma de las estructuras subterráneas es consistente con los procesos conceptuales de deposición cuaternaria existentes. Sin embargo, la consistencia visual no es suficiente para confiar en un modelo.

Vista 3D de cortes del modelo de fracción de arcilla, con la incertidumbre asociada. La escala z está exagerada.

Histograma de los errores y el error normalizado para las dos variables principales: Resistividad (Rho) y fracción de arcilla (CF).

El error y las distribuciones de error normalizado [ver Ecs. (6) y (7)] calculados con validación cruzada se muestran en la Fig. 7. Durante la simulación, las dos variables principales son la resistividad y la fracción de arcilla. La estimación del error durante el paso de validación cruzada se realiza en ambos. El error se centra en cero para las dos variables, lo que sugiere que, en promedio, tendemos a predecir un valor correcto sin sesgo. Además, en el 44% y el 37% de los casos, respectivamente, logramos predecir un valor que estaba en el rango del 10% del valor real. En términos de error normalizado, la media del error normalizado para la resistividad es 0,98 (\(\sigma = 0,56\)). Esta distribución sugiere que las desviaciones estándar de las simulaciones son del mismo orden de magnitud que el error real. Por tanto, la incertidumbre está bien prevista. Por otro lado, la media del error normalizado del CF es 0,89 (\(\sigma = 0,5\)), lo que significa que tendemos a sobreestimar ligeramente el error en los datos simulados en comparación con los reales. Pero en general, la validación cruzada indica que, en promedio, se simula el valor correcto con una desviación estándar que refleja bien la posible incertidumbre. Un resultado sólido es que la resistividad se predice con precisión, lo que demuestra que este método puede ser de gran utilidad incluso para interpolar únicamente mapas de resistividad.

En un esfuerzo por homogeneizar y digitalizar todos los datos geológicos de las formaciones cuaternarias en Suiza, el proyecto GeoQuat realizó un prototipo y un estudio de demostración en el valle del Aare. En este contexto, se realizó un modelo geológico determinista del área mediante la interpretación manual de sondeos, de datos geofísicos y conocimiento geológico mediante el uso de secciones transversales geológicas39. Al comparar los dos modelos, debemos tener en cuenta que en el modelo estocástico se incluyen nuevos datos que no estaban disponibles en el momento de la construcción del modelo determinista: incluso si la base de datos de los pozos es la misma, los datos tTEM solo se adquirieron en 2020. Después de incorporar los pozos, correlacionaron manualmente unidades y facies y utilizaron la interpolación del vecino más cercano para extender el modelo a un volumen 3D completo39.

Comparación entre el modelo geológico existente y este modelo de estudio. Se describe la posición de la intersección de las secciones 3 y 4.

La Figura 8 muestra la comparación entre algunas secciones transversales en el modelo determinista y en el modelo generado con MPS. El modelo geológico determinista y los pozos se muestran utilizando sus componentes primarios USCS. La transparencia del modelo MPS refleja la incertidumbre. Cuanto mayor es la transparencia, mayor es la incertidumbre asociada. Ambos modelos concuerdan con los datos del pozo, al menos para estructuras de gran escala. Podemos ver que el modelo geológico generado por MPS muestra estructuras que son más realistas geológicamente hablando. Se pueden hacer algunos comentarios para las secciones transversales:

Las secciones 1 están de acuerdo. El modelo está bien limitado por dos pozos profundos que atraviesan todo el modelo. Incluso si la interpolación del vecino más cercano del modelo determinista muestra algunas formas de bloques poco realistas, las estructuras globales son similares entre los dos modelos.

Las secciones 2 son drásticamente diferentes entre el modelo determinista y el estocástico. La ausencia de perforaciones que limiten el modelo determinista provoca la propagación de una capa de grava desde perforaciones distantes, que probablemente no esté presente en este lugar. La resistividad invertida calculada a partir de los datos de tTEM es baja y probablemente esté asociada con cuerpos arcillosos gruesos. Esta sección muestra una mejora en los modelos, debido a los datos adicionales.

La sección 3 ilustra la importancia de la incertidumbre. No se han adquirido datos tTEM en el área y ambos modelos se basan en los mismos datos. Los metros superiores de los modelos están bien identificados como el cuerpo de grava superior en ambos modelos. Sin embargo, el uso del método propuesto reveló un posible cuerpo resistivo (arena o grava) más profundo. La incertidumbre es alta, como lo pone de manifiesto la gama cromática. No se puede establecer con seguridad la presencia de esta capa. Aunque es probable que sea una posibilidad que haya que tener en cuenta a la hora de evaluar la zona.

Las secciones 4 muestran cómo una estructura Cuaternaria local puede afectar drásticamente el modelo determinista. El pozo de 100 m de profundidad en el centro de la sección transversal indica una capa de grava que es anormalmente gruesa en comparación con los otros pozos cercanos en la sección. Este sobreespesamiento es extremadamente local y probablemente se deba a un movimiento de masa gravitacional hace unas decenas de miles de años, según comunicación personal con un geólogo del Cuaternario. Debido a esta perforación especial, el modelo geológico determinista sobreestimó el espesor del acuífero superior en un área muy amplia.

Finalmente, pudimos probar la calidad del modelo comparando la geología observada en varios pozos recién realizados que no fueron tenidos en cuenta en ninguno de los dos modelos porque fueron perforados recientemente. Los nuevos pozos poco profundos (máximo 10 m) concuerdan bien con ambos modelos en un 85% y 94% respectivamente para los modelos determinista y estocástico. Estos resultados no son sorprendentes ya que la capa superior de grava está más o menos presente en todo el dominio. La mayoría de los errores se deben a alguna sobreestimación o subestimación menor de la profundidad de la transición. Sin embargo, también se perforaron un pequeño número de pozos profundos. En la Fig. 9, se compara un nuevo registro de pozo con los dos modelos. El pozo se encuentra en el centro del pueblo de Wichtrach, donde no se pueden registrar datos tTEM. El modelo estocástico refleja bien la geología observada en el pozo. El modelo estocástico predice la presencia del acuífero de grava en la parte superior y altas resistividades como el modelo determinista hasta una profundidad de alrededor de 25 m. Pero luego el modelo estocástico predice la presencia de una capa de arcilla masiva hasta una profundidad de 80 m antes de indicar una gran incertidumbre a mayor profundidad. Esta predicción coincide bien con las observaciones geológicas, que también muestran una enorme capa de arcilla. Por el contrario, el modelo determinista predijo una enorme capa de grava justo debajo de una fina capa de arcilla, seguida de limos. Esto no es lo que se ha observado en el pozo. Ambos modelos están limitados por los mismos datos de pozos alrededor, pero mientras el modelo estocástico infiere la tendencia regional integrando todos los datos del área, el modelo determinista solo está limitado por puntos y secciones transversales cercanos.

Comparación entre un nuevo pozo en el pueblo de Wichtrach (no incluido en ninguno de los dos modelos), y la celda correspondiente en el modelo estocástico y el modelo geológico determinista.

En este estudio, presentamos una metodología que permite propagar la incertidumbre desde los datos geofísicos y los registros de pozo hasta el modelo final de la fracción de arcilla. Además, todos los datos de entrada se utilizan como datos concretos dentro de su rango de incertidumbre. Un punto de datos ruidoso en el campo dará lugar a una celda más incierta en el modelo. Esta característica es un punto clave para obtener una estimación robusta del error en el modelo final. Además, ser capaz de integrar perforaciones inciertas debido a descripciones deficientes o ausentes, posicionamiento inseguro o datos faltantes es una gran ventaja. El uso de MPS permite generar modelos que presentan patrones complejos que son imposibles de reproducir mediante métodos como el kriging. Este método también tiene la ventaja de depender de los propios datos para deducir patrones espaciales sin ningún conocimiento previo. Por supuesto, esto sólo es aplicable a regiones donde la cobertura de datos es suficiente. Sin embargo, en una región con menor cobertura de datos o para aumentar la calidad del modelo, también se puede agregar conocimiento previo mediante el uso de variables auxiliares. Después de la preparación de los datos, la generación de modelos a escala de valles de alta resolución se puede lograr en unos pocos días, en comparación con meses para los deterministas.

En términos de técnica geofísica, la aplicación de tTEM en el alto valle del Aare ha demostrado que la técnica es fiable y capaz de proporcionar información relevante sobre la estructura 3D del subsuelo. Esto confirma la conclusión anterior de Sandersen et al.46. Pero, además, sostenemos que una forma adecuada y eficiente de integrar estos tipos de grandes conjuntos de datos 3D con datos de pozo sólo puede lograrse mediante métodos automáticos. También esperamos que la metodología propuesta pueda ser muy útil para la integración de conjuntos densos de datos aéreos.

Sin embargo, la estimación de la litología a partir de la resistividad también tiene sus limitaciones. En primer lugar, si el área presenta variaciones significativas de salinidad, puede introducir no unicidad. La misma resistividad puede corresponder a material resistivo saturado o arcilla. Entonces sólo dependeríamos de las perforaciones para identificar la naturaleza del subsuelo. El método propuesto debería ampliarse para incluir una estimación 3D de la salinidad como variable auxiliar adicional. Esto puede ser posible utilizando información hidrogeológica adicional. Además, el método propuesto es eficiente para distinguir entre capas resistivas (arena o grava) y capas conductoras (arcilla), pero la diferencia de resistividad entre arena y grava no permite que el método las distinga claramente. Esa es la razón principal por la que decidimos centrarnos únicamente en la identificación de las capas de arcilla.

Incluso sin simular el valor de CF, nuestro enfoque también debe considerarse al interpolar modelos de resistividad 1D o 2D a 3D, especialmente cuando se trata de inversiones pronunciadas. Tradicionalmente, esto se hace mediante Kriging. Pero al ser kriging el mejor estimador lineal imparcial, conduce a una interpolación suave entre modelos nítidos e inevitablemente agrega artefactos en los modelos de resistividad 3D. En el contexto de los depósitos del Cuaternario, el rango del variograma a veces puede ser menor o equivalente a la distancia entre dos áreas de adquisición, lo que lleva a una interpolación errónea o incompleta. El método propuesto tiene la ventaja de depender directamente de los sondeos 1D en lugar de una cuadrícula realizada mediante kriging e interpolará entre los modelos con una nitidez coherente. Una metodología similar propuesta por Vilhemlmsen et al.15 se basa en un indicador kriging sobre el valor del grupo con un rango corto para llenar los espacios entre las líneas de adquisición en el TI, provocando un posible suavizamiento excesivo. La metodología presentada aquí tiene la ventaja de no depender de ninguna interpolación previa de los datos antes de la simulación.

Por supuesto, la inversión en sí misma ya es una interpretación de los datos y, además, la incertidumbre asociada con los modelos geofísicos invertidos a menudo es difícil de estimar con exactitud. Una mejora adicional podría ser el desarrollo de un algoritmo que trabajaría directamente en las mediciones geofísicas, sin ninguna inversión, en lugar de utilizar modelos invertidos. Al hacerlo, podríamos comparar directamente los errores entre los datos de campo reales y el modelo CF final. También simplificaría el flujo de trabajo al evitar una doble inversión, una para el modelo de resistividad en sí y otra para el algoritmo CF.

La comparación entre el modelo existente y el modelo automático (Fig. 8) muestra el riesgo de utilizar un enfoque determinista. En la sección transversal 2, el modelo determinista predice un cuerpo de grava gruesa a unos 40 m de profundidad. Este cuerpo fue colocado allí basándose únicamente en pozos que se encuentran a cientos de metros de distancia (cerca de la sección transversal 1), y después de incluir datos geofísicos ahora estamos seguros de que no existe tal cuerpo allí. Por supuesto, se trataba de una interpretación basada en los datos disponibles en ese momento, pero no había indicios de incertidumbre. En este contexto, tomar decisiones basadas en este tipo de modelos es arriesgado y debe evitarse.

Un número cada vez mayor de países está desarrollando bases de datos centralizadas para albergar mediciones relacionadas con la geología. En este contexto, el uso de un algoritmo de agregación de datos ágil y confiable es un enfoque prometedor. Ser capaz de deducir patrones espaciales a partir de datos sin o con poco conocimiento previo evita la infusión de estructuras a partir de elecciones arbitrarias, lo que generalmente se hace durante una interpretación manual de modelos geofísicos. La estandarización de los métodos de descripción a través de diferentes contratistas de perforación y los datos de acceso abierto son los puntos clave para que la aplicación de este método sea rápida y más adaptable a nuevas áreas de campo.

Una mejora del método podría ser pronosticar los valores de permeabilidad. Si se realizan mediciones de permeabilidad, ya sea mediante pruebas de empaque o de bombeo, se podría agregar una nueva variable a la simulación para generar directamente campos de permeabilidad. Otra posibilidad sería la definición de una función de densidad de probabilidad acoplada entre resistividad, fracción de arcilla y permeabilidad a partir de un conjunto de conocimientos previos o de campo. Entonces podríamos generar directamente campos paramétricos. Finalmente, como se mencionó anteriormente, pensamos que la metodología de interpolación introducida aquí podría aplicarse a otros modelos geofísicos 3D, para completar mapas parciales. A menudo, debido a limitaciones de campo (áreas inaccesibles, datos corruptos, adquisiciones diferentes,...) algunas áreas están menos cubiertas que otras. Utilizando MPS no deterministas, con pirámides gaussianas, se podrían generar modelos paramétricos completos y homogéneos con una cuantificación adecuada de la incertidumbre.

En este artículo, demostramos que un flujo de trabajo novedoso que combina el algoritmo de estimación de fracción de arcilla14 con nuestro algoritmo estadístico de puntos múltiples modificado es un método sólido para la generación automática de un modelo de fracción de arcilla 3D. Los modelos 3D resultantes pueden ser utilizados por las autoridades locales o los directores de proyectos para planificar mejor el desarrollo y el uso subterráneo de la región. Estos modelos se pueden utilizar, por ejemplo, para localizar zonas potencialmente altamente permeables para la explotación de aguas subterráneas o desarrollos geotérmicos, o para evaluar la posible presencia de materiales de construcción en el subsuelo. Dado que el método proporciona estimaciones de incertidumbre, también puede ayudar a diseñar la adquisición de más datos.

Nuestro método tiene la ventaja de estar basado en datos y no depender de la interpretación manual de las estructuras. El flujo de trabajo es automático e incluye: (1) la inversión de la función traductora en los modelos de resistividad y los datos del pozo, (2) el ajuste automatizado del variograma utilizado para el modelo GRF y (3) la generación de diferentes TI y HD para cada simulación MPS. El usuario sólo tiene que seleccionar algunos parámetros, como el número n de vecinos o el umbral t para las simulaciones MPS, pero se pueden utilizar valores predeterminados. Además, reproduce estructuras que no se pueden modelar mediante la interpolación clásica de dos puntos, como el kriging. Al comparar la geología observada a lo largo de los pozos recién realizados con nuestro modelo, este estudio siempre predice bien la tendencia general como se muestra, por ejemplo, en la Fig. 9. El flujo de trabajo se puede incorporar fácilmente con bases de datos públicas, lo que permite a las autoridades actualizar automáticamente el modelo 3D con regularidad. cuando estén disponibles nuevos datos de campo.

Finalmente, el artículo también demuestra la eficiencia del método electromagnético transitorio (TEM) remolcado. Cuando la mayoría de las publicaciones anteriores utilizaron datos TEM aéreos, integramos datos TEM remolcados. Esto nos permitió alcanzar una resolución en el espacio y la profundidad que es inalcanzable con cualquier otro método geofísico en este momento. Esta adquisición de datos debería considerarse cada vez más, incluso para problemas de mediana o pequeña escala. La integración fácil, rápida y económica de múltiples tipos de datos, a menudo existentes, sólo puede ser beneficiosa sea cual sea el propósito.

Los datos tTEM utilizados en este artículo están disponibles en el archivo de acceso abierto Zenodo en la siguiente URL: https://doi.org/10.5281/ZENODO.4269887. El conjunto de datos está documentado en Neven et al.26. El conjunto de datos del pozo está disponible en la Oficina Federal Suiza de Topografía (Swisstopo), pero se aplican restricciones a la disponibilidad de estos datos, que se utilizaron bajo licencia para el estudio actual y, por lo tanto, no están disponibles públicamente. Sin embargo, los datos están disponibles a través de Swissopo previa solicitud razonable.

Los códigos geoestadísticos utilizados en este estudio se implementan en el paquete Python Geone disponible gratuitamente en https://github.com/randlab/geone. Para utilizar todas las funciones de DeeSse, se puede obtener gratuitamente una licencia académica. El modelo CF estará disponible gratuitamente.

Ringrose, P. y Bentley, M. Diseño de modelos de yacimientos (Springer, 2015).

Google Académico

Pyrcz, MJ y Deutsch, CV Modelado de yacimientos geoestadísticos (Oxford University Press, 2014).

Google Académico

Jørgensen, F. et al. Un método para el modelado cognitivo de vóxeles geológicos 3D de datos AEM. Toro. Ing. Geol. Reinar. 72, 421–432. https://doi.org/10.1007/s10064-013-0487-2 (2013).

Artículo de Google Scholar

Jørgensen, F., Høyer, A.-S., Sandersen, PB, He, X. & Foged, N. Combinación de técnicas de modelado geológico 3D para abordar variaciones en geología, tipo de datos y densidad: un ejemplo del sur de Dinamarca. Computadora. Geociencias. 81, 53–63. https://doi.org/10.1016/j.cageo.2015.04.010 (2015).

Artículo de Google Scholar

Wellmann, JF, Varga, MDI, Murdie, RE y Gessner, K. Estimación de la incertidumbre para un modelo geológico del cinturón de arenisca verde, Australia Occidental: conocimientos de la inversión geológica y geofísica integrada en un marco de inferencia bayesiano. Geol. Soc. 453, 41–56. https://doi.org/10.1144/sp453.12 (2017).

Artículo de Google Scholar

Sófocleo, M. et al. Modelado numérico integrado para la gestión del agua en toda la cuenca: el caso de la cuenca de Rattlesnake Creek en el centro-sur de Kansas. J. hidrol. 214, 179–196. https://doi.org/10.1016/S0022-1694(98)00289-3 (1999).

ADS del artículo Google Scholar

Henriksen, HJ y cols. Metodología para la construcción, calibración y validación de un modelo hidrológico nacional para Dinamarca. J. hidrol. 280, 52–71. https://doi.org/10.1016/S0022-1694(03)00186-0 (2003).

ADS del artículo Google Scholar

Kollet, SJ y Maxwell, RM Modelado integrado del flujo de aguas superficiales y subterráneas: una condición de límite de flujo superficial libre en un modelo de flujo de aguas subterráneas paralelo. Adv. Recurso Acuático. 29, 945–958. https://doi.org/10.1016/j.advwatres.2005.08.006 (2006).

ADS del artículo Google Scholar

Lemieux, J.-M., Sudicky, EA, Peltier, WR y Tarasov, L. Dinámica de la recarga y filtración de aguas subterráneas sobre el paisaje canadiense durante la glaciación de Wisconsin. J. Geophys. Res. 113, F01011. https://doi.org/10.1029/2007JF000838 (2008).

ADS del artículo Google Scholar

Linde, N., Binley, A., Tryggvason, A., Pedersen, LB y Revil, A. Caracterización hidrogeofísica mejorada mediante la inversión conjunta de la resistencia eléctrica de los orificios transversales y los datos del tiempo de viaje del radar de penetración terrestre. Recurso Acuático. Res. 42, 1–10. https://doi.org/10.1029/2006WR005131 (2006).

Artículo de Google Scholar

Archie, G. El registro de resistividad eléctrica como ayuda para determinar algunas características del yacimiento. Trans. OBJETIVO 146, 54–62. https://doi.org/10.2118/942054-G (1942).

Artículo de Google Scholar

Knight, R., Gottschalk, I. y Dewar, N. Física de rocas a escala de campo para aplicaciones cercanas a la superficie. Encíclica. Geol. 1, 884–899. https://doi.org/10.1016/B978-0-12-409548-9.12514-X (2021).

Artículo de Google Scholar

Caballero, R. et al. Mapeo de sistemas acuíferos con electromagnetismo aéreo en el valle central de California. Agua subterránea 56, 893–908. https://doi.org/10.1111/gwat.12656 (2018).

Artículo CAS Google Scholar

Foged, N., Marker, PA, Christansen, AV, Bauer-Gottwein, P. y Jørgensen, F. Modelado 3D a gran escala mediante la integración de modelos de resistividad y datos de pozo mediante inversión. Hidrol. Sistema Tierra. Ciencia. 18, 4349–4362. https://doi.org/10.5194/hess-18-4349-2014 (2014).

ADS del artículo Google Scholar

Vilhelmsen, TN y cols. Combinando métodos de agrupamiento con MPS para estimar la incertidumbre estructural de modelos hidrológicos. Frente. Ciencia de la Tierra. 7, 181. https://doi.org/10.3389/feart.2019.00181 (2019).

ADS del artículo Google Scholar

Matheron, G. Principios de geoestadística. Economía. Geol. 58, 1246–1266. https://doi.org/10.2113/gsecongeo.58.8.1246 (1963).

Artículo CAS Google Scholar

De Marsily, G. et al. Afrontar la heterogeneidad espacial. Hidrogeol. J. 13, 161–183. https://doi.org/10.1007/s10040-004-0432-3 (2005).

ADS del artículo Google Scholar

Chiles, J.-P. & Delfiner, P. Geoestadística: modelado de la incertidumbre espacial vol. 497 (Wiley, 2009).

MATEMÁTICAS Google Scholar

Strebelle, S. Simulación condicional de estructuras geológicas complejas utilizando estadísticas de puntos múltiples. Matemáticas. Geol. 34, 1–21. https://doi.org/10.1023/A:1014009426274 (2002).

Artículo MathSciNet MATEMÁTICAS Google Scholar

Mariethoz, G. & Caers, J. Geoestadística de puntos múltiples: modelado estocástico con imágenes de entrenamiento (Wiley, 2014).

Reservar Google Académico

Strebelle, S., Payrazyan, K. & Caers, J. Modelado de un yacimiento de turbidita en aguas profundas condicionado a datos sísmicos utilizando geoestadísticas de puntos múltiples. SPE anual. Tecnología. Conf. Anexo.https://doi.org/10.2118/77425-MS (2002).

Artículo de Google Scholar

Pirot, G., Straubhaar, J. & Renard, P. Simulación de series temporales de modelos de elevación de ríos trenzados con estadísticas de múltiples puntos. Geomorfología 214, 148-156 (2014).

ADS del artículo Google Scholar

de Carvalho, PRM, da Costa, JFCL, Rasera, LG & Varella, LES Simulación de facies geoestadísticas con patrones geométricos de un yacimiento de petróleo. Estoco. Reinar. Res. Evaluación de riesgos. 31, 1805–1822. https://doi.org/10.1007/s00477-016-1243-5 (2017).

Artículo de Google Scholar

Dall'Alba, V. et al. Simulaciones estadísticas de puntos múltiples en 3D del acuífero del Plioceno continental del Rosellón utilizando DeeSse. Hidrol. Sistema Tierra. Ciencia. 24, 4997–5013. https://doi.org/10.5194/hess-24-4997-2020 (2020).

ADS del artículo Google Scholar

Neven, A., DallAlba, V., Juda, P., Straubhaar, J. & Renard, P. Estimación del volumen de hielo y topografía basal mediante métodos geoestadísticos y mediciones de radar de penetración terrestre: aplicación a los glaciares Tsanfleuron y Scex Rouge, Suiza Alpes. Criosfera 1, 5169–5186. https://doi.org/10.5194/tc-15-5169-2021 (2021).

ADS del artículo Google Scholar

Neven, A., Maurya, PK, Christiansen, AV y Renard, P. ttem20aar: Un conjunto de datos geofísicos de referencia para sedimentos fluvioglaciares no consolidados. Sistema Tierra. Ciencia. Datos 13, 2743–2752 (2021).

ADS del artículo Google Scholar

Christiansen, AV, Foged, N. y Auken, E. Un concepto para calcular el espesor de la arcilla acumulada a partir de registros litológicos de pozo y modelos de resistividad para la evaluación de la vulnerabilidad a los nitratos. J. Aplica. Geofís. 108, 69–77. https://doi.org/10.1016/j.jappgeo.2014.06.010 (2014).

ADS del artículo Google Scholar

Mariethoz, G., Renard, P. & Straubhaar, J. El método de muestreo directo para realizar simulaciones geoestadísticas de múltiples puntos. Recurso Acuático. Res. 46, 1-14. https://doi.org/10.1029/2008WR007621 (2010).

Artículo de Google Scholar

Straubhaar, J. DeeSse guía del usuario. Tecnología. Rep. (Centro de Hidrogeología y Geotermia (CHYN), Universidad de Neuchâtel: Neuchâtel, 2019).

Meerschman, E. y col. Una guía práctica para realizar simulaciones estadísticas de múltiples puntos con el algoritmo de muestreo directo. Computadora. Geociencias. 52, 307–324. https://doi.org/10.1016/j.cageo.2012.09.019 (2013).

ADS del artículo Google Scholar

Mariethoz, G., McCabe, MF y Renard, P. Reconstrucción espaciotemporal de brechas en campos multivariados utilizando el enfoque de muestreo directo. Recurso Acuático. Res. 48, 12115. https://doi.org/10.1029/2012WR012115 (2012).

Artículo de Google Scholar

Straubhaar, J., Renard, P. y Chugunova, T. Estadísticas de puntos múltiples utilizando imágenes de múltiples resoluciones. Estoco. Reinar. Res. Evaluación de riesgos. 1, 1–23. https://doi.org/10.1007/s00477-020-01770-8 (2020).

Artículo de Google Scholar

Oriani, F., Borghi, A., Straubhaar, J., Mariethoz, G. y Renard, P. Simulación de datos faltantes dentro de series temporales de caudal utilizando estadísticas de puntos múltiples. Reinar. Modelo. Software. 86, 264–276. https://doi.org/10.1016/j.envsoft.2016.10.002 (2016).

Artículo de Google Scholar

Dietrich, CR & Newsam, GN Un método rápido y exacto para simulaciones estocásticas gaussianas multidimensionales. Recurso Acuático. Res. 29, 2861–2869. https://doi.org/10.1029/93WR01070 (1993).

ADS del artículo Google Scholar

Branch, MA, Coleman, TF & Li, Y. Un método de gradiente conjugado, interior y subespacial para problemas de minimización con restricciones limitadas a gran escala. SIAM J. Ciencias. Computadora. 21, 1–23. https://doi.org/10.1137/S1064827595289108 (1999).

Artículo ADS MathSciNet MATH Google Scholar

Juda, P., Renard, P. y Straubhaar, J. Un marco para la validación cruzada de simulaciones geoestadísticas categóricas. Ciencia espacial terrestre. 7, 1152. https://doi.org/10.1029/2020ea001152 (2020).

Artículo de Google Scholar

Hartigan, JA & Wong, MA Algoritmo AS 136: Un algoritmo de agrupamiento de K-medias. Aplica. Estadística. 28, 100. https://doi.org/10.2307/2346830 (1979).

Artículo MATEMÁTICAS Google Scholar

Kellerhals, P., Haefeli, C. & Tröhler, B. Hidrogeología Aaretal, entre Thun y Berna. agua y Oficina de la Industria Energética del Cantón de Berna (WEA) (Estela de documentación geológica suiza, Hidrología y Geología del Estado, 1981).

Google Académico

Volken, S., Preisig, G. & Gaehwiler, M. GeoQuat: Desarrollo de un sistema para la gestión sostenible, modelado 3D y aplicación de datos de depósitos cuaternarios. Boletín suizo de geología aplicada. https://doi.org/10.5169/SEALS-658182 (2016).

Graf, HR y Burkhalter, R. Depósitos cuaternarios: concepto de clasificación y nomenclatura estratigráfica: un ejemplo del norte de Suiza. El suizo J. Geosci. 109, 137-147. https://doi.org/10.1007/s00015-016-0222-7 (2016).

Artículo CAS Google Scholar

Schlüchter, C. El registro cuaternario más completo del antepaís alpino suizo. Paleogeogr. Paleoclimatol. Paleeco. 72, 141-146. https://doi.org/10.1016/0031-0182(89)90138-7 (1989).

Artículo de Google Scholar

Casagrande, A. Clasificación e identificación de suelos. Trans. Soy. Soc. Ing. Civil. 113, 901–930. https://doi.org/10.1061/TACEAT.0006109 (1948).

Artículo de Google Scholar

Auken, E. y col. Una descripción general de un algoritmo inverso estable y directo altamente versátil para datos electromagnéticos y eléctricos aéreos, terrestres y de pozos. Explorar. Geofís. 46, 223–235. https://doi.org/10.1071/eg13097 (2015).

ADS del artículo Google Scholar

Christiansen, AV y Auken, E. Una medida global para la profundidad de la investigación. Geofísica 77, 171-177. https://doi.org/10.1190/geo2011-0393.1 (2012).

Artículo de Google Scholar

Alumbaugh, DL & Newman, GA Evaluación de imágenes para inversión electromagnética 2-D y 3-D. Geofísica 65, 1455-1467. https://doi.org/10.1190/1.1444834 (2000).

ADS del artículo Google Scholar

Sandersen, PBE y cols. Utilizando el método electromagnético transitorio remolcado (tTEM) para lograr detalles cercanos a la superficie sin precedentes en el mapeo geológico. Ing. Geol. 288, 106125. https://doi.org/10.1016/j.enggeo.2021.106125 (2021).

Artículo de Google Scholar

Descargar referencias

Esta investigación fue apoyada por la Fundación Nacional Suiza para la Ciencia a través del proyecto Phenix (Subvención No. 182600). Los autores desean agradecer a Nikolaj Foged por su experiencia en la aplicación del algoritmo ACT.

Centro de Hidrogeología y Geotermia, Universidad de Neuchâtel, Neuchâtel, Suiza

Alexis Neven y Philippe Renard

Departamento de Ciencias de la Tierra, Universidad de Aarhus, Aarhus C, Dinamarca

Anders Vest Christiansen

Departamento de Geociencias, Universidad de Oslo, Oslo, Noruega

Philippe Renard

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

También puedes buscar este autor en PubMed Google Scholar.

AN realizó el análisis de datos, diseñó e implementó la metodología geoestadística. AVC participó en la aplicación del algoritmo de fracción de arcilla a los datos y supervisó el trabajo. PR obtuvo financiación para el estudio. Supervisó el trabajo, participó en el diseño del flujo de trabajo y participó en la preparación de datos, redacción y edición del artículo. Todos los autores revisaron el manuscrito.

Correspondencia a Alexis Neven.

Los autores declaran no tener conflictos de intereses.

Springer Nature se mantiene neutral con respecto a reclamos jurisdiccionales en mapas publicados y afiliaciones institucionales.

Acceso Abierto Este artículo está bajo una Licencia Internacional Creative Commons Attribution 4.0, que permite el uso, compartir, adaptación, distribución y reproducción en cualquier medio o formato, siempre y cuando se dé el crédito apropiado a los autores originales y a la fuente. proporcione un enlace a la licencia Creative Commons e indique si se realizaron cambios. Las imágenes u otro material de terceros en este artículo están incluidos en la licencia Creative Commons del artículo, a menos que se indique lo contrario en una línea de crédito al material. Si el material no está incluido en la licencia Creative Commons del artículo y su uso previsto no está permitido por la normativa legal o excede el uso permitido, deberá obtener permiso directamente del titular de los derechos de autor. Para ver una copia de esta licencia, visite http://creativecommons.org/licenses/by/4.0/.

Reimpresiones y permisos

Neven, A., Christiansen, AV y Renard, P. Modelo estocástico automático de fracción de arcilla en 3D a partir de datos de sondeo y sondeo tTEM. Informe científico 12, 17112 (2022). https://doi.org/10.1038/s41598-022-21555-z

Descargar cita

Recibido: 03 de noviembre de 2021

Aceptado: 28 de septiembre de 2022

Publicado: 12 de octubre de 2022

DOI: https://doi.org/10.1038/s41598-022-21555-z

Cualquier persona con la que comparta el siguiente enlace podrá leer este contenido:

Lo sentimos, actualmente no hay un enlace para compartir disponible para este artículo.

Proporcionado por la iniciativa de intercambio de contenidos Springer Nature SharedIt

Al enviar un comentario, acepta cumplir con nuestros Términos y pautas de la comunidad. Si encuentra algo abusivo o que no cumple con nuestros términos o pautas, márquelo como inapropiado.