Story Transcript
Modelación Matemática y Computacional Aplicada al Diseño Óptimo de Redes de Monitoreo del Agua Subterránea Graciela Herrera Zamarrón Instituto de Geofísica, UNAM
25 de septiembre de 2009
Créditos ► George Pinder, University of Vermont ► Joe Guarnaccia, Ciba Geigy ► Yingqi Zhang, Lawrence Berkeley National
Laboratory ► Edgar Yuri Mendoza, Instituto Mexicano de tecnología del Agua ► Jessica Briseño Ruiz ► Hugo Júnez Ferreira ► Roel Simuta Champo Estudiantes de doctorado Facultad de Ingeniería, UNAM
Guía de la presentación ► ► ► ► ► ►
Redes de monitoreo del agua subterránea Diseño óptimo de redes de monitoreo Diseño de redes de monitoreo con modelos estocásticos Caso diseño de red de calidad del agua para el acuífero Toms River Redes de monitoreo basadas en geoestadística Desarrollos en proceso para diseño de redes de monitoreo basados en modelos estocásticos
Redes de monitoreo de prevención contra la sobreexplotación y contaminación del agua subterránea ► En trabajos de saneamiento de acuíferos sistema para verificar la efectividad de los métodos ► Monitoreo de ► Sistemas
los niveles del agua (cantidad del agua) la calidad del agua
Diseño de redes de monitoreo del agua subterránea ► Consiste
en la elección de sitios y frecuencia de monitoreo ► Métodos geohidrológicos
Sitios en los que se diseña una red por primera vez Sitios en los que aún no se han presentado problemas de contaminación
► Métodos geoestadísticos Sitios en los que se espera que las condiciones cambien en forma paulatina ► Modelos estocásticos de flujo y transporte Sitios en los que se espera que las condiciones cambien en forma paulatina o en los que las condiciones cambiarán abruptamente y es necesaria la predicción
Diseño óptimo de redes de monitoreo requiere optimizar algún criterio para elegir las posiciones de los pozos y/o las frecuencias de monitoreo ► Minimizar ► Se
incertidumbre de estimación costo
Pasos en el diseño de una red de monitoreo óptima ►
Estudio hidrogeológico del acuífero
►
Estudio hidrogeoquímico del acuífero
Modelo conceptual del acuífero
Entendimiento de la problemática de calidad del agua del acuífero
Revisión en campo de los pozos que pueden conformar la red de monitoreo ► Establecer objetivos de la red de monitoreo ► Optimización de la red de monitoreo ►
Traducir objetivos en términos matemáticos Diseño de la red de monitoreo
►
Evaluación de la red de monitoreo
►
Seguimiento, evaluación y modificación de la red de monitoreo
Criterios hidrogeológicos Criterios de calidad del agua Validación estadística
Método para la optimización de redes de monitoreo de la calidad del agua subterránea (1998) Programa de cómputo GWQ-Monitor
Método de Herrera y Pinder Desarrollado en mi tesis doctoral Propósito: Diseño óptimo espacio-temporal de redes de monitoreo de calidad del agua subterránea
Diseño óptimo de redes de monitoreo de calidad del agua subterránea Lenguaje Fortran Cuenta con una interfaz ArgusONE desde 2001
Métodos: Modelación – Modelo estocástico de flujo y transporte Estimación - Filtro de Kalman Optimización sucesivas
Implementa método propuesto en mi tesis doctoral
-
Método
de
adiciones
gráfica
en
Propósito de las redes de monitoreo ► Estimar
el valor de la variable de interés y su evolución en los pozos de monitoreo ► Con base en datos en estos pozos, estimar la variable en puntos y tiempos en dónde no se ha medido ► El monitoreo de las variables se requiere por un tiempo largo
M étodos de dise ño Métodos diseño ► Estimación ► Variables
Espacio Espacio
de decisión
Tiempo Tiempo
Espacio Espacio yy tiempo tiempo Tesis Tesis doctoral doctoral (1998) (1998)
Método de Herrera y Pinder Tres procedimientos 1) Predice la incertidumbre de la estimación de la variable cuando se muestrea en posiciones dadas
2) Utiliza esta predicción de la incertidumbre como criterio para escoger las posiciones y en su caso los tiempos de muestreo, los cuales definen la red de monitoreo y su programa de muestreo
3) Calcula una estimación de la variable a analizar y la actualiza con datos disponibles
Predicción de la incertidumbre cuando se muestrea en posiciones dadas (Método de estimación)
Métodos utilizados estocástico de flujo y transporte ► Análisis geoestadístico ► Filtro de Kalman ► Modelo
Ecuaciones para modelar el flujo del agua subterránea ►
Ecuación en 3D
∂h ⎞ ∂ ⎛ ∂h ∂ ⎛ ∂h ⎞ ∂ ⎛ ∂h ⎞ ⎟⎟ + ⎜ K z −R ⎟ + ⎜⎜ K y ⎟ = Ss ⎜Kx ∂y ⎠ ∂z ⎝ ∂z ⎠ ∂t ∂x ⎝ ∂x ⎠ ∂y ⎝ ►
h -carga hidráulica
►
K x -conductividad hidráulica (capacidad del medio para
►
Ss
►
conducir agua)
-almacenamiento específico (depende de la elasticidad del medio)
R -fuentes o sumideros (extracción de agua por pozos)
Ecuación de transporte de solutos con advección y dispersión ⎞ ∂ ⎛ ∂c ∂c ∂ ⎛ ∂c ⎞ ∂ ⎛ ⎞ ∂c ⎜ ⎟ − V y c ⎟ + ⎜ Dz − Vx c ⎟ + ⎜ D y − Vz c ⎟ = ⎜ Dx ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂t ⎠ ∂z ⎝ ∂z ∂h ∂h 1 ∂h V = − ( K x , Ky , K z ) θ ∂x ∂y ∂z - concentración del soluto
►
c
►
Vx - velocidad efectiva
►
Dx - dispersión hidrodinámica (depende características medio poroso)
►
θ
- porosidad
Resolución de las ecuaciones ► Solución
analítica para problemas simples ► Métodos numéricos Diferencias finitas Elemento finito
Modelo estocástico de flujo y transporte ► Parámetros
aleatorios
Conductividad hidráulica Concentración de la variable en la fuente del contaminante ► Resolución
de las ecuaciones
Simulación estocástica o Monte Carlo
Simulación estocástica realizaciones de los parámetros aleatorios ► Resolver el modelo de flujo y de transporte para cada realización ► Calcular la distribución probabilística o los estadígrafos relevantes de las concentraciones ► Obtener
Filtro de Kalman una estimación lineal sin sesgo y de varianza mínima del estado de un sistema utilizando datos con ruido
► Proporciona
un método secuencial para actualizar las estimaciones cuando se proporcionan datos nuevos, sin necesidad de hacer referencia a datos anteriores
► Establece
Fórmulas del filtro de Kalman ► Variable
a estimar C , datos z ► Ecuación de muestreo
z n = H n C + vn ► Ecuaciones
{(
(
Cˆ n +1 = Cˆ n + K n +1 z n +1 − H n +1Cˆ n
)
Pn +1 = Pn − K n +1 H n +1 Pn K n +1 = Pn H
(H
)(
Pk = E C − Cˆ k C − Cˆ k
para la actualización
T n +1
Cˆ k = E {C / z1 ,..., z k }
PH
n +1 n
T n +1
+ Rn +1
)
−1
{
Rk = E v k v k
T
)} T
}
Requiere de una estimación inicial para C y la matriz de covarianza del error de esta estimación
Cálculo de la matriz de covarianza inicial ► En
el método de Herrera y Pinder la matriz se calcula para todas las posiciones y tiempos de estimación o monitoreo
► Modelo
estocástico de flujo y transporte
Parámetros aleatorios: conductividad hidráulica y contaminación en la fuente de contaminante Caracterizados con un análisis geoestadístico ► Resolución
de las ecuaciones
Simulación estocástica o Monte Carlo
Puntos de estimación y de monitoreo T1
X4
X4
Posibles pozos de monitoreo Puntos de estimación
T2
T3
Matriz de covarianza espacio-temporal Matriz de covarianza espacial
5X5 5X5
C1
C1-C2
C1-C3
X1
σ12
X2
C2-C1
C2
C2-C3
σ22
X3
σ32
X4
σ42 σ52
X5
C3-C1
C3-C2
C3
X1
X2
X3
X4
X5
Selección de las posiciones y los tiempos de muestreo
Métodos utilizados ► Método
de estimación del paso 1 ► Optimización con método de adiciones sucesivas
Función objetivo ►
Varianza en los puntos de estimación sumada sobre todas las posiciones y tiempos de estimación
X4
X4
Posibles pozos de monitoreo Puntos de estimación
Matriz de covarianza del filtro de Kalman
Método de optimización de adiciones sucesivas ► La
estrategia de monitoreo se construye secuencialmente ► Se elige un sitio espacio-temporal de monitoreo a la vez, el sitio seleccionado en cada paso es aquél que minimiza la variancia total ► La búsqueda se realiza por inspección completa de todos los posibles pozos ► Hay que definir un criterio para determinar el número total de pozos
Cálculo de la estimación de la variable
Métodos estocástico ► Filtro de Kalman ► Modelo
Estimación toma la media espacio-temporal de la variable calculada con la simulación estocástica ► Se actualiza con datos usando el Filtro de Kalman ► Se
Resultados
Análisis Geoestadístico
Protocolo 1.1.1.- Aná Análisis geoestadí geoestadístico para obtener mediante simulació simulación secuencial gaussiana realizaciones de conductividad hidrá hidráulica
Configuración en ARGUS-ONE del modelo determinista
Simulación Secuencial Gaussiana SGSIM
2.2.2.- En la interfaz grá gráfica del Princeton Transport Code (PTC), se configura el modelo determinista de flujo y/o transporte
Configuración de la red de monitoreo en ARGUSONE
3.3.- En la interfaz grá gráfica del GWQMonitor se configura el diseñ diseño de la red de monitoreo, monitoreo, se corre el modelo estocá estocástico con las realizaciones de conductividad, se genera la matriz de covarianza en espacio y tiempo, y se se ponderan las posiciones de monitoreo propuestas
+
3.-
4.-
Postprocesamiento y validación
Redes de Monitoreo Piezométricas y/o de Calidad del Agua
4.4.- Se realiza el postpost-procesamiento, procesamiento, se propone el número total de pozos a incluir en la red de monitoreo y se valida la red propuesta
POZOS DE MONITOREO
5.-
CONTAMINANTE SUBSUELO
5.5.- Diseñ Diseño de la red óptima de monitoreo
AGUA SUBTERRANEA
PLUMA
Aplicación para probar el método ► Sitio
del superfondo de la EPA de los EUA
Localizado en el estado de Nueva Jersey
► Contaminante
Clorobenceno (CB)
► Un
sistema de saneamiento funciona en el sitio desde 1985 ► Propósito optimizar la red para un periodo de once meses Seleccionar posiciones y tiempos de monitoreo Minimizar la incertidumbre en dos zonas de riesgo en octubre de 1986
Área de estudio Parque ecuestre Parque Pine Lake
N
Parque Winding River
Parque Industrial
Límite de la Propiedad
Parque Winding River Comunidades de Retirados Ruta 37 Residencial Manchester Residencial Manchester propiedad de Ciba Residencial/Conservation
Area Oak Ridge
Residencial/Conservation propiedad de Ciba Residencial R-90 y R-150 Comercios de la Carretera Rural Industrial 1” =2000 ft.
Industrial Propiedad de Ciba 0
500 1000 1500 2000
Modelo conceptual ► Nueve
miembros geológicos Arenas no consolidadas, limos, y arcillas
► Contaminantes
detectados en las primeras cinco unidades geológicas
redirección Grosor aproximado [pies]
fuente de agua de precipitación escurrimiento
zona vadosa
Pozo de bombeo
fuente de agua de la zona vadosa agua colgada 50
infiltración
arcilla fuente de agua de la zona colgada capas del modelo 1 Cohansey Primario
30
2 3
30
4
10
5
Cohansey Inferior
30
6
Kirkwood Superior
10
7
Kirkwood No. 1
zona saturada (Acuífero de Arena Superior)
Transición Cohansey/Kirkwood
Río Toms
Dominio y malla del modelo
Posibles pozos de monitoreo
62 pozos
Red de monitoreo y número de muestras Red óptima 23 pozos 40 muestras 1
Red original 322 muestras de 62 pozos
Reducción en costo del 87%
2
1
2
1
1
1
2 3
1 1
11 7 11 2 1
2 3
1 3 1
Número de muestras por tiempo 7
6
Número de muestras
5
4
3
2
1
0 1 2 3
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 Tiempo de muestreo
Necesidades tiempo de cómputo ► Método de estimación de parámetros para el modelo estocástico ► Incluir costo en la función objetivo ► Reducir
Desarrollos Posteriores
Método geoestadístico ► Tesis
maestría Hugo Júnez
Adaptación del método a problemas sin modelo ►Matriz
de covarianza de error de la estimación espacial obtenida por análisis geoestadístico
Diseño de redes de monitoreo para varios parámetros ► Aplicaciones
Acuíferos: Irapuato Valle (2004), Querétaro (2004), Pátzcuaro (2006, 2007), Chalco-Amacameca, Zona Metropolitana de la Ciudad de México y Texcoco (2008), Acuíferos somero y profundo de San Luis Potosí (2006-2009) Otros autores 5 acuíferos de Chihuahua
Desarrollos actuales método geoestadístico ► Tesis
doctoral Edgar Yuri Mendoza (terminada) Método geoestadístico espacio-temporal
► Tesis
doctoral Hugo Júnez
Diseño óptimo redes de monitoreo espacio temporales
Desarrollos modelación estocástica ► Tesis
Yingqi Zhang (dirigida por el Dr. Pinder) Método Latin Hypercube Sampling en 3D (para reducir tiempo de cómputo) Optimización algoritmo genético con posibilidades de incluir costo Función objetivo para estimar las tendencias de concentraciones en algunos pozos
Desarrollos actuales modelación estocástica ► Tesis
doctoral Jessica Briseño
Modificación para diseño de redes de la cantidad del agua Método de estimación de parámetros para el modelo estocástico
Método de estimación de parámetros para el modelo estocástico ► Filtro
de Kalman ► Aplicado a la matriz de covarianza de la concentración aumentada con la conductividad hidráulica y la carga hidráulica Matriz Matriz de de covarianza covarianza H-lnK-C H-lnK-C H
H-ln K
H-C
ln K-H
ln K
ln K-C
C-H
C-ln K
C
Matriz de covarianza de H-lnK-C
Método propuesto para la estimación de parámetros
Con datos de K
Análisis geoestadístico
ln K, H, C, matriz de covarianza de H-lnK-C como parámetros a priori
Datos de ln K, H y C
ln K y su semivariograma
Realizaciones de ln K
Modelo estocástico de flujo y transporte
Filtro de Kalman
Realizaciones de H y C
Cálculo de H, ln K, C, y matriz de covarianza de H-lnK-C en espacio y tiempo
ln K actualizada
Caso de estudio Estimación de parámetros de un modelo estocástico de flujo y transporte. Este Este caso caso de de estudio estudio se se basa basa en en una una representación representación simplificada simplificada del del acuífero acuífero de de Querétaro. Querétaro.
Se considera un derrame en la zona centro del Valle de Querétaro. Suponemos que el derrame se tiene una fuente constante de emanación de 1 y se dispersa durante 50 años en los cuales no tenemos datos. Modelo completo
Posteriormente, consideramos que el segundo periodo se cuenta con datos de ln K, H, y C en algunos puntos cada 6 meses.
8073 nodos 15860 elementos
Modelo reducido
Condiciones Condiciones de de frontera frontera Se establecieron como condiciones de frontera de carga asignada a la media de las 4000 realizaciones de H
Matriz Matriz de de covarianza covarianza
Convergencia de las realizaciones de H, ln K y C Datos de convergencia de las matrices de covarianza de H, ln K y C
Modelo completo Gráficas de convergencia de las matrices de covarianza de H, C y ln K
Datos de convergencia de las matrices de covarianza de H, ln K y C
Modelo reducido Gráficas de convergencia de las matrices de covarianza de H, C y ln K
Comparación de los modelos (completo vs reducido) Comparación de la media de H del modelo Completo vs. Reducido Caso
EM
ECM
Min Err
Max Err
Modelos C vs R
0.00144
0.00015
-0.02420
0.02130
Comparación de la estimación de H del modelo completo vs. reducido en las estimaciones con el filtro de Kalman en los casos de estudio.
Estimación de H
EM
ECM
Min Err
Max Err
Caso A
0.11355
0.05459
-0.33000
0.49300
Caso B
0.12600
0.02464
-0.05900
0.29400
Caso C
-0.09742
0.01511
-0.20700
0.08900
Avances
Si ln K y Cov ln K no son correctas
LOCALIZACIÓN GENERAL 112 °
104 °
9 6°
88°
E S T A D O S U N I D OS D E A M É R I C A
32°
B A JA C A L I FO R N I A NTE
32 °
No rte
S O N OR A CH IH UA H U A
C O A H U I LA B AJA C A L I FO R N I A S UR 2 4°
S I NA LO A
N U E V O LE Ó N G O L FO D E M É X I C O
D U RA N GO
24°
TA M A U L I P A S
ZA C A T E C A S
S LP NAY A RIT
A GS Y U C A TA N
N ort e
Estimación Estimación de de parámetros parámetros de de un un modelo modelo estocástico estocástico de de flujo. flujo.
Si Si ln ln K K yy Cov Cov ln ln K K son son correctas, correctas, empleando empleando el el filtro filtro de de Kalman, Kalman, se se estimó estimó el el ln ln K, K, con con datos datos de de ln ln K K yy H. H. El El campo campo distribuido distribuido de de ln ln K K estimado estimado tiene tiene poco poco error error aunque aunque el el efecto efecto de de los los datos datos de de H H en en la la estimación estimación es es reducido. reducido.
J A LI S C O
CO LI MA
GT O
QRO H GO
Q U I N TA N A R O O
MÉ X D . F. M OR P UE B LA GU ER R E RO
CA MP E C H E
MI C H
TA B A S C O B E L I ZE
1 6°
OA X A C A
C H IA P A S
1 6 °
G U A TE M A L A HON D UR A S
San Luis Potosi
E L S A LV A D O R
112°
1 04 °
96°
88°
Area de Estudio
Guanajuato
Queretaro Qto.
Hidalgo
Michoacan
Edo. de México
Se requiere estimar K y Cov K, con las cuales al realizar la estimación de parámetros se logré menor error, y que hagan que el FK de resultados apropiados.
Desarrollos actuales modelación estocástica ► Tesis
doctoral Roel Simuta
Método Latin Hypercube Sampling en 3D Optimización algoritmo genético y comparaciones con método secuencial Extensión al muestreo en 3D
Error en la covarianza del ln K (2D) 0.12
0.1
LHS_axy=2500 LHS_axy=672 SGS_axy=2500 SGS_axy=672
0.06
0.04
0.02
Número de realizaciones
1600
1500
1400
1300
1200
1100
1000
900
800
700
600
500
400
300
200
100
0 0
RECM
0.08
Error en la covarianza del ln K (3D) 0.12
0.1
RECM
0.08
LHS_axyz=2500 LHS_axyz=672 SGS_axyz_2500 SGS_axyz=672
0.06
0.04
0.02
0 0
200
400
600
800
1000
1200
Número de realizaciones
1400
1600
1800
2000
Desarrollos futuros ► Paralelización
del programa de cómputo ► Unificar los métodos geoestadístico y de modelación