Moisès Jaén (1)

Enric Vázquez (1)

Jesús Carrera (1)

Adolfo Castro (2)

J.Mª. Salvany (1) (1)

Departamento de Ingeniería del terreno de la Escuela Técnica Superior de Caminos Canales y Puertos de Barcelona, ETSCCPB. Universidad Politécnica de Cataluña, UPC. (2)

Instituto de materiales y suelos. Universidad Nacional de San Juan, UNSJ.

The hydrogeological model presented in this paper is included in the characterization of the accident in the mine tailing reservoir at Aznalcóllar (Seville) from April 25,1998. The study area is located in a zone where the alluvium aquifer of the Guadiamar River overlays discordantly the Almonte-Marismas aquifer.

In that location the Almonte-Marismas aquifer have a transition from siltly to sandy sediments with gravel levels. Consequently, there is a continuity in the groundwater flow and transport of the contamination. This affects the quality of the water for irrigation. The process has three main steps: 1) model conceptualization, 2) discretization, and 3) building numerical model. Discretization requires building a tridimensional finite element mesh, which represents the complex geometry of the sedimentary system, with erosive surfaces and discordances among formations, as well as the heterogeneity in the materials. The numerical flow and transport model reflects the contaminant movement. The result of the model confirms the groundwater connection, but shows that contaminant will spend more than ten years to arrive to the end of the alluvial and to the pumping wells.

El presente modelo se engloba en los estudios de evaluación del impacto sobre las aguas subterráneas del vertido de lodos mineros debido a la rotura de la balsa de lodos de la Mina de Aznalcóllar (Sevilla), el 25 de Abril de 1998. El área de estudio se sitúa en una zona en la que el acuífero aluvial del Río Guadiamar se dispone discordante sobre el acuífero de la formación Almonte-Marismas, que pasa de sedimentos limosos a sedimentos arenosos con niveles de gravas. Esta situación revela una continuidad en el flujo subterráneo y, por consiguiente en el posible transporte de solutos, pudiendo afectar a la calidad de las aguas de extracción de las comunidades de regantes colindantes. Tras establecer el modelo conceptual del sistema que se pretende simular, se ha seguido un proceso de construcción de una malla de elementos finitos tridimensional, que represente la geometría compleja que alberga un sistema sedimentario donde se observan distintas superficies de erosión y discordancias entre formaciones, así como los diferentes materiales y sus relaciones laterales que forman el sistema de estudio. Una vez discretizado el dominio, se ha simulado el flujo de agua y el transporte de solutos, determinándose tanto las posibles vías de circulación de los contaminantes, como los tiempos de recorrido. Los resultados del modelo confirman la conexión, pero ponen de manifiesto que los contaminantes que lleguen al final del aluvial aun tardarán más de diez años en contaminar los pozos de bombeo.

Este modelo se enmarca en el estudio que ha sido elaborado por el Grupo de Hidrología Subterránea (Figura 1) dentro del proyecto de caracterización y seguimiento de la contaminación en los acuíferos afectados por el vertido de Aznalcóllar, en las proximidades del Parque Nacional de Doñana.

Los objetivos por los que se realiza un modelo hidrogeológico de detalle de esta zona son entre otros, validar el modelo conceptual del área de estudio y cuantificar el riesgo de contaminación de los acuíferos que están en contacto con el acuífero aluvial del Río Guadiamar, para evaluar el transporte de los contaminantes provenientes de la incorporación a las unidades hidrogeológicas de los solutos derivados del vertido de lodos tóxicos de la rotura de la presa de la mina de Aznalcóllar, en Abril de 1998.

La metodología consiste en modelar el flujo y transporte de contaminantes en las aguas subterráneas, mediante el código TRANSIN III (Galarza et al. 1996), discretizando el dominio mediante una malla de elementos finitos tridimensional, de manera que permita representar las heterogeneidades del medio

La zona que se pretende estudiar en detalle (Figura 2), presenta una serie de características estratigráficas que son de gran importancia a la hora de evaluar el grado de afectación del transporte de contaminantes en el sistema hidrogeológico general.

En el marco geológico donde se encuadra el presente modelo se produce un contacto entre varios niveles transmisivos que produce una continuidad en el flujo, pudiendo provocar una continuidad en el transporte de contaminantes provenientes del aluvial del Río Guadiamar, donde se extendieron los lodos. Esta conexión hace posible que lleguen los contaminantes a los pozos de riegos más cercanos que extraen aguas subterráneas de las formaciones acuíferas aquí tratadas.

La extensión (Figura 2) del modelo viene condicionada por dos factores; el dominio de interés de cara a definir el alcance de la contaminación, y el dominio físico y geológico.

El dominio de interés es la zona que comprende el contacto de la formación Marismas con la formación Almonte (Figuras 2,3 y 4). Para modelar, reproduciendo la realidad lo mejor posible, se zonificará el dominio según la estratigrafía y geometría de la zona, siempre que sea posible La fm. Almonte está constituida por materiales de edad Pliocena- Pleistocena, estando formada por dos litologías principales, unos limos basales, y un cuerpo en la parte superior arenoso con niveles de gravas de unos 60m de potencia. El buzamiento general es hacia el S. La fm. Marismas es de edad Cuaternaria formada por acumulaciones de arcillas con intercalaciones de gravas en disposición horizontal. Su potencia es variable, llegando en algún punto a los 80m. Además en la zona que nos ocupa se da un tránsito gradual entre los sedimentos aluviales del río Guadiamar con los de la formación Marismas. La Fm Marismas se dispone discordantemente sobre la fm.

Almonte que se manifiesta por una superficie de erosión (Figura 3) definida a partir de datos de sondeos y de los numerosos cortes realizados por J.Mª Salvany (Figura 6).

En esta superficie se producen intersecciones en profundidad de las capas de la fm. Almonte que han sido determinadas por métodos gráficos a partir de las direcciones y buzamientos de las capas y de la topografía de la superficie erosiva, resultando la Figura 4. En ésta se muestra la cartografía tanto en superficie topográfica (cuando no se encuentra la fm Marisma), como en superficie erosiva (cuando ésta se dispone discordante sobre la fm. Almonte) de las capas consideradas. Esta información está extraída de las columnas litológicas de los abundantes pozos y sondeos que existen en la zona (Figura 5).

En el modelo hidrogeológico se consideran las formaciones constituidas por una serie de capas de diferente transmisividad (Figura 7). Cada capa tendrá un comportamiento hidráulico distinto, debido a las características hidrogeológicas. Se producirán flujos de unas capas a otras, tanto más verticales cuanto mayor sea el contraste entre permeabilidades de capas en contacto. Se pretende representar este efecto adoptando una discretización tridimensional del problema. Los flujos más verticales se producirán en las capas tanto de arcillas de la fm. Marismas como de limos de la fm.

Almonte, que actúan como semiconfinantes de las capas adyacentes.

Como se refleja en la Figura 8, se produce una conexión de las dos formaciones a lo largo de la superficie erosiva. El flujo en las capas de arcillas y de limos, tendrá una componente ascendente, ya que el acuífero entra en carga a una presión mayor que la superficie del terreno.

3.1. Parámetros hidráulicos.

Se adoptarán los parámetros hidráulicos que se han obtenido con motivo de este estudio, de la reinterpretación de los ensayos de bombeo existentes en los pozos de la zona (ITGE 1998). Estos se asignan a las diferentes capas, como se muestra en la Figura 9.

La asignación de las condiciones de contorno se basa en del modelo conceptual de flujo que tiene en cuenta la situación de los niveles piezométricos abatidos en el entorno inmediato de la zona reproducida por este modelo. La asignación es la siguiente:

Condición A: Se asigna condiciones de nivel fijo conocido a los contornos en los que la fm. Almonte no esta recubierta por la fm. Marismas. El nivel piezométrico se asignará como un nivel que integra en vertical los niveles de cada capa que intersecte el contorno, a excepción de la capa de limos de base.

Condición B: Se asigna condición de nivel fijo a los contornos en los que la fm. Marismas se dispone discordante sobre la fm. Almonte, y que dado los límites de este modelo se corresponden con el  nivel de arenas S3.

Condición C: Se asigna caudal nulo o muy bajo a los contornos verticales de las capas de la fm. Almonte en las que, coincidiendo con los límites del modelo, cambian lateralmente a limos.

Condición D: Se asigna caudal nulo o muy bajo a los contornos horizontales inferiores de las capas de la fm. Almonte, correspondiendo con los limos de base.

Condición E: Se asigna caudal nulo o muy bajo a los contornos verticales de las capas de arcillas la fm. Marismas, dado que se considera que los flujos horizontales son despreciables frente a los restantes.

Condición F: Se asigna caudal nulo o muy bajo a los contornos verticales de las capas de gravas la fm. Marismas, dado que se considera que estas se acuñan llegando a quedar confinadas por las arcillas.

Esta asignación de condiciones se ilustra en la Figura 10.

El borde norte del modelo es la zona donde los depósitos cuaternarios de origen aluvial empiezan a tener una génesis de interacción con los sedimentos de marisma, que son de tipo arcilloso-limoso como corresponde a un ambiente de este tipo. Estos sedimentos de marisma dejan a las gravas, provenientes del aluvial del Guadiamar, confinadas en ellos. Se da también una situación morfológica de abertura, en extensión hacia el sur, de la superficie erosiva que viene encajando a las gravas GM sobre los limos de la fm. Almonte, coincidiendo con un cambio de litología de la formación Almonte de limos a arenas y gravas.

En este límite superior del modelo, pues, se va a considerar la entrada de la contaminación a través

de la capa de gravas GM, que se dispone discordante sobre limos. Esta capa es la que se considera que transporta los elementos contaminantes que se vertieron directamente sobre las terrazas aluviales que afloran a lo largo del curso del río Guadiamar aguas arriba del límite norte del modelo que nos ocupa. La zona física donde se simulará la entrada del contaminante se muestra en la Figura 11. Se simulará una entrada continua a lo largo del tiempo de simulación, ya que hay evidencias de metales pesados que permanecerán en los suelos con una elevada concentración que se irán lixiviando a lo largo del tiempo.

Un modelo numérico requiere discretizar el dominio en una serie de elementos finitos de manera que cubran todo el espacio sin solaparse. El método de elementos finitos resuelve el problema de flujo y transporte de aguas subterráneas en los vértices de los elementos (nudos) y dentro de cada elemento interpola linealmente. El error cometido será pequeño si el tamaño del elemento es pequeño en relación a la variabilidad de la función interpolada (niveles y concentraciones). Por otro lado tampoco conviene hacer elementos muy pequeños debido a que el tiempo de cálculo depende del número total de estos nudos.

Se utiliza un esquema de discretización mediante tetraedros y prismas triangulares. Se han tenido que construir dos mallas tridimensionales que se corresponden a las formaciones de Marismas y de Almonte, ya que el contacto entre ellas es discordante y no se ha encontrado en la revisión bibliográfica un generador de malla tridimensionales que contemple intersección entre capas como una opción explícita.

Una vez realizadas y con la topografía de la superficie de erosión como nexo común, estas dos mallas han sido acopladas nudo a nudo de manera que resulte una malla única. La generación de la malla se ha hecho con los programas 2DUMG (Bugeda, 1990), y 3D TRIDI (Vives, 1996) atendiendo a que la malla se adapte a los contornos y geometrías que reproducen la geología del dominio modelado y que se utilizarán para delimitar las zonas de parámetros.

La malla resultante está formada por elementos finitos tridimensionales en forma de prismas de base triangular o de tetraedros, como se ilustra en la Figura 12. La capa de gravas G1 se ha

tratado como elementos bidimensionales situados entre las capas S1 y S2. Esto se justifica porque debido a su alta transmisividad la componente principal del flujo está contenida en el plano del elemento. Este tratamiento permite ahorrar elementos en la discretización y por tanto reduce los tiempos de cálculo. La malla resultante consta de 6387 nudos y 11744 elementos tridimensionales.

Se da una situación que es común en este tipo de formaciones sedimentarias; las capas aquí tratadas tienen una longitud media en planta de 10km de, mientras que en vertical estas tienen un espesor medio de 10 a 20m. Este hecho hace que en la discretización resulten elementos que proyectados en planta tengan longitudes medias entre nudos de unos 200m, mientras que sus alturas medias sean de unos 10m. Esta relación base-altura de 1:20 daría problemas en la resolución numérica de las ecuaciones que gobiernan el sistema, sobre todo en relación a los elementos tetraédricos. Para evitar problemas numéricos, se ha optado por exagerar la escala vertical en un factor de 20. Este cambio de escala conduce a un replanteo de las ecuaciones de flujo y transporte, introduciendo el cambio de variable z' = 20*z. El efecto que produce en el planteamiento del problema es una modificación de los parámetros identificativos del funcionamiento hidráulico y de transporte, que lleva a trabajar con estas modificaciones. Con esto se evitan problemas numéricos, y los resultados se modifican adecuadamente para deshacer el escalado de coordenadas, para así obtener la solución correspondiente al problema original.

Se ha simulado el problema en situación de flujo estacionario y en régimen transitorio de transporte, con un tiempo total de simulación de 20 años. Los intervalos de calculo son de 1 día durante el primer mes y mensuales hasta el final.

Se han adoptado los parámetros de conductividad hidráulica y condiciones de contorno que se muestran en la Figura 13 para cada capa, con el fin de simular el flujo en régimen estacionario según el modelo conceptual anteriormente descrito:

La piezometría en planta, Figura 14, muestra un gradiente que disminuye hacia el sur, mostrando una inflexión hacia el contorno CH 73 que es el que tiene el nivel piezométrico mas bajo dado por la condición de nivel fijo impuesta en él. Se observa en planta un aporte de aguas desde la formación Marismas hacia las capas inferiores dado por la inflexión de las curvas piezométricas en los contactos de las dos formaciones. Este hecho también se pone de manifiesto en las piezometrías

en sección. Se observa el aporte hacia la capa S3, que es la que tiene continuidad, dado que las inferiores pasan a limos en el límite sur.

El balance de masas (Tabla 1) muestra el caudal de agua que entra o sale por cada contorno, viéndose que los contornos CH 51 y CH 73 concentran la mayoría del caudal de salida:

Tabla 1: Balance de masas en régimen estacionario. Se muestran los caudales de agua que entran (caudal positivo) o salen (negativo) por cada zona de nivel fijo.

Como resultado de la simulación en el periodo de 20 años considerado, se presentan los gráficos de la Figura 15, donde se observa el avance del penacho de contaminante. El camino preferente de avance en los primeros meses es a través de la capa Gm. Para un tiempo de 2 años de avance se observa que el penacho se introduce en la formación Almonte, a través de las arenas S1, situándose el frente de avance aproximadamente a 3km de la sección de inyección. En la situación final pasados los 20 años de simulación, se observa que el frente de avance a través de las gravas de la Marisma tiene una tendencia a frenar. Esto se corresponde con el modelo conceptual de flujo que contempla las gravas acuñándose en el límite sur de la malla, aplicándose en el modelo una condición de caudal nulo en este límite, haciendo que el flujo se dirija hacia las capas inferiores de la formación Almonte, tal como se observa en la sección de la piezometría en el estacionario. En cambio en sección vertical se observa que el frente de avance se divide en dos, a través de las arenas S1 empieza a tener relevancia un frente diferenciado, aunque este se encuentra retrasado respecto del superior, pudiendo ser esta vía la que conduzca a los contaminantes a continuar su transporte.

Los resultados obtenidos nos llevan a validar el modelo conceptual utilizado para esta hipótesis de flujo. Este modelo se podría utilizar como herramienta de control y gestión de los acuíferos aquí tratados, pudiendo ser capaces de determinar el alcance del vertido en un intervalo de tiempo determinado. Se está trabajando con otras hipótesis de flujo considerando los bombeos de riego que extraen agua de las formaciones aquí tratadas, así como condiciones de contorno mas precisas, extraídas del modelo de flujo regional realizado por la UPC para esta ocasión. Con estas nuevas hipótesis y a la espera de datos de concentraciones reales, se simulará en etapas posteriores un funcionamiento hidráulico que reproduzca lo mas fielmente posible los datos de niveles y concentraciones que se dispongan, y así poder determinar el alcance del vertido tóxico y diseñar medidas correctoras en medida de lo posible.

Este estudio ha sido posible gracias al patrocinio del ITGE (Instituto Tecnológico y GeoMinero de España) a través de su financiación y de aporte de datos. El autor quiere agradecer a todas las personas que han colaborado en la recopilación y tratamiento de información y que sin su apoyo e interés mostrado este trabajo no hubiera sido posible.

Bugeda G., 1990 "Utilización de técnicas de estimación de error y generación automática de mallas en procesos de optimización estructural ". Tesis Doctoral UPC.

Galarza, G. , Medina, A. , Carrera, J. – "TRANSIN III – Fortran code for solving the coupled non-linear flow and transport inverse problem. 1996. ETSICCP- UPC.

ITGE, 1998. Modelo regional de flujo subterráneo del sistema acuífero Almonte-Marismas y su entorno. Trabajo inédito.

Salvany, Josep Mª. 1999. Estudio sedimentológico de detalle del contacto aluvial del Guadiamar – Marismas. Trabajo inédito. Vives L., 1996. "TRIDI, manual del usuario". Informe interno UPC.