Modelos matemáticos hidrogeológicos en diferencias finitas (MDF) y elementos finitos (MEF)

Definición

El método de las diferencias finitas (MDF) ha sido ampliamente usado en los estudios de aguas subterráneas desde tempranamente los años 60. El MDF fue estudiado por Newton, Gauss, Bessel y Laplace (Pinder y Gray 1970).

La exactitud del método depende del tamaño de la red y la uniformidad. La aproximación de la derivada mejora a medida que el espaciamiento en la malla se acerca a cero; sin embargo, la dispersión numérica y el error de truncamiento se incrementa. Hay tres métodos diferentes de aproximación de diferencias finitas: adelante, atrás y la diferencia central, según la forma en que se aplica las diferencias finitas. La diferencia central da mejores resultados dado que el error de truncamiento es de segundo orden o (Δx)2.

La ecuación general que rige para las condiciones transitorias, heterogéneas y anisotrópicas está dada por:

Para simplificar, se consideró un caso unidimensional y se resolvió para h mediante el método de diferencias finitas, dando como resultado:

Las ventajas del método de diferencias finitas son que este es fácil de aplicar, bien documentado y produce unos resultados bastante buenos. Sin embargo, el método de diferencias finitas tiene algunas desventajas. La principal desventaja es que este no se ajusta adecuadamente a una frontera irregular del modelo (Fig. 6). Además, la distribución de las redes, su tamaño, y si ellas son de igual tamaño, afectan elevadamente la exactitud y el esfuerzo de cálculo.

Las más ampliamente usadas diferencias finitas, basadas en modelos de aguas subterráneas es el MODFLOW (Harbaugh y McDonald 1996).

De otro lado, la base del método de elementos finitos es la solución de ecuaciones integrales sobre el dominio del modelo. Cuando el método de elementos finitos se sustituye en las ecuaciones de diferenciales parciales, se produce un error residual. El método de elementos finitos fuerza este residuo hasta ir a cero (Konig L. y Weiss J., 2009).

El método de elementos finitos discretiza el dominio del modelo en elementos. Estos elementos pueden ser triangulares, rectangulares o bloques prismáticos. El diseño de la malla es muy importante en el método de elementos finitos ya que afecta significativamente a la convergencia y la exactitud de la solución. El diseño de la malla en el método de los elementos finitos es un arte más que una ciencia, pero hay reglas generales para una buena configuración. Esto es muy recomendable para asignar los nodos en los puntos importantes, como una fuente o sumidero, y para refinar la malla en áreas de interés donde las variables cambian rápidamente.

El método de los residuos ponderados está siendo ampliamente utilizado en los problemas de aguas subterráneas de elementos finitos. El ingeniero ruso B. G. Galerkin introdujo este método en 1915 (Pinder y Gray 1970). Para ilustrar la aproximación de los residuos ponderados, se considera un problema de agua subterránea o transporte de soluto. Este problema en un dominio B se puede simplificar como:

Para solucionar la ecuación, la función de ponderación W (x,y,z) necesita ser identificada.

Las ventajas de este método incluyen: una mejor configuración de la malla, lo cual se ajusta a las fronteras irregulares del modelo, la anisotropía es bien incorporada, el sistema de gobierno de las ecuaciones es simétrica y las formas irregulares se pueden utilizar para representar los elementos.

El método de los elementos finitos tiene algunas desventajas. La malla del elemento finito no es fácil de construir y consume mucho tiempo, especialmente en problemas complicados. También, a diferencia del método de diferencias finitas, el balance de masa en el método de elementos finitos se puede alcanzar para todo el dominio, pero no para cada elemento (Konig L. y Weiss J., 2009).

El más conocido elemento finito basado ​​en los modelos de aguas subterráneas son Feflow, Femwater y MODFE (Harbaugh y McDonald 1996).

Investigación

Consultoría  para la actualización del modelo matemático de los acuíferos Rímac y Chillón. Servicio de Agua Potable y Alcantarillado de Lima (SEDAPAL)

Ubicación: Distritos de Miraflores, San Isidro, La Molina, San Miguel, Santa Anita, Ate, Lurigancho, San Martin de Porres, Los Olivos, Ancón, Puente Piedra, Carabayllo y Comas (Lima – Perú)

Responsable: Equipo Multidisciplinario

Duración: 1 año y 4 meses (2010)

Los acuíferos del Rímac y Chillón constituyen una de las fuentes de agua de importancia en la historia de la ciudad de la gran Lima. El crecimiento poblacional de la ciudad capital se halla en incremento, consecuentemente las demandas de agua con fines de abastecimiento poblacional, tienen el mismo ritmo. Frente ante esta realidad, SEDAPAL viene implementando diversos proyectos en la cuenca alta del rio Rímac a fin de almacenar y aprovechar en forma racional el recurso hídrico superficial conjuntamente con la explotación de las reservas del acuífero de Rímac.   En ese mismo contexto, en el valle Chillón, la empresa Agua Azul, viene haciendo uso conjunto de las aguas subterráneas y superficiales, para ello ha construido pantallas de recarga con 100 ha de superficie, a través de las cuales el acuífero renueva sus reservas durante los meses Diciembre a  Abril, y explotar agua subterránea mediante 28 pozos tubulares durante los meses de Mayo a Noviembre, concretando el abastecimiento de agua en forma sostenida al cono norte de la ciudad de Lima.

La información hidrogeológico posterior a 1997,  ha sido revisado en 259 estudios orientados a la perforación de nuevos pozos, en muchos de ellos hacen uso de información ya existente y en otros registran pruebas hidrodinámicas orientadas a conocer la conductividad hidráulica, mas no así el rendimiento especifico, por adolecer de piezómetros cercanos a los pozos. Así mismo, se ha revisado información de registros de 103 SEVs, georeferenciados, en otros casos, han sido tomados de información existente para sus propósitos.

Conjuntamente, con la información recopilada y los nuevos escenarios de bombeo de las aguas subterráneas, se elaborara el modelo matemático de los acuíferos Rímac y Chillón, la misma que será calibrada y validada, a fin de simular escenarios que permitan tener los elementos de juicio y tomar las decisiones pertinentes relacionadas con el manejo racional de las reservas de las aguas subterráneas.

Recomendación

Simular el flujo de las aguas subterráneas usando el programa Visual Modflow, en régimen estacionario y no estacionario.

El desarrollo del modelo conceptual es uno los pasos de importancia en el desarrollo de los modelos de simulación de acuíferos. Por definición, el modelo conceptual es la representación simplificada de las características hidrogeológicas del sistema acuífero, así como los detalles del comportamiento hidrológico del acuífero y sus externalidades.

El modelo conceptual del acuífero ha sido elaborado acorde con la información geológica, límites permeables e impermeables, marco hidrológico, sus variaciones naturales, y los factores antropogénicos preponderantes.

Selección del código de modelamiento

El código de modelamiento, es el programa de cómputo que contiene los  algoritmos para resolver numéricamente las ecuaciones de flujo que rige en medios porosos saturados.

El modelo conceptual y el esquema hidrogeológico ayudaran a definir las condiciones de frontera para la solución del modelo matemático para el acuífero de Rímac y Chillón. Con este propósito, se ha optado por hacer uso del software Visual ModFlow Premium 4.2

Diseño del modelo

La discretización del sistema acuífero en los modelos que anteceden a la presente fueron realizados con mallas de 1000 m x 1000 m y en otro caso con mallas de 100 m x 100 m, en vista de esto se ha uniformizado las mallas en diferencias finitas, empleando mallas cuadradas Δx = 100 m  y  Δy = 100 m, habiéndose inicializado la discretización del sistema con 450 filas y 443 columnas, haciendo un total de 199350 mallas, de ellas 57272 mallas corresponden al acuífero propiamente cubriendo una  superficie estimada en 572.72 Km2.

Las capas son empleadas en los modelos para representar las unidades hidroestratigráficas, las cuales son unidades geológicas con similares o diferentes propiedades. Acorde con los perfiles estratigráficos se conoce que en el acuífero de Rímac y Chillón es monocapa, consecuentemente el modelo de acuífero es monocapa, con espesores y cotas definidas, tanto en el basamento y el techo, definidas en base a la topografía establecidas y las cotas de los SEV´s, este último principalmente para conocer la topografía subterránea.

El tipo de acuífero (Layer Tipe 2) considerado acorde con las opciones del modelo corresponde a acuífero libre/confinado con Sy (rendimiento especifico) y T (transmisibilidad) variable.

La información analizada ha establecido 15 rangos de valores de conductividades hidráulicas.

Los valores se hallan comprendidos entre 4 m/d y 172.8 m/d, caracterizando a un acuífero muy heterogéneo.

De otro lado, la primera condición de frontera está referida a las mallas con potencial impuesto correspondiente al contacto hidráulico entre el acuífero y el Océano Pacifico. La carga hidráulica asignada a estas mallas es  h= 0 msnm.

La segunda condición de frontera o frontera tipo 2 (Tipo Neumann), denominado también frontera con flujo especificado, corresponde básicamente a las zonas irrigadas a través de las cuales se produce la recarga del acuífero.

Acorde con los cálculos efectuados por Cruz (2009), la recarga subterránea a través de la sección de flujo en la zona Huampani, es estimado en 143317 m³/día (1.66 m³/s).

De las evaluaciones efectuadas de la morfología del agua subterránea en la zona de Punchauca por Rafael (2004), para el periodo Noviembre 2001 hasta Marzo del 2004, a nivel mensual, se tiene que el caudal medio de flujo subterráneo a través de la sección Punchauca, en la inmediaciones del área del desarrollo del proyecto de recarga inducida del acuífero Chillón, el caudal promedio estimado es 6503 m3/día, que hacen un total de 2.30 MMC/año.

El tercer tipo de condición de frontera corresponde al comportamiento del hidrograma de los ríos Rímac y Chillón, así como de sus cargas hidráulicas y del acuífero. En base a la información de las descargas medias mensuales estimadas del registro de caudales (1965-2010) del rio Rímac, y considerando la eficiencia de conducción del rio en 90%, los caudales netos de recarga corresponden al 10% del caudal medio,

Similar procedimiento de cálculo efectuado a las descargas medias mensuales del rio Chillón (1969-2009), la intensidad de recarga del acuífero a nivel mensual, equivale a 15.20 MMC/año.

La explotación del agua subterránea de los acuíferos Rímac y Chillón se realiza mediante pozos tubulares y tajo abierto, desde que el modelo en régimen estacionario será calibrado para Julio del año 1985, se ha elaborado la base de datos de 876 pozos en explotación con -616811.507 m3/d de caudal total, lo que equivale a una descarga constante de 7.13 m3/s.

Para el proceso de calibración del modelo se hará uso de las cargas hidráulicas correspondiente a Julio de 1985 de 63 pozos de observación los cuales son representativos para los acuíferos en estudio.

Calibración del modelo en régimen estacionario

El proceso de calibración del modelo ha consistido en:

  1. Selección de los parámetros de entrada
  2. Simulación de flujo de las aguas subterráneas mediante el software Visual ModFlow Premium 4.2
  3. Comparación entre las cargas observadas y calculadas
  4. Selección de nuevos valores de los parámetros de entrada orientados a minimizar la diferencia entre los valores de las cargas calculadas y observadas.
La comparación de los residuales y la curva de distribución normal, y con evidencia sabemos que los residuales se aproximan a la distribución normal, lo cual garantiza que el modelo se encuentra en condiciones de ser empleadas para procesos de simulación del modelo en régimen de flujo transitorio.

Los resultados encontrados obedecen a la variación de la conductividad hidráulica en sus valores iniciales.

Ha sido necesaria adicionar las zonas con mayor intensidad de recarga principalmente en la zona media del acuífero Rímac, por cuanto en las corridas se generaban mallas secas, habiéndose optado por modificar el balance hídrico realista en los acuíferos. Las intensidades de recarga neta inicial, Ri(mm/año), y final, Rf(mm/año)

Balance hídrico en régimen estacionario

Análisis de sensibilidad

El análisis de sensibilidad se ha sido desarrollado para las conductividades hidráulicas calibradas, el criterio ha consistido en disminuir e incrementar en 50% del valor calibrado. La sensibilidad ha sido cuantificada mediante la variación de la raíz media cuadrática normalizada.

Las zonas con más alta sensibilidad al incremento de los valores de la conductividad hidráulicas corresponde a las zonas 8 y 9 con  ∆RMS-N que varía entre 4.85% y 5.73%

Calibración del modelo en régimen no permanente

El modelo en régimen de flujo transitorio, está orientado a simular el comportamiento de los niveles de agua subterránea en el pasado y el futuro.

Reproducir la fluctuación de los niveles de agua en el espacio y el tiempo, requiere disponer de información histórica de las variables actuantes sobre el acuífero tales como: el régimen de bombeo de los pozos, intensidades de recarga en el espacio y el tiempo, las cargas hidráulicas observadas en forma periódica y finalmente la distribución espacial del coeficiente de almacenamiento y/o rendimiento especifico, todos ellos consistentes.

Con el objeto de simular el flujo del agua subterránea en régimen transitorio, se ha establecido 06 zonas con valores de rendimientos específicos, con valores comprendidos entre 0.01 a 0.20.

Cargas iniciales y fronteras del modelo

Los resultados de las cargas hidráulicas del modelo calibrado en régimen estacionario para Julio de 1985 han sido utilizados como cargas iniciales para el proceso de simulación del modelo en régimen transitorio. Estas cargas constituyen el escenario base o la línea base del proyecto, sobre el cual actuarán las descargas posteriores a esta fecha, así como la información de los escenarios a ser simulados en el futuro.

El proceso de calibración del modelo de flujo del agua subterránea en régimen transitorio consiste en minimizar  la diferencia entre los datos históricos de las cargas hidráulicas con las calculadas por el modelo.

Las diferencias ocurren en todo modelo, las cuales provienen de la interpretación geológica del acuífero, registros históricos (caudales de bombeo, intensidades de recarga), por lo que es permitido variar los valores inicialmente ingresados al modelo, hasta lograr la menor discrepancia, mediante el proceso iterativo de prueba y error.

Se observa que los hidrogramas calculados siguen la tendencia de los observados, sin embargo el análisis por zonas permitirá cuantificar las discrepancias entre ambas cargas hidráulicas.

Los Hidrogramas Calculados y Observados de los piezómetros, permiten en todos los casos se ha logrado simular la tendencia entre ambas cargas hidráulicas, y con discrepancias permisibles, con medias cuadráticas normalizadas inferiores al 2%, por lo que el modelo se encuentra en condiciones de ser utilizadas para simular los escenarios requeridos. 

Escenarios de simulación

La simulación de escenarios es una etapa de importancia por cuanto, se ingresa a la etapa de uso del modelo con el objeto de formular interrogantes y que solo el modelo puede dar respuesta con cierta aproximación. Se reitera que las variables actuantes sobre el sistema acuífero difieren acorde con las características del año o ciclo hidrológico, en virtud a ello, simulaciones de hasta 10 años se consideran aceptables.

Escenario 1

Considera constante los valores de las variables actuantes sobre el sistema acuífero considerados en el año 2009 (caudales de bombeo de los pozos de SEDAPAL, pozos de terceros con NIS, pozos de terceros sin NIS, recarga a través del rio, recarga a través de las áreas verdes), hasta el año 2035.

Escenario 2

Considera la puesta en operación de Ampliación de la Planta Chillón y paralización de la batería de 28 pozos y los 83 pozos localizados en la respectiva  área de influencia.

Escenario 3

El tercer escenario considera la puesta en marcha de proyectos de los esquemas de abastecimiento de agua en la zona de  Ate-Vitarte y paralización de pozos en Proyecto Lima Norte II.

 

Escenario 4

Contempla la puesta en marcha de los  proyectos de  los esquemas de abastecimiento  de agua para Pariachi programado para el año 2012 con una descarga de 0.235 m3/s, de manera similar la puesta en marcha la explotación del agua subterránea para abastecer Santa Anita a partir del 2013 con una descarga de 0.295 m3/s, posteriormente en el 2014 serán implementados loes esquemas de abastecimiento de agua a través de pozos para Cajamarquilla, Carapongo y Nicolás de Piérola. 

Escenario 5

El escenario contempla la puesta en operación del proyecto de la Planta Alto Chosica.  Los pozos mostrados en el Cuadro №3.19, en total 106 serán paralizados  a partir del año 2015 al 2035 y la incorporación de 30 pozos a ser perforados en la margen derecha del rio Rímac con 1m³/s de descarga, los cuales fueron mostrados en el cuadro №3.16, estos pozos iniciaran su operación en el año 2024 hasta el 2035.

Escenario 6

Contempla las acciones tomadas en consideración  en los escenarios 2, 3, 4 y 5, esto es la puesta en marcha la operación de la ampliación de la Planta Chillón y la paralización de los 28 pozos y los del área de influencia, puesta en marcha de los proyectos considerados en los esquemas de abastecimiento de Ate Vitarte y la paralización de los pozos para el Proyecto Lima Norte II, puesta en marcha de los proyectos considerados en los esquemas de abastecimiento de Ate Vitarte y la paralización de los pozos para el Proyecto Lima Norte I y la puesta en operación de los proyectos de la Planta Alto Chosica y la operación de la batería de pozos en el mediano y largo plazo.

Análisis de intrusión marina

En base al modelo formulado para el escenario 6, se ha analizado la problemática de la intrusión marina en el acuífero del Rímac. El análisis se inicia desde el año 1985, pasando por los años críticos del 1997 – 1998, para llegar a los años 2003, 2009 y 2035.