Adolfo Castro Ledesma (*)

Enric Vázquez Suñé (**)

Jesús Carrera Ramírez (**)

Moisès Jaén Dupont (**)

Josep Mª Salvany (**)

(*) Instituto de Materiales y Suelos, Univ. Nacional de San Juan – UPC

(**)Dept. Enginyeria del Terreny. Universitat Politècnica de Catalunya (UPC)

Jordi Girona 1-3, Edificio D2. 08034 Barcelona, España

The fault of the mining waste disposal in Aznalcóllar (Andalucía, Southern Spain) in April 1998 produced the contamination of the affected soils. Among other tasks, the risk of impact on groundwater in the zone, motivated the making of a regional hydrogeological model which gathers together the available hydrogeological knowledge about the zone and acts also as a frame for two detail models focused towards the most affected areas. The regional conceptual model arises from the study of the previous information and recent studies. It considerates two-dimensional flow in porous satured media, in an aquifer with free behavior in the North area and confined behavior under marsh deposits (Marismas del Guadalquivir) in the Southwest area. The inverse problem is solved with automatic calibration of parameters using the program TRANSIN III (Galarza et al, 1996). The problem is solved in transient state for 28 years of simulation. A two dimensional finite element mesh is used, with an element size which is smaller at the interest zone, where the detailed models will be located. The role of the drainage system (rivers and brooks) is taken into account by carefully adjusting the mesh to those entities. The results validate the conceptual model in terms of hydrogeological parameters, mass balance and calculated piezometric head. The recharges are consistent too, with a trend of the model to calibrate low values in zones of intensive extraction. The results strongly suggests that the extractions are underestimated both in quantity, duration and spatial distribution.

La rotura de la balsa de lodos de la mina de Aznalcóllar (Andalucía, España) produjo un gran vertido con la consiguiente contaminación de los suelos en una extensa área. Ante el riesgo de posible afección a las aguas subterráneas, se ha emprendido, entre otras tareas, la elaboración de diversos modelos numéricos de flujo y transporte de contaminantes, dos de ellos de detalle en la zona del aluvial del Guadiamar y en el contacto de este aluvial con las marismas del Guadalquivir y uno del funcionamiento hidrogeológico regional que enmarque a los anteriormente mencionados. En este artículo se resumen las características principales y las conclusiones de este modelo regional. Se ha resuelto el problema inverso con calibración automática de parámetros, en régimen transitorio, con un tiempo de simulación de 28 años que abarca desde enero de 1970 hasta diciembre de 1997. Para ello se ha aplicado el código TRANSIN III (Galarza et al, 1996), que usa el método de los elementos finitos. La discretización espacial está condicionada por la geometría del modelo conceptual, estableciéndose una malla de elementos finitos triangulares. Los resultados validan en general el modelo conceptual en términos de parámetros hidrogeológicos y de balance de masas. El ajuste de los niveles es satisfactorio, sobre todo en las zonas de interés, por lo que este modelo regional podrá utilizarse como marco para los modelos de detalle. El ajuste de los niveles podría haberse logrado, ateniéndose a la información previa de extracciones, con parámetros de escasa coherencia. Se ha preferido priorizar la coherencia de los parámetros con la información previa, lo cual ha conducido a la necesidad de aumentar la cuantía de las extracciones.

Como consecuencia de la rotura, la noche del 24 al 25 de abril de 1998, de la balsa de lodos mineros de la mina de Aznalcóllar, unos cuatro o cinco millones de metros cúbicos de lodos piríticos y agua se extendieron sobre una superficie de unas 5000 ha a lo largo del valle del río Guadiamar, llegando a las proximidades del Parque Nacional de Doñana. Diversas instituciones, en parte coordinadas por el ITGE, están realizando tareas a fin de definir el alcance de la posible contaminación de los suelos y las aguas subterráneas como consecuencia del accidente. En este contexto, la Universidad Politécnica de Cataluña está realizando trabajos de síntesis de información hidrogeológica, hidroquímica y geoquímica. El objetivo final de estos trabajos es integrar dicha información con el fin de predecir la posible evolución futura de la contaminación.

Resulta evidente que para conseguir el objetivo señalado es preciso enmarcar la hidrogeología del valle del Guadiamar en su contexto regional, lo que constituye la motivación para la realización de un modelo de flujo regional.

El área que comprende el trabajo aquí descrito (Figura 1) ha sido objeto de diversos estudios, relevamientos, campañas de control y muestreo cuyo examen crítico e integración ha sido el punto de partida para la conceptualización del funcionamiento hidrogeológico del sistema.

En 1976 se realizó el primer modelo matemático del acuífero Almonte-Marismas, si bien el antecedente más inmediato lo constituye el modelo bidimensional (ITGE, 1992) y su actualización en 1994, utilizándose un esquema de diferencias finitas.

Por otro lado es necesario destacar los trabajos llevados a cabo tanto en la caracterización general del sistema Almonte Marismas, como muy especialmente en el área del Abalario (Custodio y Palancar, 1994; Salvany y Custodio, 1995) y en cuanto a caracterización de la recarga y funcionamiento hidrogeológico general (De Haro et al., 1998). También es de destacar la labor en lo tocante a hidrogeoquímica e hidrología isotópica (Manzano et al., 1991, Iglesias et al 1996 y 1998) y más recientemente en cuanto a la hidrogeología general del Abalario (Trick, 1998; Iglesias, 1999), entre otros.

A continuación se sintetizan las características más relevantes del modelo conceptual a partir del cual se va realizar la modelización numérica:

• Geológicamente se trata de un sistema sedimentario complejo, en parte deltáico, formado por paquetes permeables de arenas, limos y gravas. Hidráulicamete puede definirse como acuífero libre en el sector N para luego pasar a quedar confinado por la marisma en el sectorSE (Figura 1). El sustrato de baja conductividad hidráulica, y a efectos prácticos impermeable, lo constituyen las margas azules (mioceno). En la zona de El Abalario, el retrabajamiento marino y eólico dan lugar a una formación arenosa y de dunas que en parte recubren a las arenas y gravas de Almonte-Marismas dando lugar a la aparición de niveles piezométricos bien diferenciados.

• Los ríos actuan como elementos drenantes de los acuíferos condicionando fuertemente la piezometría (Figura 2), sobre todo en el sector N y en la zona del arroyo de La Rocina y de los ecotonos de La Vera y La Retuerta. Por su importancia estos deberán tenerse en cuenta en el modelo numérico.

• La recarga se produce en los sectores donde el acuífero se comporta como libre, vale decir en los limos y arenas del N y en el manto eólico superficial hacia el S.

Las descargas se producen a través de los ríos, arroyos y ecotonos mencionados, hacia el mar en la zona de El Abalario (aguas de infiltración relativamente reciente en el manto eólico), por extracciones y por flujos verticales a través de la marisma. En cuanto a la salida al mar a través de la estrecha franja que rodea a la desembocadura del Guadalquivir (flecha litoral), los estudios isotópicos y de salinidad parecen indicar que este flujo es muy lento y complejo.

En todo caso parecería que el flujo vertical ascendente a través de la marisma es muy relevante frente al anterior.

• La evolución de las extracciones a partir de la segunda mitad de la década del ’70 ha ido causando un sistemático descenso de los niveles con la aparición de conos de bombeo considerables. Esto ha motivado el desecamiento del ecotono de La Retuerta (Ecotono Norte en el contexto del modelo numérico) y la inversión del flujo las zonas de marisma expuestas a este efecto.

• El tipo de datos de que se dispone no permite un tratamiento tridimensional a nivel general ya que muchas de las medidas piezométricas corresponden a pozos ranurados en todos sus tramos permeables, lo cual no permite hacer consideraciones sobre los flujos verticales. Es por eso que se ha decidido hacer un modelo bidimensional, aunque multicapa, dada la estructura geológica y la disponibilidad de datos de los acuíferos en la zona del Abalario, donde existe constancia de diferencias de nivel del orden de uno a varios metros.

Se ha dividido la superficie del dominio en 1804 nudos y 3671 elementos triangulares. Estos elementos son más pequeños, del orden de 300 m. de lado, en la zona de contacto entre el río Guadiamar y las marismas, con el fin de acotar una zona con un margen de error moderado, ya que los resultados en esta zona se utilizarán como condiciones de contorno para el modelo de detalle 3D (1)

. El tamaño de los elementos es mayor hacia la zona W del modelo (área Moguer, Palos), llegándose a un tamaño de unos 3000 m. El tamaño medio de los elementos es de unos 1500 m.

En El Abalario, existe un manto eólico y formaciones arenosas del Plio-cuaternario con una intercalación de uno o más niveles limo-arcillosos que actúan como un acuitardo separando dos acuíferos (Trick,1998). Diferencias piezométricas significativas entre ambos acuíferos justifican un tratamiento especial. El esquema seguido para separar estos acuíferos ha sido convertir el modelo en bicapa, de manera que cada capa corresponda a un acuífero, y los elementos unidimensionales que las unen correspondan al acuitardo.

El tiempo inicial se ha tomado en 1970, dado que en esa época empieza a disponerse de una buena red de medidas piezométricas en todo el dominio y que la explotación en los acuíferos es muy incipiente por lo que puede asumirse que la superficie piezométrica corresponde prácticamente a la situación natural. La elevada cantidad de medidas piezométricas disponibles (más de 65.000 en unos 350 pozos) ha permitido realizar una discretización temporal con tiempos de observación mensuales y tiempos de cálculo semanales. Se inicia en Enero de 1970 y finaliza en Diciembre de 1997.

En el contexto del programa TRANSIN III, la variabilidad espacial de un determinado parámetro hidrogeológico Pi se introduce dividiendo el área del modelo en subdominios dentro de los cuales el parámetro es constante. El parámetro físico Pi, que en general se considera variable

siendo:

Pz , parámetro de zona

fe(x), coeficiente de elemento (o de nudo)

fz(t), función de tiempo

En cuanto a fe(x) es el coeficiente de nudo o elemento según que el parámetro esté asociado a uno u otro ente geométrico. Este coeficiente sirve para introducir la variabilidad espacial, supuesta conocida, dentro de la zona. Un ejemplo son las cotas de los ríos (introducidas como coeficientes de nudo) o la variación del espesor saturado (introducidas como coeficiente de elemento).

El parámetro de zona es el parámetro que calibra el programa, y está asociado en general a formaciones geológicas, usos del suelo, etc.

En cuanto a la función de tiempo, su uso más común está relacionado con variables que evolucionan en el tiempo, tales como la recarga, las extracciones, etc. El tratamiento de estas variables se realiza mediante una zonificación espacial como la descrita, y dentro de cada zona se calibra un factor multiplicador (parámetro de zona) Pz que afecta en esta caso a la función de tiempo.

En el modelo que aquí se describe se usaron cuatro funciones de tiempo para la recarga, calculadas mediante un balance de agua en el suelo para cuatro formaciones geológicas básicas: arenas, dunas, aluvial y limos. Estas cuatro funciones de tiempo se asignaron a las 25 zonas de recarga (Figura 3) condicionadas por la propia formación geológica y los diferentes usos del suelo.

Las particularidades de cada zona vienen dadas por el factor multiplicador a calibrar.

Las extracciones reciben un tratamiento similar. La zonificación (Figura 4) se hace de acuerdo a la información previa. En algunos casos se conoce el caudal medio extraído en toda un

área y en otros se conoce explícitamente la localización y el caudal de las extracciones. En el primer caso la extracción se asigna en forma areal, mientras que en el segundo caso la extracción se asigna en forma puntual.

A diferencia de tiempo, conocida distribución temporal.

El código parámetros basándose del valor del parámetro. Durante la fase de calibración a la información previa se le debe especificar la desviación tipo del parámetro, lo cual permite acotar los márgenes de variación permitidos al calibrar. Así, para sitios en donde determinado parámetro es bien conocido, la incertidumbre asociada al mismo será pequeña y también lo será la desviación. Lo contrario ocurrirá cuando exista escasez o ausencia de datos. Dado que la transmisividad es uno de los parámetros físicos de la naturaleza que varía en más órdenes de magnitud, y que su variabilidad espacial (aun en medios porosos "homogéneos") es grande, la desviación tipo será grande comparada con otros parámetros de menor variabilidad, como por ejemplo el coeficiente de almacenamiento. Para evaluar la calidad de los parámetros calibrados (en este caso T) se muestra la Figura 5, en la que se comparan los parámetros calibrados con la estimación previa, afectada de la desviación tipo.

A pesar de algunas diferencias, es innegable la coherencia general de los parámetros, sobre todo teniendo en cuenta que se trata de un modelo de magnitud regional, donde se están promediando características generales de las unidades. Con todo ello, hay que destacar que la distribución de zonas de transmisividad, su número y densidad obedecen de forma muy ajustada a los datos previos

6. NIVELES CALCULADOS

En la figura 6 se presenta el mapa de errores medios (diferencia entre nivel medido y calculado) del modelo. Este mapa es muy explícito por cuanto ayuda a valorar el grado de ajuste alcanzado por el modelo. La figura 7 es el mapa de desviaciones tipo. En general el ajuste es bueno en el área de interés afectada por los vertidos mineros (aluvial del Guadiamar, contacto de éste con la Marisma) y en la zona de El Abalario y regadíos.

7. RECARGAS Y EXTRACCIONES CALCULADAS

Para la recarga y las extracciones el valor de estimación previa del parámetro de zona (factor multiplicador que afecta a las funciones de tiempo) en principio se fija igual a la unidad. Al dejar que el modelo calibre las recargas, el análisis del balance de masas global (integrando los 28 años de simulación) denota una tendencia del modelo a disminuir el valor de dicho parámetro especialmente en zonas afectadas por intensas extracciones, e incluso en algunas zonas (zona de Pinares de Hinojos, al N de El Rocío) donde no hay evidencia cuantitativa acerca de extracciones y que por tanto no se han incluido en el modelo. Esto lleva a la conclusión de que puede existir una subestimación de los caudales de extracción.

Las estimaciones previas en cuanto a extracciones medias en el período de estudio son del orden de los 37 hm³ /año (ITGE, 1992). Este valor corresponde a las extracciones tomadas por el modelo matemático del ITGE, adoptada como estimación previa para la presente modelación.

El parámetro de zona (inicialmente igual a la unidad) ha sido, en muchos casos, aumentado considerablemente. Esto ha llevado a que en el balance global se obtenga un valor para las extracciones medio de 65 hm³ /año considerando el período efectivo de bombeo desde el año 1974. Las zonas del Rocío, San Ramón, Hatos, Pichardo, Pescantes, Almonte y los sectores I y II son las que han debido ser más aumentadas.

A lo dicho hay que agregar que las funciones de tiempo podrían ser realmente más largas de lo que indican los datos iniciales, ya que se consigue un mejor ajuste alargando algunas de ellas (Zona del Plan Regable, Almonte y Meneu Export). Esto indica que no sólo hay subestimación en volúmenes extraídos sino también en los tiempos de duración de la extracción.

En la figura 8 se comparan las extracciones totales (media anual) de la información previa y las extracciones calibradas a partir del modelo. Al obtenerse estos resultados se ha consultado a los organismos de cuenca la situación

actual de las extracciones, y a falta de la finalización del nuevo inventario de extracciones, los trabajos en curso confirman que hasta la fecha las extracciones habían sido subestimadas.

Se ha logrado ajustar muy bien tanto la variación espacial como la temporal de niveles, habiéndose conseguido este objetivo con parámetros coherentes con las características de hidrogeológicas de las formaciones y ajustados a la información previa, validando así el modelo conceptual existente.

La calibración y el ajuste de niveles mejora aumentando el caudal de algunas extracciones y si además extienden en el tiempo. Como conclusión de lo expuesto, se puede decir que la calibración del modelo mediante el problema inverso ha permitido la detección de errores en los datos de partida y la estimación de la magnitud de ese error.

La influencia de bombeos en esta área será decisiva a la hora de modelar flujo y transporte en detalle, los buenos ajustes de este modelo de regional permiten utilizarlo como encuadre a los diversos modelos de flujo y transporte de contaminantes que se están realizando.

A lo largo de los casi 30 años modelados se ha calculado una clara disminución del almacenamiento en los acuíferos debido básicamente a las extracciones de aguas subterráneas.

Custodio, E., Palancar Sánchez, M., 1994 "Las aguas subterráneas en Doñana". Revista de Obras Públicas, Madrid. 142 (3340): 31-53

J. M. De Haro, J. V. Giráldez, R. Ordóñez, E. Custodio, M. Iglesias, M. Manzano, J. J. López. 1998. "Variación temporal de la recarga en el Parque Natural de Doñana". EN PRENSA.

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

M. Iglesias, 1999 "Características hidrogeoquímicas del flujo del agua subterránea en El Abalario, Doñana, Huelva. Tesis Doctoral, UPC.

M. Iglesias, E. Custodio, J. V. Giráldez, M. Manzano, R. Ordóñe. 1996 "Caracterización química de la lluvia y estimación de la recarga en el área de El Abalario, Doñana, Huelva". 4ª SIAGA , Almeria, 1996. pp. 99-121.

M. Iglesias, T. Trick, E. Custodio, M. Manzano, J. V. Giráldez., 1998. "Estado actual de conocimiento de la recarga en el acuífero de Doñana". 1ª Asamblea Hispano-Portuguesa de Geodesia y Geofísica. Almeria. 25 pp.

ITGE, 1983. "Síntesis hidrogeológica de la cuenca del Guadalquivir". Colección Informes. Plan Nacional de Investigación de Aguas Subterráneas.

ITGE, 1992. Hidrogeología del Parque Nacional de Doñana y su entorno. Publ. Instituto

Tecnológico Geominero de España, Col. Informes aguas subterráneas y geotecnia, 64 pp, 2 mapas

Manzano, M., Custodio, E., Poncela, R. 1991. Contribución de la hidrogeoquímica al conocimiento de la hidrodinámica de los acuíferos en la área de Doñana. 1991. El agua en Andalucía – ITGE

Salvany, J. M.; Custodio, E., 1995. Características litoestratigráficas de los depósitos plio-cuaternarios del bajo Guadalquivir en el rea de Doñana: implicaciones hidrogeológicas. Rev. Soc. Geol. España, 8 (1-2), pp 21-31

T. Trick, 1998 "Impactos de las extracciones de agua subterránea en Doñana. Tesis Doctoral, UPC. pp. 6-15 a 11-6