Evaluación de la incertidumbre de las predicciones de los niveles piezométricos del acuífero de la fábrica de uranio de Andújar mediante el uso de filtros de Kalman de conjuntos

Título del Trabajo Fin de Máster: EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS Intensificación: RECURSOS HÍDRICOS Autor: GUTIÉRREZ ESPARZA, JULIO CÉSAR Director/es: DR. GÓMEZ HERNÁNDEZ, J. JAIME Fecha: ABRIL 2012 TITULO: EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL ALUMNO: JULIO CÉSAR GUTIÉRREZ ESPARZA DIRECTOR: J. JAIME GÓMEZ HERNÁNDEZ RESUMEN: Este trabajo se ha centrado en el área que es inmediata a la Uranio de Andújar (FUA), localiza, en el cual, tras la clausura de la fábrica lixiviados, que se desplazan lentamente desde los diques de estériles. Posterior al descubrimiento se han creado varios modelos matemáticos de flujo y de transporte de solutos, seleccionándose para este estudio el último de los modelos mencionados. Para evaluar la incertidumbre en las predicciones de los niveles piezométricos se hizo un análisis geoestadístico de las conductividades de la zona a estudiar, corroborando los patrones observados en trabajos anteriores; a partir del a obteniéndose mapas de la conductividad hidráulica de los diferentes materiales existentes en la zona. Obtenidos los mapas se desarrolló un código numérico que acopla la metodología de los filtros de Kalman de co para la simulación numérica de procesos de flujo de agua, transferencia de calor y transporte de solutos REactivos). Este código utiliza los valores medidos de piezometría para condicionar los campos de conductividad. El anál isis de los campos de piezometría resultantes nos permiten hacer una evaluación de la incertidumbre en las predicciones obtenidas con el código CORE. ABSTRACT: This work focuses on the aquifer unerlying the Andújar Uranium Mill (FUA). After the mill shutd the mill's tailings. Several mathematical groundwater flow and solute transport models have been created, choosing for this study the latest one. In order to evaluate the piezometric head's uncertainty analysis of the conductivity data in the studied area has been done, verifying the observed patterns in earlier assessments; from this analysis, different simulations were made, obtaining hydraulic conductivity maps for the study area. A fter obtaining the maps, a numerical code has been written that couples CORE (COde for flow, heat transfer and solute transport numerical simulation of REactive solutes) with the ensemble Kalman Filter EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS JULIO CÉSAR GUTIÉRREZ ESPARZA LUGAR DE REALIZACION J. JAIME GÓMEZ HERNÁNDEZ FECHA DE LECTURA Este trabajo se ha centrado en el área que es inmediata a la Uranio de Andújar (FUA), más concretamente en el acuífero que ahí se localiza, en el cual, tras la clausura de la fábrica se observó la presencia de lixiviados, que se desplazan lentamente desde los diques de estériles. Posterior al descubrimiento se han creado varios modelos matemáticos de flujo y de transporte de solutos, seleccionándose para este estudio el último modelos mencionados. Para evaluar la incertidumbre en las predicciones de los niveles piezométricos se hizo un análisis geoestadístico de las conductividades de la zona a estudiar, corroborando los patrones observados en trabajos anteriores; a partir del a nálisis, se realizaron diversas simulaciones obteniéndose mapas de la conductividad hidráulica de los diferentes materiales existentes en la zona. Obtenidos los mapas se desarrolló un código numérico que acopla la metodología de los filtros de Kalman de co njuntos al código CORE (CÓdigo para la simulación numérica de procesos de flujo de agua, transferencia de calor y transporte de solutos REactivos). Este código utiliza los valores medidos de piezometría para condicionar los campos de conductividad. El isis de los campos de piezometría resultantes nos permiten hacer una evaluación de la incertidumbre en las predicciones obtenidas con el código This work focuses on the aquifer unerlying the Andújar Uranium Mill (FUA). After the mill shutd own, this aquifer shows presence of contamination from the mill's tailings. Several mathematical groundwater flow and solute transport models have been created, choosing for this study the latest one. In order to evaluate the piezometric head's uncertainty analysis of the conductivity data in the studied area has been done, verifying the observed patterns in earlier assessments; from this analysis, different simulations were made, obtaining hydraulic conductivity maps for the study area. fter obtaining the maps, a numerical code has been written that couples CORE (COde for flow, heat transfer and solute transport numerical simulation of REactive solutes) with the ensemble Kalman Filter I EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS VALENCIA MAYO, 2012 Este trabajo se ha centrado en el área que es inmediata a la Fábrica de concretamente en el acuífero que ahí se se observó la presencia de lixiviados, que se desplazan lentamente desde los diques de estériles. Posterior al descubrimiento se han creado varios modelos matemáticos de flujo y de transporte de solutos, seleccionándose para este estudio el último Para evaluar la incertidumbre en las predicciones de los niveles piezométricos se hizo un análisis geoestadístico de las conductividades de la zona a estudiar, corroborando los patrones observados en trabajos nálisis, se realizaron diversas simulaciones obteniéndose mapas de la conductividad hidráulica de los diferentes Obtenidos los mapas se desarrolló un código numérico que acopla la njuntos al código CORE (CÓdigo para la simulación numérica de procesos de flujo de agua, transferencia de calor y transporte de solutos REactivos). Este código utiliza los valores medidos de piezometría para condicionar los campos de conductividad. El isis de los campos de piezometría resultantes nos permiten hacer una evaluación de la incertidumbre en las predicciones obtenidas con el código This work focuses on the aquifer unerlying the Andújar Uranium Mill (FUA). own, this aquifer shows presence of contamination from the mill's tailings. Several mathematical groundwater flow and solute transport models have been created, choosing for this study the latest one. In order to evaluate the piezometric head's uncertainty , a geostatistical analysis of the conductivity data in the studied area has been done, verifying the observed patterns in earlier assessments; from this analysis, different simulations were made, obtaining hydraulic conductivity maps for fter obtaining the maps, a numerical code has been written that couples CORE (COde for flow, heat transfer and solute transport numerical simulation of REactive solutes) with the ensemble Kalman Filter methodology. This program generates conductivity maps observed piezometric heads. The analysis of the ensemble of piezometric head maps allows analyzing their uncertainty RESUM: Aquest treball s'ha centrat en l'àrea que és immediata a la Fabrica d'Urani de Andújar (FUA), més concretament qual, després de la clausura de la fàbrica es va observar la presència de lixiviats, que es desplacen lentament des dels dics d'estèrils. Posterior al descobriment s'han creat diversos models matemàtics de flux i de de soluts, seleccionant Per a avaluar la incertesa en les prediccions dels nivells piezométrics es va fer una anàlisi geoestadístic de les conductivitats de la zona a estudiar, corroborant el es van realitzar diverses simulacions obtenint hidràulica dels diferents materials existents en la zona. Obtinguts els mapes es va desenvolupar un codi numèric q metodologia dels filtres de Kalman de conjunts al codi CORE (Codi per a la simulació numèrica de processos de flux d'aigua, transferència de calor i transport de soluts REactius). Aquest codi utilitza els valors mesurats de piezomteria per a c camps de piezometria resultants ens permeten fer una avaluació de la incertesa en les prediccions obtingudes amb el codi CORE. PALABRAS CLAVE: I ncertidumbre Conjuntos/ methodology. This program generates conductivity maps observed piezometric heads. The analysis of the ensemble of piezometric head maps allows analyzing their uncertainty Aquest treball s'ha centrat en l'àrea que és immediata a la Fabrica d'Urani de Andújar (FUA), més concretament en l'aqüífer que ací es localitza, en el qual, després de la clausura de la fàbrica es va observar la presència de lixiviats, que es desplacen lentament des dels dics d'estèrils. Posterior al descobriment s'han creat diversos models matemàtics de flux i de de soluts, seleccionant - se per a aquest estudi l'últim dels models esmentats. Per a avaluar la incertesa en les prediccions dels nivells piezométrics es va fer una anàlisi geoestadístic de les conductivitats de la zona a estudiar, corroborant el s patrons observats en treballs anteriors; a partir de l'anàlisi, es van realitzar diverses simulacions obtenint - se mapes de la conductivitat hidràulica dels diferents materials existents en la zona. Obtinguts els mapes es va desenvolupar un codi numèric q metodologia dels filtres de Kalman de conjunts al codi CORE (Codi per a la simulació numèrica de processos de flux d'aigua, transferència de calor i transport de soluts REactius). Aquest codi utilitza els valors mesurats de piezomteria per a c ondicionar els camps de conductivitat. L'anàlisi dels camps de piezometria resultants ens permeten fer una avaluació de la incertesa en les prediccions obtingudes amb el codi CORE. ncertidumbre / Analisis Geoestadistico/ Filtro de Kalman Conjuntos/ Conductividad Hidraulica/ Nivel Piezométrico II methodology. This program generates conductivity maps conditioned to the observed piezometric heads. The analysis of the ensemble of piezometric Aquest treball s'ha centrat en l'àrea que és immediata a la Fabrica d'Urani en l'aqüífer que ací es localitza, en el qual, després de la clausura de la fàbrica es va observar la presència de lixiviats, que es desplacen lentament des dels dics d'estèrils. Posterior al descobriment s'han creat diversos models matemàtics de flux i de transport se per a aquest estudi l'últim dels models esmentats. Per a avaluar la incertesa en les prediccions dels nivells piezométrics es va fer una anàlisi geoestadístic de les conductivitats de la zona a estudiar, s patrons observats en treballs anteriors; a partir de l'anàlisi, se mapes de la conductivitat hidràulica dels diferents materials existents en la zona. Obtinguts els mapes es va desenvolupar un codi numèric q ue acobla la metodologia dels filtres de Kalman de conjunts al codi CORE (Codi per a la simulació numèrica de processos de flux d'aigua, transferència de calor i transport de soluts REactius). Aquest codi utilitza els valors mesurats de ondicionar els camps de conductivitat. L'anàlisi dels camps de piezometria resultants ens permeten fer una avaluació de la incertesa en les prediccions obtingudes amb el codi CORE. Filtro de Kalman de Nivel Piezométrico / III Agradecimientos Quiero agradecer en primera instancia a la Conselleria de Educación, Formación y Empleo, en particular a Vicente José Ruiz Sánchez por su apoyo brindado en los diversos procesos que ha conllevado el ser becario del programa Santiago Grisolia. A Mi director de Tesis J. Jaime Gómez Hernández por el apoyo y por brindarme parte de su experiencia cuando he tenido dudas. También quiero agradecer a todos los miembros del grupo de Hidrogeología de la Universidad Politécnica de Valencia que, de una u otra manera me han aconsejado o ayudado a lo largo de mi estancia en España, nombrando de manera particular a Oscar, quien ha sido un buen compañero y amigo tanto en lo laboral como fuera de la oficina. Quiero agradecer a mis padres por darme una buena educación, a mi hermana por darme buenos consejos, quiero también agradecer a mi esposa, quien al igual que yo ha estudiado en esta universidad y que me da su punto de vista sobre el trabajo realizado, pero sobre todo, este documento se lo dedico a Natalia mi hija quien nacerá en este verano. V Índice general Índice general…………………………………………………………………………………………………………………………………V Índice de figuras………………………………………………………………………………………………………………….VII Índice de tablas………………………………………………………………………………………………………………….…XI 1 Introducción ........................................................................................................................ 1 1.1 Objetivos ...................................................................................................................... 5 1.2 Organización del documento ...................................................................................... 7 2 Estado del arte .................................................................................................................... 9 2.1 El uso de los filtros de Kalman ..................................................................................... 9 2.2 Aplicación de los filtros de Kalman de conjuntos en hidrogeología ......................... 11 2.3 Filtros de Kalman de conjuntos ................................................................................. 13 3 Modelo del FUA 2004 ....................................................................................................... 17 3.1 Gestión de datos del GIS............................................................................................ 19 4 Análisis geoestadístico ...................................................................................................... 23 4.1 Modelo de continuidad espacial ............................................................................... 23 5 Estimación ......................................................................................................................... 29 5.1 Krigeado ordinario ..................................................................................................... 29 6 Simulación ......................................................................................................................... 33 6.1 Simulación gaussiana secuencial ............................................................................... 33 7 CORE .................................................................................................................................. 41 7.1 Aspectos teóricos del código ..................................................................................... 41 7.1.1 CORE acuífero confinado .................................................................................... 41 VI 7.1.2 CORE acuífero libre ............................................................................................ 43 7.1.3 CORE zona no saturada ..................................................................................... 44 7.2 Solución numérica del flujo en CORE ........................................................................ 47 7.2.1 Solución del flujo en el acuífero ......................................................................... 47 7.2.2 Solución del flujo en la zona de transición ......................................................... 55 7.3 Datos de entrada al modelo CORE ................................................................................. 59 7.4 Simulación de alturas piezométricas modelo CORE ....................................................... 61 8 Filtro de Kalman de conjuntos .......................................................................................... 67 8.1 Filtro de Kalman de conjuntos aplicado al FUA ......................................................... 67 8.2 Descripción del programa.......................................................................................... 69 9 Resultados y conclusiones ..................................................................................................... 71 9.1 Resultados ...................................................................................................................... 71 9.2 Análisis de sensibilidad ................................................................................................... 79 9.2 Análisis de incertidumbre ............................................................................................... 83 9.3 Conclusiones ................................................................................................................... 87 10 Bibliografía........................................................................................................................... 89 VII Índice de figuras Figura 3.1 Malla de elementos triangulares creada a partir de la información de entrada utilizada por el código CORE……………………………………………………………………………………………...20 Figura 3.2 Cuadricula de 887040 elementos que envuelve el área de estudio………………..……20 Figura 3.3 Mapa de localización de los 154 datos de Conductividades.…………………..…………..21 Figura 4.1 Distribución espacial de los valores de conductividad hidráulica……………..….….….23 Figura 4.2 Histograma y curva de probabilidad acumulada…………………………………………………24 Figura 4.3 Variograma omnidireccional………………………………………………………………………….…...25 Figura 4.4 Variograma unidireccional en la dirección de los 30 grados azimut…………………….26 Figura 4.5 Variograma unidireccional en la dirección de los 120 grados azimut……………..…..26 Figura 4.6 Ejes de máxima y mínima continuidad de los datos de conductividad hidráulica..27 Figura 5.1 Mapa de conductividades obtenido mediante el Krigeado Ordinario………………….31 Figura 5.2 Mapa de la varianza en la conductividad hidráulica de la zona de estudio….………32 Figura 6.1 Mapa de conductividades hidráulicas (m/d) resultado de la simulación 22………..38 Figura 6.2 Mapa de conductividades hidráulicas (m/d) resultado de la simulación 35…..……38 Figura 6.3 Mapa de conductividades hidráulicas (m/d) resultado de la simulación 43………..39 Figura 6.4 Mapa de los valores esperados de las conductividades hidráulicas (m/d)………….39 Figuras 6.5 Varianza Condicional en del Área de estudio ……………………………………….……………40 Figura 7.41.1 Mapa con 12254 zonas de conductividades a partir de la simulación 11….……62 Figura 7.41.2 Mapa con 12254 zonas de conductividades a partir de la simulación 61…..…..63 Figura 7.41.3 Mapa de valores esperados de las 12254 zonas de conductividades a partir de el conjunto de simulaciones………………………………………………………………………………………………..63 VIII Figura 7.41.4 Mapa de varianza condicional de las 12254 zonas de conductividades a partir de el conjunto de simulaciones……………………………………………………………………………………………64 Figura 7.4.1.5 Alturas Piezométricas observadas en el PC1 comparadas con los valores que reproduce el CORE para el periodo 2003-2020………………..………………………………………………….64 Figura 7.4.1.6 Alturas Piezométricas observadas en el PC5 comparadas con los valores que reproduce el CORE para el periodo 2003-2020……………………………………………………………….…..65 Figura 7.4.1.7 Alturas Piezométricas observadas en el pozo 611 comparadas con los valores que reproduce el CORE para el periodo 2003-2020……………………………………………………..……..65 Figura 7.4.1.8 Alturas Piezométricas observadas en el pozo 472 comparadas con los valores que reproduce el CORE para el periodo 2003-2020…………………………………………………………….66 Figura 8.2.1 Diagrama de flujo del programa de Filtros de Kalman de Conjuntos añadido al código CORE………………………………………………………………………………………………………….……….…..70 Figura 9.1.1 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el pozo 611 realización 2 …………………………………………………………………………………….………..71 Figura 9.1.2 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el pozo 505 realización 2 . ………………………………………….……………………………………….………..72 Figura 9.1.3 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC6 realización 2 . ……………………………………………….….……………………………………….………..72 Figura 9.1.4 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC5 realización 2 . ………………………………………….……………………………………….…………….…..73 Figura 9.1.5 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC3 realización 2 . ………………………………………….……………………………………….…………….…..73 Figura 9.1.6 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC2 realización 2 . ………………………………………….……………………………………….…………….…..74 Figura 9.1.7 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC1 realización 2 .………………………………………….……………………………………….…………….…..74 EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 3 observadas, el valor máximo de la concentración y el valor de estabilización definido por el año en el que alcanza el valor asintótico de la concentración (Enresa 2004). Se obtuvieron varias hipótesis a partir de los resultados obtenidos, entre las que se mencionan que existe la posibilidad de que el uranio quede retenido o precipite debido a procesos geoquímicos en el acuífero, la posible salida de uranio a través del río debido a vías no consideradas en el modelo, o la posibilidad de que el uranio quede atrapado en la zona no saturada del aluvial debido a procesos hidrodinámicos. También se descartaron algunas hipótesis debido a que los resultados obtenidos al tomar en cuenta ciertos factores no aportaba una base clara de las diferencias tan marcadas entre las mediciones y el modelo. Entre estas hipótesis se encuentra por ejemplo la precipitación del uranio, la cual no guardaba una relación clara con las concentraciones de uranio. Los resultados mostraron que existía una fuerte relación entre los procesos de dilución causados hidrológicamente, y los valores medidos en los pozos de observación, por lo que se ajustó el modelo para reproducir este efecto. El hecho de que no existiese una razón clara para que se reprodujesen de manera congruente los valores llevó a la modificación del término fuente, empleando algunas hipótesis hidrogeoquímicas, de las cuales solamente se trataron a grandes rasgos algunas posibilidades debido a la gran exigencia de cálculo y la complejidad que se tenía. En cuanto a la heterogeneidad espacial de la zona, el estudio introdujo cambios ligeros en la estructura del modelo de flujo, para representar adecuadamente los datos de las zonas próximas a la instalación; también se recalibraron los parámetros del modelo de flujo y transporte siendo similares a los obtenidos en el modelo de 1994. Posterior a este reporte, en 2011, se realizó un análisis de incertidumbre del modelo FUA04 variando la conductividad hidráulica que se introduce en los ficheros de entrada del modelo. Para ello se llevo a cabo un análisis estocástico por medio del uso de técnicas de simulación de Monte Carlo (Gómez, 2011). El análisis requirió generar un modelo de continuidad de las variables aleatorias mediante una función multigaussiana a fin de crear múltiples realizaciones del campo de EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 4 conductividades que estuviesen espacialmente correlacionadas, posteriormente el modelo del FUA04 era dotado del nuevo campo heterogéneo de conductividades. Los resultados obtenidos por este análisis muestran que, en promedio, los niveles piezométricos no se ven alterados en gran medida al introducir heterogeneidades a pequeña escala, no obstante existe una gran discrepancia entre las observaciones y el modelo tanto al tomar en cuenta la heterogeneidad como cuando no es tomada en cuenta. El análisis señala que aun teniendo en cuenta la heterogeneidad del medio a pequeña escala, no se logra explicar en su totalidad la gran variación de concentraciones que se tiene en el medio, en parte por las pocas realizaciones que se llevaron a cabo, las cuales no logran ajustar adecuadamente las concentraciones. Es a partir de estos resultados que se decide emplear los filtros de Kalman de conjuntos para, intentar reducir esa discrepancia entre los valores observados y obtenidos, ajustando primeramente la piezometría para después utilizar las herramientas creadas a fin de acoplarlas al problema de transporte en un futuro estudio. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 5 1.1 Objetivos El presente estudio plantea observar la capacidad que tienen los filtros de Kalman de conjuntos para asimilar información existente a fin de ajustar las predicciones del modelo del FUA. Simultáneamente preparar el modelo que utiliza el código CORE para reproducir la piezometría y posteriormente en trabajos subsecuentes al alcance de este estudio, intentar reproducir las concentraciones de contaminantes observadas en los distintos puntos de la red de vigilancia. En resumen los objetivos del trabajo son: • Aplicación de los filtros de Kalman de Conjuntos a un caso real complejo. • Generar una serie de campos heterogéneos de conductividades que sean representativos de las zonas de materiales utilizadas en el anterior modelo del FUA, a fin de reproducir la heterogeneidad del medio a pequeña escala. • Modificar el CORE para poder utilizar un mayor número de zonas de material y así tomar en cuenta la heterogeneidad del terreno a pequeña escala, que, aunque para resolver el problema de flujo no es tan significativo, sí lo es para resolver el problema de transporte de un contaminante reactivo. • Hacer uso de la metodología de los filtros de Kalman de conjuntos para ajustar las conductividades hidráulicas utilizando conjuntos de campos de conductividades que, una vez ajustados por el filtro reproduzcan satisfactoriamente la piezometría de la zona. • Finalmente obtenidos los campos ajustados medir la incertidumbre que introduce el uso de diferentes conductividades hidráulicas al modelo. Entre los objetivos secundarios de este trabajo están:  Crear subrutinas para mayor facilidad en el manejo de los datos.  Crear o modificar según sea el caso los archivos de entrada para cada programa utilizado.  Generar los formatos necesarios para introducir los resultados y observarlos visualmente en el programa SGEMS. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 7 1.2 Organización del documento Este trabajo trata de tomar en cuenta la heterogeneidad a pequeña escala para poder reproducir adecuadamente los valores observados en campo en el área circundante a la FUA. Para esto se han empleado filtros de Kalman de conjunto, una herramienta que ha probado ser robusta en el campo de la ingeniería del petróleo. Por tanto se realiza una revisión del estado del arte de los filtros de Kalman de conjuntos, sus aplicaciones en diversos campos, los estudios realizados en hidrogeología más actuales que usan esta metodología, así como un ejemplo ilustrativo que intenta describir su funcionamiento de manera sencilla al lector de este documento, todo esto se verá más a fondo en el capítulo 2. El capitulo 3 nos habla acerca del modelo del FUA, la extracción de datos realizada y sus diversas modificaciones. El capitulo 4 hace hincapié al análisis geoestadístico realizado. Para ello se hizo uso de un modelo de continuidad espacial a fin de poder representar de manera adecuada la variación de valores que tienen el área de estudio en cuanto a conductividades se refiere. Luego de haber obtenido el modelo de continuidad espacial, este es usado en el capítulo 5 para realizar una estimación del campo de conductividades hidráulicas mediante un krigeado ordinario. En el mismo capítulo se describe como se realiza el krigeado, así como los resultados del mismo. El capitulo 6 muestra cómo, a partir del modelo de continuidad, se pueden recrear diversas realizaciones de campos equiprobables, con los cuales poder realizar el análisis de sensibilidad del modelo posteriormente. En este capítulo también se realiza un desarrollo del método de simulación que ha sido realizado con el programa SGEMS (Stanford Geolostatistical Modeling Software); por último se observan algunas de las simulaciones creadas, asi como los valores esperados del conjunto de simulaciones. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 8 El capitulo 7 hace una revisión al CORE, código que resuelve las ecuaciones de flujo, transporte de solutos y transferencia de calor en contornos irregulares. Se detallan los aspectos teóricos del CORE, la resolución de las ecuaciones de flujo en acuífero confinado, acuífero libre, así como en la zona con saturación variable. También en el capítulo 7 se describe la metodología que utiliza CORE para aproximar la solución real por medio de la solución numérica del flujo en CORE. Esta se describe para acuífero confinado y libre, y además la metodología iterativa que utiliza CORE para resolver el problema numérico en la zona de transición. Finalmente se describen los datos de entrada al modelo, así como la respuesta del modelo antes de introducir el uso de los filtros de Kalman de conjuntos. El capitulo 8 nos habla de los filtros de Kalman de conjuntos, como fueron aplicados al FUA, así como el programa que fue utilizado para ajustar las conductividades por medio de los archivos de salida provenientes de CORE. El capitulo 9 nos habla de los resultados y conclusiones que se obtuvieron al utilizar los filtros de Kalman de conjuntos. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 9 2 Estado del arte 2.1 El uso de los filtros de Kalman El filtro de Kalman fue introducido en 1960 por Rudolph Kalman quien en su publicación (Kalman, 1960) hizo una gran contribución a la teoría de control estocástica, hoy conocida como el filtro de Kalman. Es un método bayesiano secuencial de asimilación de datos, que en cada intervalo de asimilación combina una estimación previa con la medida observada para dar un estimado posterior, condicionado a que toda la información disponible en ese momento sea utilizada. Ambos el valor estimado previo y el posterior, como también la medición, contienen cierta incertidumbre. La solución posterior es generalmente considerada la solución óptima para el caso de procesos lineales con variables Gaussianas (Thulin, Skaug, Aanonsen, Nævdal, 2012). El filtro fue originalmente diseñado para solucionar el problema lineal clásico de la estimación de mínimos cuadrados en el procesamiento de señales y la teoría de control. Posteriormente el filtro de Kalman extendido (EKF) fue introducido para estimar el estado de un sistema en un proceso no linear. Básicamente el filtro de Kalman extendido lineariza la función para acercarse a la media estimada actual y calcula la matriz Jacobiana del estado de transición y de la función de observación en cada intervalo de cálculo. Cabe enfatizar que aun las aplicaciones tempranas del filtro de Kalman extendido introdujeron la estimación de parámetros del modelo desconocidos (Cox,1964;Kopp y Orford, 1963) añadiéndolos a la estimación del estado de la variable dinámica. Esto es conocido como la estimación dual. Sin embargo el EKF no es apropiado para modelos muy grandes o con no linearidades importantes. El filtro de Kalman de conjunto (EnKF) fue introducido como una alternativa por Evensen (1994). El EnKF usa un conjunto de realizaciones para representar los estadísticos del estimador actual, y todas las realizaciones son propagadas hacia adelante y analizadas de acuerdo a las ecuaciones del filtro de Kalman. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 10 Este método ha probado tener aproximaciones prometedoras, además de ser eficiente y robusto aun cuando se utiliza en modelos no lineales muy grandes y se ha vuelto muy popular entre aplicaciones como los modelos oceánicos (Evensen 2006), modelos de reservorios (Haugen et al., 2008; Skjervheim et al., 2007; Nævdal et al. 2005) y modelos atmosféricos (Sun et al. 2009, Kepert, Sun y Steinle, 2003). EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 11 2.2 Aplicación de los filtros de Kalman de conjuntos en hidrogeología Como se comentó, el uso de los filtros de Kalman de conjuntos ha ido en aumento en las ciencias atmosféricas, en el estudio de la interacción suelo-atmosfera, en la ingeniería del petróleo y en hidrogeología. Mientras que en las ciencias atmosféricas y en modelos que estudian la interacción suelo-atmosfera son usados generalmente solo para modelar el estado del sistema y actualizarlo, en la ingeniería del petróleo y en la ingeniería se han utilizado para actualizar tanto las variables de estado como los parámetros del sistema (Naevdal et al., 2005). Los filtros de Kalman de conjuntos se han aplicado satisfactoriamente para asimilar la altura piezométrica y así mejorar la capacidad de predicción de modelos de manera dinámica. (Chen y Zhang, 2006; Hendricks Franssen y Kinzelbach, 2008; Li et al., 2011; Zhou et al., 2011). Hendricks Franssen et al (2011) en su artículo publicado en Water Resources Reserch (WRR), utilizó información y modelos en tiempo real de la ciudad de Zurich (Suiza) para simular el flujo dando como resultado que la asimilación diaria de datos de alturas piezométricas por medio de filtros de Kalman de conjuntos logra una mejor caracterización de las alturas piezométricas que la obtenida por medio de calibración inversa con datos históricos pero sin asimilar nueva información. Li et al. (2011) mostró que los filtros de Kalman de Conjuntos utilizando variables normalizadas NS-EnKF (Zhou et al. 2011) pueden caracterizar distribuciones hidráulicas no multigaussianas aun y cuando se tiene un modelo de distribución erróneo como distribución de partida. Nowak y Hendricks Franssen (2012) utilizan los filtros de Kalman de Conjuntos para estimar parámetros a partir de una tomografía hidráulica de un campo logmultigaussiano de conductividades en 3D. Combinando los filtros de Kalman de conjuntos con la anamorfosis gaussiana (GA). Gráficamente, la anamorfosis consiste en deformar el histograma de los datos en un histograma no gaussiano, de modo que la variable transformada, denotada Y(x), tenga una distribución gaussiana estándar (media 0 y varianza 1). Los resultados obtenidos logran un mejor acercamiento en la predicción de flujo y transporte. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 19 3.1 Gestión de datos del GIS Se comenzó a trabajar con los datos del modelo FUA haciendo uso de los archivos de entrada que contienen la información necesaria para que el modelo estime las concentraciones y niveles piezométricos en los distintos nodos del modelo. Además se utilizó un proyecto GIS el cual contiene información útil del área de estudio, dicho proyecto fue creado en 1994, para ARCView3.2, así que lo primero fue, a partir de la información existente, crear un proyecto semejante en ArcGIS, el cual se fue actualizando con la información más reciente y además se le fue añadiendo información paulatinamente a medida que se iba obteniendo. Como se comentó en la introducción del documento se trabajó con el modelo del FUA 2004 que utiliza el CORE (CÓdigo para la simulación numérica de procesos de flujo de agua, transferencia de calor y solutos REactivos) para obtener las concentraciones y alturas piezométricas en cada uno de los nodos en los que esta discretizada el área de estudio. El modelo consta de 12254 elementos triangulares y 6218 nodos, derivados de la subdivisión de la malla utilizada en el modelo de 1994; esto con el fin de tener un mayor número de elementos en las zonas más cercanas a la fábrica de Uranio del FUA, así como en algunos otros lugares de interés, como son los lugares en donde las concentraciones tenían gran variabilidad entre las calculadas por el modelo y las observadas durante las campañas de muestreo. La figura 3.1 nos muestra la malla creada a partir de la información disponible en los archivos de entrada de CORE, la cual fue añadida al SIG del FUA. Posteriormente de obtener la malla se creó una cuadricula de 1344 filas por 660 columnas, con tamaño de celdas de 5x5 metros que abarca toda el área de estudio, iniciando en la coordenada UTM 30 401920 Este, 4208450 Norte y finalizando en la en la coordenada UTM 30 408640 Este, 4211750 Norte. Esto da como resultado una cuadricula de 887040 celdas. La Figura 3.2 nos muestra como la cuadricula abarca toda el área de estudio. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 20 Figura 3.1 Malla de elementos triangulares creada a partir de la información de entrada utilizada por el código CORE. Figura 3.2 Cuadricula de 887040 elementos que envuelve el área de estudio. Una vez creada la cuadricula, se obtuvieron los valores de conductividades hidráulicas de los puntos de muestreo de la red de Vigilancia. A partir de la información disponible en el SIG. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 21 En total se tienen 154 puntos con información de conductividades que fueron extraídos del SIG para utilizarlos como información de partida para el programa SGEMS. Figura 3.3 Mapa de localización de los 154 datos de Conductividades. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 23 4 Análisis geoestadístico 4.1 Modelo de continuidad espacial Para poder realizar un análisis Geoestadístico es necesario contar con información útil que nos ayude a comprender la variabilidad espacial de la propiedad que nos interesa conocer, ya sea obteniendo valores medidos de esta propiedad, o bien, a través de otra propiedad que esté correlacionada con la primera, en este caso se creó el campo de conductividades hidráulicas directamente del modelo FUA04 mediante el programa SGEMS introduciendo las conductividades previamente obtenidas de los puntos de observación de la red de vigilancia. La figura 4.1 nos muestra la distribución espacial de dichos puntos, así como el valor de la conductividad hidráulica para cada uno de los 154 puntos en m/d. Figura 4.1 Distribución espacial de los valores de conductividad hidráulica. Se muestra la distribución espacial de los valores de conductividad hidráulica. Cabe señalar que se tiene un aumento en la conductividad hacia el oeste de la zona de la FUA. En la figura se logra apreciar que se tiene una cierta continuidad cercana a los 120° azimut, por lo que se generó un variograma omnidireccional para apreciar mejor el alcance que se tiene en los EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 24 datos, además se generaron diversos variogramas direccionales variando el ángulo de búsqueda en 15 grados cada vez. En cuanto a los datos estadísticos de la muestra se realizó un histograma de frecuencias para conocer los valores de la media y la desviación típica de los datos, estos se observan en la figura 4.2. Primeramente llama la atención en el histograma el que existen una gran cantidad de valores altos de la muestra, al observar la función de distribución acumulada se aprecia que de manera general se sesga hacia valores menores de 200 m/d, se observa que los datos no están distribuidos de manera uniforme en el histograma, posiblemente debido a la gran cantidad de pozos de observación cercanos a las instalaciones de la fábrica. La desviación típica de la muestra es importante, mientras que la media llega a ser de 329.18 m/d. Figura 4.2 Histograma y curva de probabilidad acumulada. Se observa una gran variabilidad en los datos, que contribuye a afirmar la hipótesis de la heterogeneidad del medio mencionada en los trabajos previos. Volviendo un poco atrás a los variogramas mencionados el variograma omnidireccional muestra que se tiene un alcance cercano a los 3000 metros de distancia, a partir del cual comenzamos a buscar la dirección de máxima continuidad. La figura 4.3 nos muestra el variograma omnidireccional. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 25 Figura 4.3 Variograma omnidireccional Los variogramas direccionales tienen una gran variabilidad espacial, por lo que se buscó los ejes de mayor y menor continuidad, observándose que entre los 120 y los 135 grados se tenía la menor variabilidad espacial, en cuanto al eje de menor continuidad se observo que estaba entre los 30 y los 45 grados. Las figuras 4.4 y 4.5 nos muestran la gran variabilidad entre el variograma en la dirección de 30 grados y el de la dirección de 127 grados. El variograma de 30 grados parece más un efecto pepita, quizás porque la gran mayoría de datos están distribuidos en una banda con una orientación de 120 grados azimut, sin embargo, se logra observar en la figura 4.5 que los datos tienen una mayor correlación aparente en el eje de los 120 grados. El variograma de 120 grados además tiene una disminución de la variabilidad pasando los 3000 metros para volver a ascender nuevamente, Pyrcz y Deutsch detallan el fenómeno en su publicación haciendo hincapié a cierto ciclamiento en la variabilidad espacial de los datos que le da esa forma característica, conocido como efecto hueco (Pyrcz, Deutsch, 2003). EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 26 Figura 4.4 Variograma unidireccional en la dirección de los 30 grados azimut. Figura 4.5 Variograma unidireccional en la dirección de los 120 grados azimut. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 27 Tras ajustar un poco en busca de los ejes de máxima y mínima continuidad se determinó que el eje de máxima continuidad estaba en los 127° azimut, mientras que el de mínima continuidad era perpendicular a este último, tal y como se muestra en la figura 4.6. Figura 4.6 Ejes de máxima y mínima continuidad de los datos de conductividad hidráulica. Ya conocidos los ejes principales se definió un modelo de continuidad para los datos, el cual está definido por un efecto pepita de 10000 unidades y una contribución de 140000 unidades añadidas con un modelo de variograma esférico: 0(ℎ)= 23 ‖ℎ‖ 2 ! −125‖ℎ‖ !6 7 8 (4.1) El modelo resultante es el siguiente: 0(ℎ)=10000+1400000 () (ℎ) (4.2) Donde el primer termino corresponde al efecto pepita puro y el segundo termino es el aporte proveniente del variograma esférico con un alcance de 100000m y un rango de 5670m en el eje de máxima continuidad α = 125° y un rango de 2940m en el eje de mínima continuidad β= 35°. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 35 • Todos los subconjuntos de la función aleatoria, por ejemplo { W (Q),Q ∈ ] ⊂S} son también normales multivariados. • Tener covarianza cero (o correlación) hace que las variables sean independientes: Si _ { W (Q),W (Q’)} =0 las dos variables aleatorias W(Q) a W(Q’) no son solo no correlacionadas sino que también son independientes. • Todas las combinaciones lineales de las variables aleatorias que componen W(Q) (en modelos univariados) son normalmente distribuidas: b =∑d B CBD W (Q e ) Tiene distribución normal, ∀E, ∀ los pesos d B siempre y cuando Q e ∈ S (6.4) • Todas las distribuciones condicionales de todos los subconjuntos de la función aleatoria W(Q), teniéndose la realización de cualquier otro subconjunto, son distribuciones normales multivariadas. Por ejemplo la distribución condicional de las # variables aleatorias { Wf (Q’ X ) f=1,…,#} siempre que {Q’ X ∈S}, dada la realización a(Q B ) = a B , I=1,…,E es variada normal de orden #,∀#,∀Q’ g ,∀E, ∀Q e ,∀a B . En el caso en que #=1. Q’ h =Q ) , donde la variable aleatoria W(V i ) modela la incertidumbre de un valor no muestreado a(Q ) ), nos es de gran interés debido a que: la función condicionada de densidad acumulada W(Q ) ), dados los E datos a B , es normal y completamente caracterizada por: • La media o valor esperado, el cual es obtenido mediante los pesos del krigeado simple. • La varianza condicional del krigeado simple. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 36 Siendo solamente necesaria la resolución del sistema del sistema matricial del Krigeado Simple para obtener una función de densidad de probabilidad que nos permita obtener el valor: ;(Q)=  {W X (Q)} (6.5) Definido esto, en nuestro caso de estudio se eligió la simulación Gaussiana secuencial, donde cada variable es simulada secuencialmente de acuerdo con la función condicionada de densidad acumulada; siendo la condicionante todos los valores de los datos originales y todos los valores previamente simulados dentro de un vecindario de búsqueda determinado. El Software SGEMS fue utilizado debido a la manejabilidad que tiene y a que cuenta con un algoritmo para realizar la simulación gaussiana secuencial bastante estable (el cual fue desarrollado previamente para la librería G s TL y posteriormente, anexado al software), el programa realiza los siguientes pasos para cada realización: Primero obtiene la función de densidad de probabilidad acumulada, j k () del área en estudio. Se define un recorrido aleatorio con el cual se visita cada nudo de la malla (la cual en principio no importa si es regular o no, aunque se observó en la realización de este trabajo que no produce los mismos resultados) una vez, esto mediante un valor “semilla” que hace que se visiten los nodos como mínimo una vez, ya que si el nodo visitado no tiene datos cercanos se pasa a un segundo nodo y el primero será visitado posteriormente en espera de tener más datos disponibles. En cada nodo Q, se retiene un cierto número de valores (especificados previamente) condicionantes incluyendo tanto datos iniciales como valores previamente simulados. En el algoritmo de la simulación gaussiana secuencial se utiliza comúnmente un krigeado simple con el modelo de un variograma normal a fin de determinar la media y la varianza de la función condicionada de densidad acumulada de la función aleatoria j(Q) en el sitio Q. En nuestro caso utilizamos un krigeado ordinario a fin de poder asimilar el valor esperado EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 37 {W(Q)} no estacional para cada locación Q, existiendo asi mayor variabilidad espacial entre datos cercanos. Después de determinar la media y la varianza en el sitio Q, se toma un valor simulado a (U) (Q) de la función condicionada de densidad acumulada. Se añade el valor de a (U) (Q) a los datos iniciales. Se procede con la siguiente locación y se comienza de nuevo el proceso hasta que todos los nodos están simulados. Por último se convierten los valores normalizados {a (U) (Q),Q∈S} de vuelta a valores congruentes a los que tiene la variable original { (U) (V) = l m (a (U) (Q)) Q ∈ S}. Para ello se interpola o extrapola dependiendo de si el valor a convertir esta en la cola inferior o superior de la función, o bien en el centro de la distribución mediante distintos modelos posibles: Cola inferior Centro Cola Superior Modelo lineal Modelo lineal Modelo lineal Modelo exponencial Modelo exponencial Modelo exponencial A partir de valores tabulados A partir de valores tabulados A partir de valores tabulados Modelo Hiperbólico Tabla 6.1 Modelos de interpolación o extrapolación utilizados al transformar los valores. Para realizar las simulaciones para el área de estudio se condiciono un vecindario de búsqueda de forma elíptica con un semieje mayor de 3000m y un semieje menor perpendicular de 3000m, además se condicionó a tener un máximo de 40 valores con los cuales realizar la estimación del valor de cada celda. Se realizaron un total de 500 simulaciones las cuales arrojan valores esperados entre los 0.00001 y 1365 m/d, con una media de 444.5 m/d y una varianza de 44495 (m/d) 2 . EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 38 Las figuras 6.1, 6.2 y 6.3 nos muestran los resultados de 3 simulaciones distintas la simulación 22, la simulación 35 y la 43 respectivamente, además la figura 6.4 nos muestra el mapa con los valores esperados. Figura 6.1 Mapa de conductividades hidráulicas (m/d) resultado de la simulación 22. Figura 6.2 Mapa de conductividades hidráulicas (m/d) resultado de la simulación 35. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 39 Figura 6.3 Mapa de conductividades hidráulicas (m/d) resultado de la simulación 43. Figura 6.4 Mapa de los valores esperados de las conductividades hidráulicas (m/d). Cabe mencionar que la figura 6.4 es en sí muy parecida a la figura obtenida con el krigeado ordinario (figura 5.1), no obstante, en las otras 3 figuras se ve claramente una heterogeneidad en el terreno introducida en cada simulación lo que le da ese coloreado punteado con un aspecto similar al que crea un crayón de cera y que es más parecido a la realidad. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 40 La figura 6.5 nos muestra el campo de Varianza entre las 500 simulaciones para cada uno de los elementos dentro de la zona de estudio. Cabe señalar que la mayor varianza está en la parte sudoriental del mapa, pues coexisten datos con conductividades muy altas a poca distancia de los valores más bajos de toda el área de estudio. Figuras 6.5 Varianza Condicional en del Área de estudio. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 41 7 CORE 7.1 Aspectos teóricos del código El código utilizado (CORE2D V2.0) es una versión mejorada del código CORE-LE2D VO (Samper et al., 1998), el cual es un código de elementos finitos que, como bien se dijo anteriormente resuelve las ecuaciones de flujo, transporte de solutos y transferencia de calor en medios con contornos irregulares y propiedades físicas y geoquímicas no uniformes. 7.1.1 CORE acuífero confinado El movimiento del agua en un medio poroso es gobernado por la Ley de Darcy, la cual relaciona el flujo de agua n con el gradiente de presión o de la siguiente manera: n= f(∇o+qr∇) (7.1.1.1) Donde: q es la densidad del agua (ML -3 )  es la viscosidad dinámica (MT -1 L) f es el tensor de la permeabilidad intrínseca (L 2 ) r es la aceleración de la gravedad (LT -2 ) Si la densidad del agua sufre cambios despreciables, entonces la ley de Darcy puede escribirse en términos de altura piezométrica como sigue: n= −# ∇ℎ (7.1.1.2) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 42 Donde # se obtiene de la siguiente manera: #= f ρ g  (7.1.1.3) Siendo # el tensor de conductividades hidráulicas (LT -1 ) Y la altura piezométrica se define como: ℎ= p qr+  (7.1.1.4)  es la altura medida a partir de un nivel de referencia (L) Combinando la ley de Darcy con la ecuación de balance de masas tenemos que: ∇∙(#∇ℎ)+&=w x yℎ yz (7.1.1.5) Luego: & representa el termino fuente/sumidero [L 3 / L 3 /T] w x es coeficiente de almacenamiento específico [L -1 ] Como nos encontramos en un acuífero confinado podemos integrar el flujo dentro de todo el espesor del acuífero L, L está definido como la diferencia entre la elevación del fondo y el techo del acuífero llamadas  U y  K respectivamente. La ecuación resultante es: ∇∙(%∇ℎ)+"=w yℎ yz (7.1.1.6) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 43 Donde: " es el termino fuente/sumidero por unidad de superficie [L 3 / L 2 /T] % es el tensor de transmisividades [L 2 T -1 ] y S es el coeficiente de almacenamiento [-] que se define como: T(x,y)=~#(,a,)$=#L   (7.1.1.7) S(x,y)=~w x (,a,)$=w x L   (7.1.1.8) w x  y # son los promedios de las conductividades hidráulicas y los coeficientes de almacenamiento en el sentido vertical. 7.1.2 CORE acuífero libre En un acuífero libre el límite superior coincide con la lamina de agua, en estos casos la trasmisividad % está dada por #(ℎ− K ) y depende de la piezometría del acuífero. Ademas, el coeficiente de almacenamiento se calcula de la siguiente manera: w=w x (ℎ− K )+w  (7.1.2.1) Donde sy es el rendimiento especifico [-]. Se debe además complementar la ecuación de flujo con condiciones iniciales apropiadas, como podría ser el caso de satisfacer las condiciones del estado estacionario: ∇∙(%∇ℎ ) )+" ) =0 (7.1.2.2) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 44 o bien conociendo a priori el estado del sistema, tomando ademas en cuenta las diferentes condiciones de contorno posibles: • Dirichlet conociendo la altura piezométrica en el contorno Γ  • Neuman teniendo un flujo prescrito en el contorno Γ  • Cauchy con un flujo dependiente de la altura piezométrica en Γ 7 7.1.3 CORE zona no saturada Si se tiene un flujo en un medio con saturación variable, la altura piezométrica se escribe como la suma de la altura de presión Ψ y la elevación z. Siendo Ψ positiva en medios saturados y negativa en regiones parcialmente saturadas. Por tanto se puede calcular el flujo por medio de la ley de Darcy y la ley de conservación de la masa, tomando en cuenta la altura de presión Ψ: ∇∙#  #∇(Ψ+)+&=yw - yΨ+w - w x yΨ yz (7.1.3.1) De esta ecuación el término de la conductividad hidráulica es el producto de la conductividad relativa K  [-] y la conductividad saturada K. El término que multiplica a la derivada parcial de la altura de presión con respecto al tiempo es la capacidad de almacenamiento, el primer término es importante en zonas no saturadas, mientras que el segundo término es relevante en zonas saturadas. Ambos términos anteriormente señalados, así como la conductividad relativa K  dependen de la altura de presión, lo que hace a esta ecuación extremadamente no lineal, tomando en cuenta que además se debe calcular el contenido de humedad de la siguiente manera: w - =∅ (7.1.3.2) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 51 El superíndice § denota el número de elemento que contiene los nodos ? y E, [ son segmentos adyacentes al nodo E que juegan un papel de frontera. La expresión de S C  , ¦ C  , $ C , S C x y $ Cx son: S C  =1 4Δ © ª%   L   L C +%   (L C «   +L   « C )+%   « C «   ¬ (7.2.1.16) ¦ C  =w  Δ  12 C Donde:  C =®2 [\ ?=E 1 [\ ? ≠E (7.2.1.17) $ C ="  ∆  3 (7.2.1.18) S C x =I x  x 6 C (7.2.1.19) $ Cx =(I x ¡ x +¢ x ) x 2 (7.2.1.20) Donde ∆  es el área del elemento, L   y «   son coeficientes de los elementos que dependen de la geometría del elemento. La ecuación 2.3.11 pueden ser escritas en forma de matriz como: Sℎ+¦$ℎ $z=$ (7.2.1.21) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 52 Donde S es una matriz de conductancias con   valores, ¦ es una matiz de capacitancias con   valores, $ es un vector de términos independientes en forma de columna, y ℎ es un vector en forma de columna de alturas piezométricas. Este sistema puede resolverse por medio del método de diferencias finitas resultando en: Sℎ X° + (1-)ℎ X ]+D ±²³´m±² µ¶ =$ X° (7.2.1.22) Donde: Δz=z X° −z X , es decir, es el resultado de aplicar el método a la derivada respecto a z $z 0≤≤1,ℎ X a ℎ X° son los vectores de los niveles piezométricas en los tiempos z X y z X° siendo el esquema de integración explicito =0 implícito =1 o con un esquema CrankNicholson =0.5 Asi el sistema se resuelve secuencialmente una vez conocido ℎ X se procede a resolver para el tiempo ℎ X° mediante la expresión: S+¦ ∆zℎ X° =$ X° +£(−1)S+¦ ∆z¤ℎ X (7.2.1.23) La trasmisividad % es el producto del espesor del acuífero L y la conductividad hidráulica #. Si el acuífero es confinado, el espesor del acuífero es constante con respeto al tiempo, variando solamente de entre elementos, CORE asume un valor de b=1 para el flujo vertical en acuíferos confinados. Para el flujo en 3 dimensiones con simetría axial en acuíferos confinados CORE calcula L  como la longitud promedio de la coordenada radial del elemento: L  =2¹   +  + 7 3 (7.2.1.24) Siendo   ,  ,  las coordenadas radiales de cada elemento. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 53 Si el acuífero es libre, entonces el espesor saturado varia tanto con la distancia como con el tiempo, por lo que es necesaria la utilización de métodos iterativos para resolverlo. El espesor saturado promedio es actualizado a partir de la iteración anterior de acuerdo a la siguiente expresión: L  =ℎ  +ℎ  +ℎ 7 3− K (7.2.1.25) Donde: ℎ  son las alturas piezométricas en los nodos y  K es la base del elemento §. Este cálculo iterativo se repite hasta que se cumpla la siguiente condición: max ¼ ℎ ¼X°,x° −ℎ ¼X°,x ℎ ¼X°,x° +ℎ ¼X°,x ≤d \=1,2,…, (7.2.1.26) Donde f y [ indican el intervalo de tiempo y el nivel de iteración respectivamente; y d, es un parámetro especifico que determina tolerancia del error en la altura piezométrica. En el primer intervalo y la primera iteración L es calculado a partir de la distribución inicial de ℎ, resolviendo primero el sistema en estado estacionario. Una vez definido ℎ se resuelven los componentes  e a del vector de velocidades de acuerdo a la ley de Darcy: n  =−½#   @ℎ ¼ L ¼ 7 ¼D +#   @ℎ ¼ « ¼ 7 ¼D ¾ (7.2.1.27) n  =−½#   @ℎ ¼ « ¼ 7 ¼D +#   @ℎ ¼ L ¼ 7 ¼D ¾ (7.2.1.28) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 55 7.2.2 Solución del flujo en la zona de transición Los métodos numéricos para la solución del flujo en zonas con saturación variable son similares a las del flujo en acuíferos, solo que dependen de la altura de presión. La ecuación del flujo (7.2.1.12), puede ser discretizada por medio del método de elementos finitos para obtener: @Ψ ¿  D A ¿Á +@dΨ ¿ dt  D D ¿Á =d Á +B Á (7.2.2.1) Donde: A ¿Á , D ¿Á y d Á son las ecuaciones 7.2.1.13, 7.2.1.14, y 7.2.1.15, en donde A ¿Á © , D ¿Á © puede ser expresado como: A ¿Á © =K © C ¿Á © (7.2.2.2) Donde C ¿Á © tienen la siguiente expresión: C ¿Á © =1 4Δ © ªK ÆÆ © b ¿ © b Á© +K ÆÈ © (b Á© c ¿ © +b ¿ © c Á© O+K ÈÈ © c Á© c ¿ ©  (7.2.2.3) Y D ¿Á © es calculado de la siguiente manera: C ¿Á © =S © Δ © 12 C =®2 [\ ?=E 1 [\ ? ≠E (7.2.2.4) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 56 De Acuerdo con la ecuación 7.2.1.12, para flujo no saturado, la capacidad de almacenamiento S © esta dada por: w  =  yw - yψË  +w -  w x (7.2.2.5) El término B Á de la ecuación 7.2.2.1 se deriva del término gravitatorio, y tiene la siguiente expresión: ] C = − #  2(#   « C +#   L C ) (7.2.2.6) Donde  e a son las coordenadas vertical y horizontal respectivamente. El valor promedio en cada elemento de la conductividad hidráulica relativa, la saturación de agua y los demás términos derivados son evaluados con la altura de presión promedio de los tres nudos, o bien: #  =#  (Ì  ) (7.2.2.7) w Í  =w Í (Ì  ) (7.2.2.8) yw - yψË  =yw - yψ(Ì  ) (7.2.2.9) Donde: Ì  =Ì  +Ì  +Ì 7 3  (7.2.2.10) EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 57 Los términos de la ecuación 7.2.2.1 correspondientes a la conductancia A ¿Á , capacidad de almacenamiento D ¿Á y el término independiente B Á dependen de la altura de presión, lo que ocasiona que sea una ecuación no lineal en gran medida. Para resolverla se utiliza el método iterativo de Newton-Rapshon modificado por Galarza (1993) mediante la formulación 3 de este último. A partir de esto se modifica la ecuación 7.2.2.1 de la siguiente manera: S X°Î Ì X°Î +¦ X°Î Ì X° −Ì X Δz =$ X° +] X°Î (7.2.2.11) Donde: S X°Î =S(Ì X°Î ) (7.2.2.12) ¦ X°Î =¦(Ì X°Î ) (7.2.2.13) Ì X°Î =Ì X° +(1−)Ì X (7.2.2.14) Y el valor de  debe ser 0≤ ≤1 la ecuación 7.2.2.11 puede ser reescrita de la siguiente manera: =S X°Î Ì X°Î +¦ X°Î Ì X° −Ì X Δz −$ X° −] X°Î =0 (7.2.2.15) Donde Ì X es conocida, Ì X° es desconocida, y Ì X°Î es una combinación de ambos. El método iterativo de Newton-Raphson está basado en las series de expansión de Taylor. Por lo que al reacomodar estas expansiones el sistema de ecuaciones queda de la siguiente manera: EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 58 @y C yÌ ¼X°  ¼D ΔÌ ¼X° =− C E=1,2,…, (7.2.2.16) Donde N es el número de nodos. En forma matricial la ecuación se reduce a la siguiente forma: ÏΔÌ X° =− (7.2.2.17) Donde Ï es la matriz Jacobiana que contiene los términos diferenciales de la ecuación 7.2.2.16 y  es el vector de residuos. Una vez se han calculado los valores de ΔÌ X° estos son utilizados para actualizar Ì X° de acuerdo a: Ì X°,x° =Ì X°,x +ΔÌ X° (7.2.2.18) Ì X°,+ =Ì X (7.2.2.19) Donde el subíndice [ es el número de iteración. El proceso iterativo es realizado hasta que el valor absoluto del incremento relativo de la altura de presión para todos los nodos es menor que el especificado con un valor de tolerancia en la convergencia (en general un número muy pequeño): max ¼ 5ÐΔÌ X° Ì X°,x Ð6≤d \=1,2,…, (7.2.2.20) La matriz Jacobiana Ï se obtiene tomando derivadas de la ecuación 7.2.2.15 con respecto a la incógnita ΔÌ X° , y tomando en cuenta las ecuaciones 7.2.2.12 a 7.2.2.14, o bien: EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 59 Ï=½y S f+ y Ì f+ Ì f+ ¾+ ∆zÑy ¦ f+ y Ì f+ Ò Ì f+1 −Ì f ÓÔ+ S f+ +¦ f+ ∆z −y ] f+ y Ì f+ (7.2.2.21) La matriz Jacobiana Ï es una matriz en banda no simétrica con  que puede resolverse al multiplicar elemento por elemento (es decir, teniendo 2 matrices (b;W) el primer elemento de la matriz b por el primer elemento de la matriz W) y luego reasignando las contribuciones de cada elemento de una manera similar a como estaba la matriz de conductancia S original. 7.3 Datos de entrada al modelo CORE El modelo CORE consta de 3 archivos de entrada, escritos en archivos de texto que, actuando en conjunto, son leídos por el código para resolver problemas de flujo, transporte de solutos reactivos y transporte de calor; combinando el flujo con el transporte de solutos reactivos o con el transporte de calor o bien todos en conjunto. El archivo ROOT.inp (siendo ROOT el nombre que se utilizara para esa simulación en especifico) tiene los parámetros de control generales del problema, como es el numero de parámetros, numero de iteraciones, los pesos iniciales con los que se comienza a ajustar la función de flujo, el valor mínimo y máximo de los parámetros, las coordenadas y nombres de los puntos de muestreo y observación, el número de especies químicas con que lidiara el código, etc. El archivo ROOT_che.inp es un archivo en donde se detallan valores de distintos complejos acuosos, gases, cationes de intercambio y complejos superficiales, así como las distintas especies químicas que tratara el problema. Este archivo debe contener la temperatura inicial del problema (por defecto 25 grados centígrados) y en general los parámetros que controlan el transporte de reactivos. También contiene los tipos de agua de contorno y de recarga con que tratara el sistema, así como la función de tiempo que corresponde a cada tipo de agua. El archivo ROOT_tra.inp es el archivo en donde se recoge toda la información espacial del modelo, coordenadas de nodos, numero de nodos, numero de elementos, numero de zonas de material que se encuentra en la zona, tipo de acuífero, condiciones iniciales y de contorno. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 60 También en este archivo se detalla la información de la discretización temporal del problema, las funciones de flujo y de transporte que utilizara el modelo, el número y nombre de cada uno de los archivos de salida que generará el modelo (máximo 21), los datos que tomara como criterio de convergencia el modelo (el cual se detalla en la sección anterior), las variables que controlan la creación de la matriz Jacobiana, los nodos de los que se van a imprimir resultados en los archivos de salida, etc. Además el código lee datos de valores observados, los espesores del acuífero, los datos de concentraciones de entrada, etc. Especificados en archivos de texto con la extensión .dat y deben ser creados para cada problema en especifico. No obstante dentro de los archivos .dat se encuentran 2, el archivo Masterte.dat y el archivo master25.dat que son los que controlan todas las reacciones e iteraciones de compuestos químicos con los que puede trabajar CORE, y que, no se contemplan en el alcance que tiene este estudio. Para mayor detalle véase el manual de usuario de CORE V 2.0 (ENRESA 200). El modelo CORE ha sido reestructurado desde el código fuente para poder utilizar una cantidad mayor de elementos que los que se tenían contemplados en el modelo de 1994, así como también se ha utilizado una mayor cantidad de zonas de material que las que tenía el modelo de 2004. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 67 8 Filtro de Kalman de conjuntos 8.1 Filtro de Kalman de conjuntos aplicado al FUA Para medir la incertidumbre del modelo en la estimación de los niveles piezométricos se hizo uso de los filtros de Kalman. El programa creado se puede subdividir en tres pasos: PASO 1. Modelo de predicción. La ecuación de flujo es resuelta en por medio de CORE, y que es representada por la ecuación: W X =(b Xm ,W Xm ) (8.1.1) Donde W X es el estado del sistema (nivel piezométrico) en el tiempo z X ,  representa el modelo de flujo (incluidas las condiciones de contorno, recargas y parámetros), y b Xm denota los parámetros del modelo (conductividades hidráulicas) después de la última actualización del modelo. Por ultimo W Xm es el estado previo del sistema. PASO 2. Análisis. Utilizando los datos de alturas piezométricas observadas se actualizan los parámetros del modelo. 1 Se construye el conjunto de vectores Ψ Ù , que incluye los parámetros b y la predicción del estado del sistema W. Este vector se puede dividirse en tantos miembros como realizaciones se tengan en el conjunto mediante: Ψ X,Ú =£bW¤ X,Ú (8.1.2) Siendo Û el numero de miembros del conjunto en el tiempo z X . EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 68 2 El conjunto de vectores es actualizado, realización por realización, asimilando las observaciones NW X)Kx O: Ψ Ù,Ü Ý =Ψ Ù,Ü © +G Ù (Ψ Ù,Ü ßàá +ϵ−HΨ Ù,Ü © ) (8.1.3) Donde el superíndice ! y § denota análisis y estimación, respectivamente; ϵ es un vector de error de observación aleatorio; H es un interpolador lineal que interpola la altura piezométrica estimada a la ubicación donde se encuentra la observación más cercana. En el caso de H se utiliza una función que sigue la siguiente forma: La ganancia del filtro de Kalman G Ù esta dada por: G Ù =P Ù© H å (HP Ù© H å +R Ù ) m (8.1.4) Donde R Ù es la matriz de covarianza de errores de medida. P Ù© contiene las covarianzas entre los diferentes componentes del vector de estado y puede ser estimado a partir de los resultados del conjunto estimado: P Ù© ≈EèNΨ Ù,Ü © −Ψ  Ù,Ü © ONΨ Ù,Ü © −Ψ  Ù,Ü © O å é ≈@NΨ Ù,Ü © −Ψ  Ù,Ü © ONΨ Ù,Ü © −Ψ  Ù,Ü © O å N  ëì ÜD (8.1.5) Donde N  es el número de realizaciones en el conjunto, y la barra denota el promedio del conjunto. PASO 3. El estado actualizado del sistema se convierte en el estado actual y el proceso comienza nuevamente. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 69 8.2 Descripción del programa El programa toma como input los archivos de salida de niveles piezométricos ROOT_HH1.out y ROOT_HH2.out provenientes de CORE, así como el conjunto de archivos ROOT_NRE.inp donde ROOT es el nombre del problema que se está tratando, y NRE es el numero de realización del campo de conductividades que utiliza el archivo de entrada en CORE para resolver las ecuaciones de flujo. Primeramente toma un archivo de texto con cada una de las realizaciones de conductividades hidráulicas del conjunto para modificar un archivo base, y crear un conjunto de archivos de entrada para CORE. Luego busca el primer valor de observaciones que se tenga en un archivo de texto que contiene los resultados de las campañas de medición de niveles piezométricos, enseguida compara la fecha de la observación con el intervalo de tiempo en CORE para modificar el conjunto de archivos de entrada de CORE y así, calcular hasta la fecha correcta. Se realiza una serie de cálculos en CORE para cada uno de los archivos de entrada. Esto genera un conjunto de conductividades hidráulicas, así como una serie de niveles piezométricos medidos en el intervalo temporal que se asignó a CORE para cada uno de los pozos de la red de vigilancia. Luego se introduce el conjunto completo en la subrutina que crea los distintos vectores que utilizan los filtros de Kalman de conjuntos, se genera la ganancia del Filtro de Kalman mediante las matrices de covarianzas que utiliza la metodología de los filtros de Kalman de Conjuntos, se actualiza el sistema y posteriormente se lee nuevamente el archivo de observaciones para saber si existen más observaciones disponibles, de ser asi entonces se actualizan los archivos de entrada a CORE con los nuevos valores de conductividades y se añade el intervalo de tiempo que se va a calcular para volver a resolver el sistema con CORE, esto se repite hasta que todas las observaciones sean tomadas en cuenta. El diagrama de flujo que representa el proceso anterior se observa en la figura 8.2.1. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONE Figura 8.2.1 Diagrama de flujo EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONE S DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS Figura 8.2.1 Diagrama de flujo del programa de F iltros de Kalman de Conjuntos código CORE. S DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE 70 iltros de Kalman de Conjuntos añadido al EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 71 9 Resultados y conclusiones 9.1 Resultados Dentro de los resultados obtenidos se aprecia que el programa va ajustando las conductividades para obtener un mejor ajuste de manera global en el modelo. Las figuras 9.1.1 a 9.1.7 muestran el ajuste gradual para los pozos de control PC1, PC2, PC3, PC5, PC6, así como los pozos 505 y 611, para la realización 2. Para cada una de las graficas se muestran los resultados que va obteniendo el filtro para cada intervalo de tiempo desde el intervalo 7 hasta el 96 con un incremento de 6 pasos de tiempo. Hay que recordar que cada intervalo de tiempo es de 15 días, por lo que la simulación 7 es el día 105 y la 13 el 195 y así sucesivamente. Figura 9.1.1 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC1 realización 2. 193 193.5 194 194.5 195 195.5 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC1 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 72 Figura 9.1.2 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC2 realización 2. Figura 9.1.3 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC3 realización 2. 190 191 192 193 194 195 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC2 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 192 192.5 193 193.5 194 194.5 195 195.5 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC3 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 73 Figura 9.1.4 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC5 realización 2. Figura 9.1.5 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el PC6 realización 2. 193.5 194 194.5 195 195.5 196 196.5 2002 2004 2006 2008 2010 2012 Altura Piezométrica m Años PC5 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 193 193.5 194 194.5 195 195.5 196 196.5 197 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC6 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 74 Figura 9.1.6 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el pozo 505 realización 2. Figura 9.1.7 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos para el pozo 611 realización 2. 193 193.5 194 194.5 195 195.5 196 196.5 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años 505 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 191.5 192 192.5 193 193.5 194 194.5 195 195.5 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años 611 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 75 En las siguientes figuras se detallan los mismos pozos de las figuras 9.1.1 a 9.1.6 pero para las realizaciones 3 y 9. Estas se muestran en las figuras 9.1.8 a 9.1.19 Figuras 9.1.8 a 9.1.13 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos con la realización 2. 192.5 193 193.5 194 194.5 195 195.5 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC1 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 190 191 192 193 194 195 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC2 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 192 192.5 193 193.5 194 194.5 195 195.5 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC3 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 193.5 194 194.5 195 195.5 196 196.5 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC5 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 193 193.5 194 194.5 195 195.5 196 196.5 197 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC6 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 193 193.5 194 194.5 195 195.5 196 196.5 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años 505 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 76 Figuras 9.1.14 a 9.1.19 Ajuste de conductividades realizado mediante Filtros de Kalman de conjuntos con la realización 9. 193 193.5 194 194.5 195 195.5 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC1 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 190 191 192 193 194 195 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC2 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 192 192.5 193 193.5 194 194.5 195 195.5 196 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC3 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 193.5 194 194.5 195 195.5 196 196.5 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC5 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 193 193.5 194 194.5 195 195.5 196 196.5 197 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años PC6 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 193 193.5 194 194.5 195 195.5 196 196.5 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 Altura Piezométrica m Años 505 observados sim 7 sim 13 sim 25 sim 31 sim 38 sim 43 sim 50 sim 56 sim 62 sim 67 sim 73 sim 78 sim 84 sim 90 sim 96 EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 83 9.2 Análisis de incertidumbre Se realizó un análisis de incertidumbre a partir de los 8 campos de conductividades que ajustan mejor los valores observados de la simulación 73, las conductividades hidráulicas para los pozos de control PC1 al PC6 se observan en las siguientes figuras, tomando como caso base el valor de conductividad que tenia la zona donde se encuentran cada uno de estos pozos en el modelo del FUA de 2004. La tablas 9.2.1 a nos muestran a que zona corresponde cada uno de los pozos de control, así como las conductividades correspondientes. Pozo Zona Conductividad PC1 2515 4.55E+01 PC2 2544 3.59E+02 PC3 3492 5.64E+01 PC4 2455 4.11E+03 PC5 2302 1.13E+02 PC6 6472 4.60E+06 Tabla 9.2.1 Pozos de control y zonas correspondientes para la realización 0. Pozo Zona Conductividad PC1 2515 4.70E+01 PC2 2544 3.70E+02 PC3 3492 6.88E+01 PC4 2455 4.08E+03 PC5 2302 2.90E+01 PC6 6472 5.00E+00 Tabla 9.2.2 Pozos de control y zonas correspondientes para la realización 1. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 84 Pozo Zona Conductividad PC1 2515 7.34E+01 PC2 2544 1.99E+02 PC3 3492 2.13E+01 PC4 2455 3.12E+03 PC5 2302 9.23E+01 PC6 6472 0.00E+00 Tabla 9.2.3 Pozos de control y zonas correspondientes para la realización 2. Pozo Zona Conductividad PC1 2515 2.12E+01 PC2 2544 3.74E+01 PC3 3492 6.31E+01 PC4 2455 5.56E+01 PC5 2302 6.48E+01 PC6 6472 2.11E+02 Tabla 9.2.4 Pozos de control y zonas correspondientes para la realización 3. Pozo Zona Conductividad PC1 2515 2.32E+01 PC2 2544 3.45E+04 PC3 3492 5.90E+01 PC4 2455 2.35E+04 PC5 2302 1.12E+02 PC6 6472 2.01E+05 Tabla 9.2.5 Pozos de control y zonas correspondientes para la realización 4. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 85 Pozo Zona Conductividad PC1 2515 1.73E+01 PC2 2544 1.43E+01 PC3 3492 1.07E+03 PC4 2455 5.00E+00 PC5 2302 1.86E+01 PC6 6472 3.20E+04 Tabla 9.2.6 Pozos de control y zonas correspondientes para la realización 5. Pozo Zona Conductividad PC1 2515 5.00E+00 PC2 2544 5.00E+00 PC3 3492 5.00E+00 PC4 2455 5.00E+00 PC5 2302 1.33E+01 PC6 6472 2.59E+03 Tabla 9.2.7 Pozos de control y zonas correspondientes para la realización 6. Pozo Zona Conductividad PC1 2515 2.43E+04 PC2 2544 8.76E+07 PC3 3492 1.22E+03 PC4 2455 5.00E+00 PC5 2302 5.10E+02 PC6 6472 7.91E+02 Tabla 9.2.8 Pozos de control y zonas correspondientes para la realización 7. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 86 Pozo Zona Conductividad PC1 2515 7.25E+00 PC2 2544 8.20E+01 PC3 3492 1.11E+02 PC4 2455 4.31E+05 PC5 2302 1.02E+02 PC6 6472 4.99E+06 Tabla 9.2.9 Pozos de control y zonas correspondientes para la realización 8. Pozo Zona Conductividad PC1 2515 4.21E+01 PC2 2544 4.07E+02 PC3 3492 3.16E+01 PC4 2455 1.93E+03 PC5 2302 5.52E+03 PC6 6472 0.00E+00 Tabla 9.2.10 Pozos de control y zonas correspondientes para la realización 9. Se concluye que los valores observados y calculados en algunas zonas tienen una gran variabilidad, por lo que aunque el uso de Filtros de Kalman de Conjunto ajusta el modelo para que reproduzca los valores observados lo mejor posible, no logra ajustarlos de manera congruente, creando así una varianza significativa en las conductividades. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 87 9.3 Conclusiones El presente trabajo se ha realizado para implementar el filtro de Kalman de conjuntos a un caso real y medir la incertidumbre que presentan los niveles piezométricos al introducir variabilidad en la heterogeneidad del terreno, esto mediante una serie de realizaciones de distintos campos de conductividades. El análisis geoestadístico corroboró las hipótesis acerca de la dirección y variabilidad espacial que se hicieron al crear los distintos modelos de la FUA. Se observó que, tal y como detalla Samper et. al (2004) se tienen bien definidas una serie de zonas con conductividades hidráulicas semejantes, y que, en su mayoría tienen un eje de máxima continuidad cercano a los 125°. Se generaron una serie de 500 simulaciones las cuales siguen la misma tendencia que los datos originales. Se observó que el solo hecho de crear campos de conductividades hidráulicas con variabilidad espacial a menor escala no es suficiente como para representar adecuadamente los niveles piezométricos observados utilizando el modelo de la FUA del año 2004. También se observó que, se tiene un cierto grado de incertidumbre al introducir distintos valores de conductividades que, si bien es cierto modifican las alturas piezométricas, aumentando o disminuyendo el nivel en ciertos puntos, no logran adecuarse a los valores observados pues la tendencia global no se ve alterada. Se realizó un programa que mediante el uso de Filtros de Kalman de Conjuntos, asimila información en cada paso iterativo modificando las conductividades para reproducir mejor los valores observados. De esto cabe señalar dos aspectos: Primero que los filtros de Kalman de conjuntos son una potente herramienta que ha demostrado su eficacia en distintos campos de aplicación, y que, como en este estudio logra reproducir de manera más eficiente el patrón observado al asimilar información reciente. En definitiva los niveles piezométricos son más parecidos con el uso de filtros de Kalman de conjuntos que sin ellos. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 88 El segundo aspecto a señalar es que, aunque los filtros de Kalman ajustan las conductividades de manera iterativa para que reproduzcan mejor los valores observados, se tienen conductividades hidráulicas muy elevadas en algunos puntos. Esto debido a que los valores de las distintas funciones de flujo del modelo del FUA del año 2004 no son suficientes para representar los valores observados hasta el año 2010, siendo medianamente ajustable la piezometría hasta el año 2005. La calibración que se hizo en el año 2004 a partir de la información disponible desde 1977 hasta 2003 y con la que se realizó la predicción para el periodo 2004-2020 no es suficiente para poder ajustar de manera satisfactoria el nuevo conjunto de valores observados, por lo que será necesario una nueva calibración del modelo para este periodo. Cabe señalar que, si bien no se logro ajustar las alturas piezométricas para el periodo completo de observaciones, si se observo la gran diferencia que supone el uso de los filtros de Kalman de conjuntos y la ventaja que se tiene al agregar información reciente. Este trabajo también sirve de antesala para en un futuro próximo verificar las herramientas creadas para el cálculo de filtros de Kalman de conjuntos en un periodo anterior con funciones de flujo ajustadas. Posteriormente se incluirán las funciones de flujo a los parámetros de ajuste de los filtros de Kalman de conjuntos. Se extenderá el uso de las herramientas creadas, para poder ajustar las concentraciones en el modelo. Finalmente mediante las herramientas creadas es posible observar la inflación de la covarianza en el área de estudio, así como la localización de áreas con una variabilidad muy grande. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 89 10 Bibliografía Almendral-Vazquez A. Randi Syversveen A. (2006). The ensemble Kalman filtertheory and applications in oil industry. Norsk Regnesentral. SAND/05/06. Chen Y., Zhang D. (2006). Data assimilation for transient flow in geologic formations via ensemble Kalman filter. Adnvances in Water Resources 29 (8), 1107-1122. Cox H. (1964). On the estimation of state variables and parameters for noisy dynamic systems. IEEE Transactions on Automatic Control 9 (1), 5-12. Deutsch C.V. , Journel A.G. (1998). GSLIB Geostatistical softwre Library and Users's Guide. S.E. Oxford University Press. New York. ENRESA. (2004). Actualización del modelo de flujo y transporte de Andújar. Tomos 1 a 4. Evensen G. (1994). Sequential data assimilation with a nonlinear quaisgeostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 99 (C5), 10, 143-10, 162. Evensen G. (2003). The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean Dynamics 53, 343-367. Gomez M. (2011). Uncertainty Analisis of groundwater flow and solute transport model FUA04. MSc Proyect, UPV. Haugen et al. (2008). History matching using the ensemble Kalman filter on a North Sea field Case. SPE Journal 13 (4), 382-391. Hendricks Franssen et al. (2011). Operational real-time modeling with ensemble Kalman filter of variably saturated subsurface flow including stream-aquifer interaction and parameter updating. Water Resources Reserch, Vol. 47 W02532, 20 PP., 2011 Hendricks Franssen, H., Kinselbach, W. (2008). Real-time groundwater flow modeling with the Ensemble Kalman Filter: Joint estimation of states and parameters and the filter inbreeding problem. Water Resources Reserch 44 (9), W09408. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 90 Isaaks E., Srivastava R. (1989). An Introduction to applied Geostatistics. Oxford University Press. New York. Kalman R. (1960). A new Aproach to Linear Filtering and Prediction Problems, Transactions of the ASME-Journal of Basic Engineering, Vol. 82 , Series D, 35-45. Kepert J. , Sun X. , Steinle P. (2003). Atmospheric data assimilation using the Ensemble Kalman Filter at BMRC. Bureau of Meteorology Research Centre. Kopp R., Orford R. (1963). Linear regression applied to system identification for adaptive control systems. AIAA Journal 1 (10), 2300-2306. Li et al. (2011). Jointly Mapping Hydraulic Conductivity and Porosity by Assimilating Concentration Data via Ensemble Kalman Filter. Naevdal et al. (2005). Reservoir monitoring and continuous model updating using ensemble Kalman filter. SPE Journal 10 (1), 66-74. Pyrcz M.J. Deutsch C.V. The whole story on the Hole Effect. Centre for computational Geostatistics, University of Alberta. Remi N. (2001). GsTL: The Geostatistical Template Library in C++. A report submitted to the department of petroleum engineering of Stanford University for the Degree of Master of Science. Remi N. (2004). Geostatistical Earth Modeling Softwarte: User's Manual. Samper et al. (1993). Calibración de los modelos de flujo y transporte de contaminantes en el agua subterránea.(Condición A.9.1, orden 1-2-1991). Informe elaborado para ENRESA. Samper et al. (1988). Application of an automatic calibration technique to modelling an alluvial aquifer. IAHS Publ. 195. 1990. Samper et al. (1991). Revision del modelo de Flujo y Transporte de contaminates en las aguas suberraneas del acuifero de la FUA. Informe CIMNE, IT-45 Samper et al. (2000). CORE 2D. A code for non-isothermal water flow and reactive solute transport. Users manual version 2. ENRESA. EVALUACIÓN DE LA INCERTIDUMBRE DE LAS PREDICCIONES DE LOS NIVELES PIEZOMÉTRICOS DEL ACUÍFERO DE LA FÁBRICA DE ANDÚJAR MEDIANTE EL USO DE FILTROS DE KALMAN DE CONJUNTOS 91 Samper J. Pisani B. (1998). CORE_LE_2D V0 Users Manual. ENRESA. Barcelona Schöninger A. Nowak W. Hendricks Franssen H.J. (2012). Parameter estimation by ensemble Kalman filters with transformed data: Approach and application to hydraulic tomography. Water Resources Research, Vol. 48, W04502, 18PP. Skjervheim et al. (2007). Incorporating 4d seismic data in reservoir simulation models using ensemble Kalman filter. SPE Journal 12 (3), 282-292, 95789-PA. Sun A.,Morris A., Mohanty S. (2009). Sequential updating of multimodal hydrogeologic parameter fields using localization and clustering techniques. Water Resources Reserch 45, 15 PP. Thulin K., Skaug H., Aanonsen S., Naevdal G. (2012). Dual Ensemble Kalman Filters. Submited to: Computers y Geoscience. January 2012. Zhou et al. (2011). Handling nongaussian distributions with Ensemble Kalman Filter. Advances in Water Resources. In press, doi:10.1016/j.advwatres.2011.04.014. Zhou et al. (2011). Pattern Recognition in a bimodal Aquifer using the Normal-Score Ensemble kalman Filter. Mathematical Geosciences. Under review.