L. Vives (1) , A. Medina (1) y J. Carrera (1)

(1) Escuela Técnica Superior de Ingenieros de Caminos, Canales y Puertos de Barcelona, Universidad Politécnica de Cataluña. Jordi Girona 1-3. Campus Nord. 08034 Barcelona, España.

El estudio de los medios fracturados de baja permeabilidad tiene un elevado desarrollo en la actualidad por su importancia como posibles receptores de residuos radioactivos. La simulación del flujo de agua subterránea y de transporte de solutos en estas formaciones, es fundamental para los análisis de seguridad de estos almacenamientos. Existen diferentes tipos de modelos, siendo los Modelos Discretos de Fracturas (MDF) y los Modelos Continuos con fracturas principales (MCFP), los más empleados.

El presente trabajo tiene como uno de sus objetivos el describir brevemente los modelos matemáticos MDF y MCFP y enumerar las ventajas e inconvenientes de ambos. Posteriormente, se describe cómo se ha acoplado la parte de generación de un modelo MDF (FRACAS) a un modelo MCFP (TRANSIN), con el objetivo de aprovechar las ventajas y reducir las carencias de ambos tipos de modelo. Este trabajo se desarrolló dentro del marco del proyecto HIDROBAP en el que se propone el definir, elaborar y evaluar una metodología alternativa y complementaria para la caracterización estructural e hidrogeológica de medios geológicos fracturados, así como la posterior modelización del flujo y transporte de radionucleidos en dicho medio.

El articulo se completa con la simulación de un ensayo de interferencia realizado en el batolito de El Berrocal (Toledo, España), empleando con ambos modelos.

The studys of low permeability fractured media is growing fast, as they are thought as possible radioactive waste repositories. The simulation of both, groundwater flow and solute transport is basic for a security analysis of those repositories. There are different possible models, being the discrete fracture models (DFM) and the continuous model with embedded main fractures (CM) the most currently used.

In this work, we briefly describe those models, as well as their main advantages and disadvantages. Later, we show how a fracture network generated by a DFM model (FRACAS) can be included into a MCFP model (TRANSIN) trying use the best of each model type. This work has been developed as part of project HYDROBAP. The objective of this project is the definition of a suitable and alternative methodology to the structural and hydrogeological characterization of fractured media, as well as the posterior simulation of flow and solute transport of radionucleids in those media.

The article is completed with the simulation of a cross hole test done at El Berrocal (Toledo, Spain), using both types of models.

El presente trabajo es fruto de algunas tareas realizadas dentro del proyecto HIDROBAP (Elorza et al., 1999), financiado conjuntamente por la Empresa de Residuos Radiactivos (ENRESA) y el Consejo de Seguridad Nuclear (CSN). Este proyecto tiene como objetivo el desarrollar y evaluar una metodología de simulación hidrogeológica basada en un modelo de fracturas discretas destinadas a ser utilizado en el análisis de seguridad de almacenamientos geológicos profundos en rocas fracturadas. Esta investigación se lleva a cabo utilizando como lugar de estudio el batolito de El Berrocal (Toledo, España).

Este proyecto, en líneas generales, consiste en el desarrollo de las siguientes etapas: 1) Caracterización del medio geológico fracturado, 2) Simulación del medio fracturado, 3) Modelización del flujo subterráneo y transporte de radionucleidos y, 4) Verificación de la metodología desarrollada.

En la etapa de modelización del flujo subterráneo y transporte de radionucleidos, se aborda la modelización del flujo hidrogeológico mediante el empleo de diversos tipos de modelos: a) Modelos Discretos de Fracturas, b) Modelos Homogéneos Continuos con las fracturas principales y c) Modelos Heterogéneos Continuos con las fracturas principales. Esta etapa es la que da origen a este artículo, donde se describe la integración de los modelos discretos de fracturas con los modelos de tipo continuo, de forma que se aprovechen las características más importantes de cada tipo.

El artículo comienza realizando un breve repaso de los diferentes tipos de modelos para la simulación de medios fracturados, describiendo el modo en que se caracteriza este tipo de medios. Posteriormente, se describe el modelo discreto de fracturas FRACAS (Cacas, 1989) y el modelo continuo TRANSIN (Medina et al., 1996), para continuar mostrando en detalle el acoplamiento de ambos. Por último, se presenta un ejemplo real de un ensayo hidráulico de interferencia interpretado con ambos modelos.

Por completitud, se presentan brevemente las diferentes alternativas para la modelización de medios fracturados. En la actualidad hay cuatro formas de atacar el problema (Carrera, 1987):

• Medio continuo equivalente. En este caso, se considera que el medio (incluidas las fracturas) se puede aproximar de forma adecuada mediante un medio poroso con las ecuaciones clásicas de flujo y transporte. La principal dificultad en este tipo de aproximaciones estriba en la asignación de valores a los parámetros del modelo (transmisividad, almacenamiento, etc), porque éstos no responden a una realidad física, sino al modo en que se está aproximando la realidad.

• Modelos discretos de fracturas. En este tipo de modelos se considera que el flujo de agua a través de la matriz rocosa es despreciable comparado con el flujo que fluye por las fracturas. En este caso se representan únicamente las fracturas y se analiza el flujo en las mismas. La principal dificultad se encuentra en la definición adecuada de la geometría de la red de fracturas. El flujo en las mismas se suele analizar mediante la ley cúbica o alguna modificación de la misma, sin resolver explícitamente la ecuación de flujo, lo que simplifica los cálculos.

• Modelos mixtos. En estos modelos se intenta recoger los aspectos de los dos tipos de modelos comentados anteriormente. Por un lado se representan de forma explícita las fracturas dominantes y por otro lado, la matriz y las fracturas menores que hay entre ellas se representa como un medio poroso. La dificultad principal se encuantra en la asignación de valores a los parámetros del medio poroso (similar al problema que se tiene en el primer tipo de modelos considerados). A diferencia de los modelos discretos de fracturas, la definición de la geometría de las fracturas no es, en principio, tan problemática, porque generalmente se conocen con un cierto grado de precisión las fracturas dominantes.

• Modelos de red de canales. La hipótesis en la que se fundamentan es que el flujo tiene lugar preferentemente a través de una porción relativamente pequeña del plano de las fracturas. El flujo se suele calcular mediante el empleo de alguna expresión del tipo de la ley cúbica. La dificultad principal estriba en la definición de la geometría del problema. De todos los tipos, es el que menos se ha empleado.

En esta fase se analizan los datos disponibles sobre la fracturación (estructurales, emanométricos, petrológicos, etc.) y de flujo hidrogeológico, para:

• Extraer la relación, si existiese alguna, entre las características del flujo hidrogeológico y los parámetros estructurales o mineralógicos individuales de cada fractura o familia de fracturas.

• Calcular las propiedades de las distribuciones estadísticas geométricas de las familias de fracturas necesarias para construir un modelo de fracturas discretas (MFD) del batolito del Berrocal. Por tanto, en esta etapa, se realizan los siguientes tipos de caracterización del medio:

• Caracterización Geométrica: Se determina el número de familias de fracturas, distribución de la orientación de cada familia de fracturas, distribución del tamaño para cada familia, distribución espacial para cada familia e la intensidad de facturación para cada familia.

• Caracterización Hidráulica: Se clasifican hidráulicamente las diferentes familias de fracturas obtenidas en la etapa de caracterización. Además, se determinan las propiedades hidrogeológicas locales (fundamentalmente la conductividad hidráulica) del medio fracturado.

Estos modelos intentan describir el medio fracturado a partir de cada una de las fracturas conductivas. Esto lleva consigo tres pasos: el primero de generación de fracturas, implica la necesidad de una definición de las características de dichas fracturas (localización, orientación, buzamiento, densidad de fracturación, propiedades hidráulicas, etc.); el segundo de tipo geométrico, se corresponde con la construcción de la red de fracturación (seleccionando sólo las fracturas conductivas); y el tercero de tipo numérico, se refiere a la solución de las ecuaciones de flujo y transporte en dichas fracturas.

EL programa FRACAS que simula el flujo subterráneo con redes de fracturas discretas, desarrollado en la École Nationale des Mines de París, es el código seleccionado. Está estructurado en siete bloques muy precisos, donde cada uno de ellos tiene la siguiente función:

Bloque 1: Definición del dominio de generación. El dominio de generación de fracturas es un paralelepípedo, definido mediante las longitudes de sus tres aristas (en metros). El dominio de simulación de flujo debe estar dentro de este dominio de generación de fracturas.

Bloque 2: Generación de fracturas. En este bloque se generan las fracturas para cada una de las familias propuesta. Cada fractura estará simulada por un disco, por lo que cada una de ella queda definida por las coordenadas del centro, por el vector normal al plano del disco, por el radio del disco y la familia a la que pertenece (ver Figura 1). El programa primero genera, a partir de una semilla, una nube de puntos aleatorios que representarán los centros de los discos. Posteriormente, se generan los discos para cada familia. Para lo cual, las orientaciones de los discos se generan de acuerdo a una función de distribución para cada familia, donde se debe especificar la dirección media de la familia direcciones. Los radios de los discos se generan asignando una densidad de facturación (en centros por metro cùbico) y aplicando una distribución logarítmica de las longitudes de las trazas.

Bloque 3: Definición de condiciones de contorno. El programa dispone de dos tipos de condiciones de contorno: caudal nulo y nivel fijo. Ellas se pueden imponer en diferentes identidades geométricas: planos, esferas, cilindros y segmentos. El bloque finaliza seleccionando todas las fracturas que interceptan cada uno de estos limites y elimina aquellas que están fuera del dominio.

Bloque 4: Construcción de la matriz conductiva. Este bloque tiene dos objetivos, uno de seleccionar el conjunto de fracturas conductivas (eliminando el resto) y otro de la reconstrucción de la red de conexiones.

Bloque 5: Definición y calculo de las conductividades. Se da como dato de entrada para cada familia de fracturas el valor medio de los radios de los discos y la desviación tipo. El programa asume que la distribución de permeabilidades tiene una ley lognormal, por lo que genera números aleatorios de distribución normal N(0,1), que aplica a una media y desviación tipo especificados en la entrada de datos. A partir de éstas, calcula las permeabilidades para cada disco.

Bloque 6: Construcción y resolución del sistema de ecuaciones. El sistema de ecuaciones que resuelve este programa se obtiene de la aplicación del principio de conservación de masas y de la ley de Darcy. Este bloque realiza el ensamblaje de la matriz de flujo y impone las condiciones de nivel prescrito. Posteriormente, resuelve el sistema por medio de gradientes conjugados, obteniendo los niveles en cada uno de los centros de las fracturas conductoras.

Bloque 7: Cálculo de balance de masas. El programa calcula el balance de masas general, calculando parcialmente los caudales por cada límite impuesto y para cada uno de los nudos que componen dichos límites.

En Medina et. al. (1996) se describe el código TRANSIN que permite resolver y estimar los parámetros de las ecuaciones de flujo subterráneo y transporte de solutos. Aquí nos limitaremos a dar una breve descripción general del programa y pasaremos a describir sólo la forma en que este programa trabaja con los elementos unidimensionales, ya que a partir de ellos se formula la analogía con el modelo de fracturas discretas FRACAS.

El programa TRANSIN emplea el método de los elementos finitos para la solución de la ecuación.

Estos elementos pueden ser segmentos, triángulos, rectángulos, tetraedros y prismas de base triangular. Es decir, el programa puede trabajar en dominios espaciales unidimensionales, bidimensionales, cuasi-tridimensionales (simetría radial y multicapas) y tridimensionales. El régimen puede ser estacionario, transitorio o estacionario para flujo y transitorio de transporte. Las condiciones de contorno se pueden prescribir en los nudos el caudal, el nivel o expresar el caudal como función lineal del nivel en el acuífero. Tanto los niveles externos o prescrito como los caudales pueden variar de forma arbitraria en el espacio y en el tiempo. La estimación de los parámetros del modelo se hace a través de la minimización de una función objetivo derivada de la Teoría de Máxima Verosimilitud que permite incluir información previa sobre el valor de los parámetros. El programa permite estimar los parámetros de flujo y/o transporte.

El acoplamiento de estos dos códigos se realizó en el marco del proyecto HIDROBAP, por la necesidad de realizar calibraciones automáticas de parámetros de flujo en los modelos de redes de fracturas discretas. Por un lado se disponía del código FRACAS que genera redes de fracturas pero no calibra y por otro el código TRANSIN que calibra automáticamente pero no puede generar redes de fracturas. De lo que se trataba a través de este acoplamiento, es aprovechar la generación de las redes de fracturas de FRACAS y la calibración automática de TRANSIN.

El código FRACAS simula el caudal entre dos fracturas conductoras que se interceptan, como un tubo que conecta ambos centros de fracturas, con una conductancia que dependerá de las conductividades de ambas fracturas y de las distancias relativas de los centros de las fracturas a los puntos centrales de la intersección entre ellas. Por otro lado, los elementos unidimensionales en el programa TRANSIN simulan el flujo en un conducto entre dos nudos dados, es decir podemos asociar cada uno de los tubos que conectan las fracturas que se interceptan, como elementos unidimensionales. De esta forma, el dominio de estudio queda representado por un conjunto unidimensionales que simulan las conexiones entre las fracturas.

El acoplamiento entre los programas FRACAS y TRANSIN se mejoró al implementar: diferentes tipos de condiciones de contorno (por ejemplo, condición de goteo), el régimen transitorio (FRACAS sólo simula el flujo en régimen permanente), los puntos de observaciones para incluir medidas observadas necesarias en el proceso de calibración automática de parámetros, la zonificación de las familias de fracturas para poder calibrar por separado cada una de las zonas o familias de fracturas, y finalmente, se adapto el programa FRACAS para que pueda leer de ficheros exteriores de generación de fracturas realizada por otros programas.

En esta sección es presenta la calibración automática de parámetros de un ensayo de interferencia realizado dentro del proyecto El Berrocal (Rivas et. al. 1997), con un modelo continuo (TRANSIN) y uno de redes de fracturas discretas (FRACAS-TRANSIN).

En el ensayo de interferencia que se interpretó intervienen los sondeos verticales S14 (250.m) y S18 (233.1m), separados a 14.8 m de distancia. En ambos sondeos se han realizado varios ensayos hidráulicos (pulso, slug y de nivel constante) en diferentes tramo (Rogers et al.,1993; Guimerà et al., 1993). Los valores de conductividad hidráulica estan entre 10-6 y 10-13 m/s .

El ensayo hidráulico consistió en un bombeo intermitente en el sondeo S14, en el tramo S14.2 obturado de 3 metros de longitud a una profundidad de 93m. Los caudales de extracción fueron de 0,34 m 3 /d y de 0.55 m 3 /d. Se midieron presiones de agua en este tramo (S14.2) y en el inmediatamente inferior, S14.1, que se extiende desde los 97 m hasta el final del sondeo (250m).

En el sondeo vecino S18 se midieron presiones en tres intervalos a 66-75, 76-111 y 112-220 m de profundidad (S18.3, S18.2 y S18.1).

Desde un punto de vista conceptual, el ensayo se ha interpretado como de flujo convergente (Guimerà et al., 1996), que presenta un flujo preferente entre el tramo de bombeo (S14.2) y los dos tramos superiores del sondeo S18, asociados a una zona de fracturación. El radio de influencia del ensayo se adoptó de 400m. y se estudió un espesor de acuífero de 350 m.

En la Figura 2 se muestra un esquema del modelo numérico. Este consistió en una malla bidimensional axisimétrica que representaba el medio, con permeabilidad anísotropa. En el eje de simetría de situó el pozo de bombeo representado por elementos unidimensionales, así como los tramos de observación del sondeo S18. Se incluyó un conjunto de elementos 1-D distribuidos horizontalmente (para respetar la simetría del problema), desde el tramo de bombeo hasta el contorno de la malla para simular el flujo a través de la fractura. Los contornos superior, inferior y el eje de la malla son de caudal nulo y el extremo es de descenso nulo.

Los resultados de la calibración automática son buenos, desde un punto de vista de los ajustes, balance de masas y parámetros estimados. El excelente ajuste en el tramo de bombeo se logró a través de simular las pérdidas de carga como un efecto piel (Carrera et al., 1988). El tramo S14.1 tiene el peor ajuste, que asumimos previamente a la calibración, penalizando las medidas registradas por su escasa fiabilidad. Los parámetros estimados son coherentes con el modelo conceptual, es decir la fractura es la vía preferente de agua (1.97 10 ^-7m ² /s). La hipótesis de una formación con anisotropía en la permeabilidad no se ve reforzada por el ajuste, ya que tanto la conductividad vertical como la horizontal tienen un valor estimado alrededor de 5.09 10 ^-10 m/s siendo la componente vertical ligeramente superior (6.02 10 ^-10 m/s), como era de esperar por el análisis estructural.

El modelo conceptual empleado es similar al descrito al comienzo del apartado anterior. Es decir, se definió un dominio de estudio con el pozo de bombeo situado en el centro. Los tramos observados se simulan explícitamente y se supone que los descensos son casi despreciables (descensos nulos) en los contorno.

Figura 2: Estructura del modelo y malla de elementos finitos (283 nudos y 545 elementos). Hidrogramas de descensos medidos (puntos) y calculados (línea continua) en los tramos S14.2 (tramo de bombeo), S14.1, S18.3, S18.2 y S18.1. El esquema muestra la situación relativa de estos tramos respecto de la fractura. Extraido de Guimerà et al. 1996.

Tabla 1: Red de fracturación generada y conductiva. Características de cada familia de fracturas, con su densidad de fracturación (superficie total de fracturas dividido el volumen del dominio).

El dominio de estudio del modelo se definió como un paralelepípedo rectángulo de 400x400x350 m. El dominio de generación de la red de fracturación también es un paralelepípedo rectángulo pero de mayor tamaño (600x600x350m) y la red de fracturación incluye cinco familias: A, B, C, D, E y H (horizontal). La red generada es de 2517 fracturas (Tabla 1 y Figuras 3) y la red conductiva queda de 447 fracturas (Figuras 3), es decir solo 17.76% de la red inicial.

La malla de elementos finitos tiene 1503 nudos y 2102 elementos unidimensionales, lo que significa que las fracturas conductivas se interceptan entre ellas y con los puntos de observación 1051 veces, lo que equivale a una media de conexiones de 2.35.

La zonificación de las transmisividades y de los coeficientes de almacenamiento coinciden con las familias de fracturas. Las transmisividades de las fracturas se trataron como valor del parámetro la unidad y se varían los coeficientes de elementos de acuerdo con la transmisividad asignada en la generación de las familias. El almacenamiento tiene como valor del parámetro la compresibilidad del agua por un espesor representativo a cada familia (1 mm) y para los coeficientes de elemento se tiene en cuenta la densidad de fracturación inicial y el número de fracturas que intercepta cada una de ellas. Las condiciones de contorno impuestas son de caudal en el tramo S14.2 (bombeo) y de nivel fijo nulo (descenso nulo) en todas las fracturas que cortan las caras laterales del dominio.

Los resultados de la calibración de este modelo es aceptable y son comparables, desde un punto de vista de los ajuste de niveles (Figura 4), al del modelo continuo. Aunque, es conveniente dejar claro que no es objetivo de este trabajo el obtener ninguna conclusión con respecto a que tipo de modelo (MFD o continuo) es más realista, robusto o tiene un mejor funcionamiento, sino el de acoplar dichos modelos para aplicar el problema inverso.

Las conductividades hidráulicas de las familias de fracturas que se emplearon en la red conductora, están dentro del rango de variación obtenido con el modelo continuo por Guimerà et al. (1995). La razón que estos valores estén en la zona inferior del rango se atribuye a que estas conductividades son calculadas a partir del conjunto de fracturas conductivas de cada familia y no de la inicial (la red de fracturas conductivas contiene el 17.76% de las fracturas generadas).

Además, se debe tener presente que la asignación de transmisividades para cada fractura de la red conductoras es un proceso aleatorio que contiene varianzas elevadas.

En este trabajo se han presentan brevemente los diferentes tipos de modelos matemáticos para la simulación del flujo de agua subterránea y de transporte de solutos en formaciones fracturadas de baja permeabilidad. Posteriormente, se describe cómo se ha acoplado la parte de generación de campos de fracturas de un modelo Modelo Discreto de Fracturas (FRACAS) a un Modelo Continuo con fracturas principales (TRANSIN), con el objetivo de aprovechar las ventajas y reducir las carencias de ambos tipos de modelo.

El estudio se completa con la simulación, con ambos modelos, de un ensayo hidráulico de interferencia realizado en el batolito de El Berrocal (Toledo, España). Se analizan los valores de conductividad estimados y se sugiere que para un contraste más robusto se empleen las conductividades equivalentes del dominio.

El proyecto HIDROBAP ha contado con la financiación del CSN y de ENRESA. En la obtención y proceso de datos han participado numerosas personas que, sin su ilusión no hubiese sido posible la consecución de este trabajo. Los autores expresan su agradecimiento a todos ellos.

Cacas, M.C. (1989). `` Developpement d’un modele tridimensionnel stochastique discret pour la simulation de l’ecoulement et des transferts de masse et de chaleur en milieu fracture´´. These l’Ecole Nationale Supérieure des Mines de Paris.

Carrera, J. (1987). ``Avances científicos y tecnológico (hidrogeología de rocas fracturadas poco permeables)´´. XIII Hidrogeología y recursos hidráulicos. IV Simposio de Hidrogeología. España.

Carrera, J., J. Samper, L. Vives y U. Kuhlmann (1988). ``Automatic inverse methods for the analisys of pulse test; Aplication to four pulse test at tha leuggern borehole´´. Technical Report. NAGRA, NTB 88-34.

Elorza F., L. Vives y 20 autores más. "Higrogeología en medios de baja permeabilidad" Proyecto HIDROBAP, Informe Final, CSN-ENRESA.1999.

Guimerà, J., Carrera, J., Holmes, D.C., Rivas, P., Tallos, A. y Bajos C. (1993). "Preliminary analysis of the hydrogeology of El Berrocal experimental site" Memoires of IAH XXIV Congres, As, Oslo. Part1, 225-238pp.

Guimerà, J., Vives, L.. y Carrera, J. (1995). ``A discussion on scale effect on hydraulic conductivity at a granitic site (El Berrocal, Spain)´´. Geoph. Res. Letters 22(11)1449-1452pp.

Guimerà, J., L. Vives, M. Saaltink, P. Tume, J. Carrera y P. Meiers (1996). ``Numerical modelling of pumping test in a fractured low permeability medium´´. Topical Report 14, volume IV, El Berrocal Project. European Commision Contrac nº FI2W/CT91/0080.

Rivas, P., P. Hernán, J. Bruno, J. Carrera, P. Gomez, J. Guimerà, C. Marín y L. Perez del Villar (1997).``El Berrocal Project. Characterization and validation of natural radionuclide migration process under real conditions on the fissured granitic environment´´. Final Report. European Commision Contrac nº FI2W/CT91/0080, Nuclear Science and technology.

Rogers, S.F., D.C. Holmes, R.S. Ward, D.E. Bailey (1993). ``El Berrocal project. Progress report August 1992-January 1993´´. Fluid Process Group, BGS, NERC. Report Nº FP80FB/DOC4-AA.