Modelo de administración de acuíferos: MADA

Groundwater management model: MADA

Eric Cabrera-Estupiñán, Félix Dilla-Salvador
Instituto Superior Politécnico José Antonio Echeverría, Cuba

Dr. Eric Cabrera Estupiñán
Dr. Félix Dilla Salvador

Instituto Superior Politécnico "José Antonio Echeverría"
Centro de Investigaciones Hidráulicas
Calle 114 núm. 11901 entre 119 y 127
Teléfonos: +53 (7) 2663 455 y 2601 416
19390 Marianao, Provincia de La Habana, Cuba
ecabrera@cih.cujae.edu.cu
fdilla@cih.cujae.edu.cu

Recibido: 17/04/09
Aprobado: 15/06/11

Resumen

El trabajo propone un modelo matemático para la administración óptima de la explotación de acuíferos llamado MADA. Dicho modelo está basado en la inclusión de un modelo de simulación del flujo de agua subterránea (AQÜIMPE) dentro de un algoritmo de programación lineal. Se plantea la concepción teórica de MADA realizando una especial distinción a la inclusión de AQÜIMPE dentro del modelo propuesto. Finalmente se emplea MADA en la evaluación de la explotación de la "Cuenca Sur de La Habana", haciendo uso de un Plan de Uso del Agua (PUA), donde queda demostrado que la distribución de la explotación empleada para garantizar el PUA no evalúa las posibilidades de intrusión salina en el acuífero y a tal efecto se brinda una solución. Además, se evalúa la posibilidad del acuífero de satisfacer la demanda bajo un patrón anual de lluvia del 75% de probabilidad, en un periodo de cinco años secos consecutivos.

Palabras clave: optimización, modelación del agua subterránea, AQÜIMPE, MADA, programación lineal, Cuenca Sur de La Habana.

Abstract

In this paper a mathematical model (MADA) is proposed for the optimal management of aquifer exploitation. This model is based on the inclusion of a simulation model for groundwater flow (AQÜIMPE) using a linear programming algorithm. The theoretical concept behind MADA is described, making a special distinction for the inclusion of AQÜIMPE in the proposed method. Finally, MADA is used to evaluate the exploitation of the South Havana Basin according to a Water Use Plan (WUP), where it is demonstrated that the distribution of the exploitation used in the WUP does not take into account possible saline intrusion into the aquifer; a solution to this effect is provided. Also evaluated is the ability of the aquifer to satisfy demand based on a yearly pattern of 75% probability of rain during a period of five consecutive years with low rainfall.

Introducción

Desde la década de los años sesenta, los sistemas de agua subterránea han sido estudiados mediante modelos matemáticos de simulación, los cuales se han venido desarrollando de una forma vertiginosa. En los últimos anos, los modelos de simulación han sido combinados con técnicas de optimización para determinar la mejor estrategia a seguir en la administración del recurso hídrico, con un objetivo específico y un conjunto de restricciones (Galloway et al., 2003).

Estas técnicas de optimización pueden ser clasificadas en dos grupos: las técnicas no determinísticas, dentro de la que se pueden mencionar a las redes neuronales y los algoritmos evolutivos fundamentalmente; y las técnicas determinísticas, dentro de las que se encuentran, entre otras, la programación lineal, la cuadrática y la programación dinámica (Marrero, 1985). Para resolver estos problemas de optimización, muchas veces se aplica de forma aislada alguno de estos algoritmos, sin embargo, cada vez es más frecuente que se realicen combinaciones entre ellos para evitar obtener una solución no deseada o para disminuir el tiempo de cómputo. Por ejemplo, Hsiao y Chang (2002) proponen un método novedoso para encontrar la distribución de pozos de bombeo, y el caudal a extraer en cada pozo y en cada intervalo de tiempo, que garantice un mínimo de costo global del sistema, donde se incluyen los costos de operación o bombeo, y los costos fijos o de emplazamiento. Este método combina un algoritmo genético para la determinación de la cantidad y ubicación de los pozos de bombeo, y otro algoritmo llamado Programación Dinámica Diferencial con Restricciones (CDDP) para la determinación de los costos dinámicos o de operación.

Otros ejemplos de aplicaciones realizadas en el mundo en la administración de los recursos hidráulicos pueden ser consultados en Gutiérrez (2008).

En Cuba se produjo un avance muy importante en esta rama del conocimiento con la creación del modelo matemático de simulación del flujo subterráneo AQÜIMPE (Martínez, 1989), en el Centro de Investigaciones Hidráulicas (CIH). Numerosas han sido las aplicaciones de este modelo en acuíferos nacionales e incluso en el extranjero (Gutiérrez, 2008).

Sin embargo, sólo con una herramienta de simulación no se pueden resolver con rigor científico un conjunto de problemas referentes a la administración de la explotación.

En Cuba, una de las primeras aplicaciones de patrones óptimos de explotación prolongada en un acuífero fue planteado por Dilla (1993). En este trabajo se emplea el método de la matriz respuesta o funciones tecnológicas algebraicas (FTA), donde finalmente se utiliza un método de optimización, la Programación Lineal, y de esta forma se resuelven varios problemas de optimización de naturaleza lineal.

Posteriormente, Mesa (2004) da solución al problema inverso de estimación de los parámetros hidrogeológicos, conductividad hidráulica (K) y coeficiente de almacenamiento (S), aplicando un algoritmo genético de convergencia global llamado Shuffled Complex Evolution (SCE) y fundamentado en el modelo AQÜIMPE. En este trabajo también se recomienda emplear el método SCE en la resolución de problemas referentes a la administración de acuíferos.

En la presente contribución se discute el modelo de administración MADA, basado en la inclusión del modelo AQÜIMPE dentro de un algoritmo de Programación Lineal, y se aplica al acuífero de la Cuenca Sur de La Habana.

Fundamentos y estructura general del modelo AQÜIMPE

El modelo matemático AQÜIMPE en esencia resuelve numéricamente la ecuación diferencial parcial de segundo orden, de tipo parabólica, que surge al analizar el flujo bidimensional, impermanente, lineal, y libre o confinado en medio poroso, con término fuente-sumidero y condiciones de contorno filtrante o impermeable (Martínez, 1989):

Donde:

h: carga hidráulica [L].

b: espesor saturado del acuífero [L].

K: conductividad hidráulica [LT^-1].

W: fuente de entrada o salida vertical al elemento, es decir, explotación o recarga.

La solución de esta ecuación puede expresarse en la obtención de la carga hidráulica (o piezométrica), como variable dependiente de las coordenadas del espacio bidimensional (x, y) y del tiempo (t).

AQÜIMPE resuelve la ecuación (1) mediante el método de los elementos finitos (MEF), utilizando el triángulo cuadrático como elemento, y aplicando la aproximación de Galerkin (Martínez et al, 2000).

Con AQÜIMPE se pueden simular acuíferos confinados y libres a escala regional, lo que incluye el tratamiento del acuífero visto en planta, aunque también se pueden resolver problemas vistos en una sección vertical. Además permite hacerlo a escala local, para simular, por ejemplo, pozos de bombeo y campos de pozos. Asimismo, es posible simular su interacción con un cuerpo de agua superficial, siempre que exista comunicación entre ellos.

La estrategia utilizada para la modelación del acuífero es subdividir el área total de la formación acuífera en distintas zonas, atendiendo a las propiedades hidrogeológicas; posteriormente estas zonas son subdivididas en elementos triangulares, donde en cada uno de éstos las propiedades hidrogeológicas son iguales. Al emplearse el triángulo cuadrático para la discretización espacial, la solución en cada elemento se obtiene en seis nodos y por lo tanto ésta puede ajustarse a una superficie cuadrática, aspecto que permite representar con menor número de elementos en la malla una región, logrando la misma precisión de los resultados que si se emplearan triángulos lineales (Martínez, 2006). Además, con el triángulo cuadrático se puede modelar mejor el medio cársico como elemento volumétrico representativo de un medio poroso equivalente. En este sentido, la estrategia de modelación es definir la hipótesis de sustitución del medio cársico real por un medio poroso equivalente, cuyas propiedades de comportamiento se determinan definitivamente durante la calibración del modelo (Cabrera, 2007).

No es objetivo plantear todo el desarrollo del MEF aplicado a la ecuación (1) (Martínez, 1989). Por eso, se ilustra el sistema matricial expuesto en la ecuación (2), obtenido después de la resolución de las integrales que se generan con la aproximación de Galerkin a la ecuación (1):

Donde:

se conoce como matriz de rigidez del elemento.

K_e: conductividad hidráulica del elemento e [LT^-1].

[K_ab]^e_6x6 matriz simétrica cuadrada de orden 6, cuyos valores son los mismos para todos los tiempos de un elemento e, pero diferentes para todos los elementos.

matriz constante.

Δt: cantidad de unidades de tiempo en el intervalo de tiempo.

[K_cte]_6x6: matriz constante para todos los elementos en todos los tiempos.

{h_e}^k_6x1: valor de las h de los nodos del triángulo o elemento e; cada término del vector está referido a uno de los seis nodos del triángulo y depende del tiempo k [L].

{h_e}^k-1_6x1: igual al anterior, pero en el tiempo previo [L].

Qs^k_e: caudal superficial, en forma de intensidad que entra al elemento e, recarga debida a la lluvia en el tiempo k [LT¹].

Qp^k_j: caudal puntual que entra o sale al nodo j de explotación en el tiempo k [L³T^-1].

L_i: longitud del lado i [L].

{000111}_6x1: vector columna, los tres primeros valores son los referentes a los nodos principales y el resto a los nodos secundarios.

{0...100}_nx1: vector columna que tiene tantos elementos como nodos de la malla y en él tiene valor uno los referentes a los nodos donde hay gasto puntual.

{114}_3x1: vector columna de tres elementos; garantiza que en los lados donde hay gasto lateral se reparte 1/6 de este caudal a cada nodo principal y los restantes 2/3 al nodo secundario.

n: número total de nodos de la malla.

nl: número de lados de triángulos con salidas o aportes laterales (zonas del acuífero que tienen frontera con un río, etcétera).

np: número de nodos con gastos puntuales.

ns: número de triángulos con entradas o salidas superficiales.

ne: número total de elementos (triángulos).

i: contador de lados con gasto lateral.

e: contador de elementos.

ω: tiene un valor de 2/3 si se obtiene por el MEF.

Cuando se aplica esta expresión a un acuífero confinado, el problema se restringe a resolver k veces el sistema matricial dado en la ecuación (2), el cual se resume en el sistema de ecuaciones lineales mostrado en la ecuación (3):

Donde [M] es una matriz cuadrada de orden n, en la que cada fila está asociada con uno de los nodos de la malla. Lo mismo sucede con las columnas. En el vector h se guardan las variables a calcular en cada tiempo y f es el vector de los términos independientes donde se almacena el efecto de los caudales extraídos e incorporados al modelo en cada nodo.

Estructura del modelo de programación lineal

El asistente matemático MATLAB (García et al., 2005) posee una función llamada Linprog, que resuelve el problema de programación lineal mediante el método Simplex. Esta función es usada en el modelo propuesto y en su forma más general presenta la siguiente estructura:

Donde c, A y r son las matrices que intervienen directamente en el problema, x es el vector de variables de decisión que se obtiene y c es un vector de costo o de ponderación de cada una de las variables de decisión. Con la ecuación (4) se resuelve el siguiente problema:

Sujeto a restricciones de desigualdad:

Restricciones de igualdad:

Con fronteras inferiores y superiores para cada variable de decisión:

El vector x_o representa una solución inicial del problema y en el término feval se almacena el valor evaluado de la función objetivo.

Funciones objetivo y restricciones del modelo de administración

En el modelo MADA se pueden tomar en cuenta un conjunto de funciones objetivo (Dilla, 1993).

Antes de analizar las funciones objetivo se define el término npa como la cantidad de nodos de la malla numérica con gasto puntual a optimizar.

a) Maximizar la suma de las cargas en los nodos con gasto puntual a optimizar del modelo en todos los tiempos:

b) Minimizar la suma de los abatimientos s^T_j calculados por la diferencia entre los niveles en el estado inicial y al final del horizonte de planeamiento o tiempo final T en todos los

Por otro lado, la tercera función objetivo concibe la máxima extracción posible.

c) La mayor capacidad de extracción de un acuífero se obtiene al maximizar los caudales en los nodos de extracción a optimizar npa para todos los tiempos:

En general, las restricciones pueden resumirse en tres tipos:

1. La suma de los caudales extraídos en cada periodo por todos los nodos de bombeo del modelo a optimizar npa debe igualar o superar una demanda prefijada en el periodo. Este conjunto de restricciones genera T ecuaciones como la (12):

]]> 2. El caudal de extracción de cada nodo de bombeo del modelo a optimizar npa no debe sobrepasar la capacidad máxima proyectada para los pozos asociados con el nodo; aunque al menos debe entregar una demanda local especificada.

Dicha serie de restricciones genera T x npa ecuaciones (por cada límite), es decir, 2T x npa ecuaciones en total:

3. Los niveles piezométricos en cada nodo de bombeo del modelo a optimizar npa para cada tiempo deben ser controlados, esto sugiere que estén por encima de un valor límite prefijado por el modelador. En este caso se generan npa ecuaciones:

Donde:

h_j^k: nivel piezométrico en el nodo de extracción a optimizar j.

h^k_{min j}: nivel crítico en el nodo de extracción a optimizar j. Normalmente es un dato, valor constante en el tiempo.

Las funciones objetivo definidas en las ecuaciones (9), (10) y (11) son lineales y en este sentido es posible emplear un algoritmo de programación lineal para resolver los problemas de optimización que sugieren dichas expresiones.

El modelo de administración MADA se basa entonces en introducir todas las ecuaciones del modelo de simulación AQÜIMPE dentro de un algoritmo de programación lineal, empleando para tal efecto la función Linprog de Matlab.

Es importante plantear que este modelo será efectivo para acuíferos confinados o libres que en su modelación puedan ser tratados como confinados (Gutiérrez, 2008).

El modelo de programación lineal quedaría de la siguiente forma:

• Función objetivo: cada vez que se ejecute el modelo se desea resolver un problema específico que tendrá que ver con algunas de las expresiones (9), (10) u (11). Además, se debe tener en cuenta que las variables de decisión del modelo son los caudales de explotación en los nodos que se quieran optimizar en cada intervalo de tiempo, Q_j^k, y también los niveles piezométricos en todos los nodos en cada intervalo de tiempo.

En esta función objetivo se debe declarar un vector de costo c (Gutiérrez, 2008), el mismo tiene una fila y ((n x T) + (T x npa)) columnas. Todos sus elementos son cero, excepto los que participen en la función objetivo que se esté analizando, que toman un valor de -1.

• Restricciones de igualdad: éstas vienen expresadas en Linprog con la ecuación (7) y tienen una gran importancia, ya que es en este sistema de ecuaciones donde se incluyen todas las ecuaciones del modelo de simulación, es decir, el modelo no encuentra los niveles en los nodos tiempo a tiempo con la ecuación (3), como se realiza con AQÜIMPE en el problema directo, sino que tiempo a tiempo el sistema (3) se va incluyendo dentro de otro sistema mayor, que es el representado en la ecuación (7), quedando una gran matriz A_eq con (n x T) filas y ((n x T) + (T x npa)) columnas, donde están los coeficientes que afectan a todas las variables de decisión. El vector r_eq queda con ((n x T) + (T x npa)) filas y una columna.

Cuando se resuelve el sistema planteado con la programación lineal se obtienen todos los niveles piezométricos en los nodos y en todos los tiempos, de una vez.

]]> • Restricciones de desigualdad: las restricciones de desigualdad son inecuaciones que permitirán introducir en el modelo las importantes restricciones que se abordaron en (12), (13) y (14). La restricción (12) se introduce mediante la matriz A y el vector r dados en la inecuación (6), la cual está concebida por la función Linprog de Matlab como una inecuación del tipo menor o igual que. Cuando las inecuaciones restrictivas de (12) son del tipo mayor o igual que, se multiplican todos los elementos de A y r por -1.

La matriz A queda con (T) filas y ((n x T) + (T x npa)) columnas y los vectores x y r quedan con ((n x T) + (T x npa)) filas y una columna, y (T) filas por una columna, respectivamente.

Para introducir las otras restricciones, dadas en las ecuaciones (13) y (14), se utilizan los vectores de fronteras inferiores y superiores [lb} y {ub}, cada uno de estos vectores tiene dimensiones de ((n x T) + (T x npa)) filas y una columna. El valor de cada elemento de los vectores es el mínimo o máximo, respectivamente, que puede tener su variable de decisión asociada.

Discusión

Las experiencias que se tenían en Cuba en la administración óptima de la explotación del agua subterránea empleaban las funciones tecnológicas algebraicas (las cuales se obtienen mediante un trabajo previo del usuario con el modelo de simulación AQÜIMPE), para posteriormente aplicar el algoritmo de programación lineal. Sin embargo, en MADA esto no es necesario, ya que en el proceso de optimización se tiene al modelo de flujo propiamente dicho, lo cual es una ventaja, ya que mejora todo el trabajo de modelación a realizar por el usuario, así como el tiempo de cómputo. El precio a pagar es la complejidad de la concepción e implementación de este modelo, que se le suele llamar modelo acoplado (Pérez, 1995; Yaragal y Nagaraj, 2007).

Otras de las virtudes del modelo MADA es que por emplear AQÜIMPE, que usa el triángulo cuadrático como elemento finito, tiene una muy buena precisión en los resultados que se obtienen. Por lo general, los modelos basados en el método de los elementos finitos que se encuentran en la literatura no usan el triángulo cuadrático como elemento finito.

Además permite:

• Declarar recargas al acuífero en todos los tiempos y cada elemento.

]]> • Declarar flujos laterales de entrada o salida.

• Declarar nodos con gastos puntuales conocidos (que no serán optimizados).

• Declarar nodos con cargas fijas en todos los tiempos.

• Declarar nodos con cargas variables en todos los tiempos.

• Limitar inferior y superiormente cualquiera de las variables de decisión (cargas y caudales), e incluso a combinaciones de ellas.

• Amplia variedad de salidas gráficas, incluidas animaciones.

Modelo regional del acuífero "Cuenca Sur" de La Habana

Características de la zona de ubicación del acuífero

La zona de estudio se encuentra en el sur de la provincia de La Habana, abarcando un área de 830 km² (ver figura 1). Existe en esta zona una gran carsificación de las rocas del complejo carbonatado de mioceno, una diversidad de las formas del carso, compleja relación con la tectónica, una abundante alimentación al acuífero por filtración y poco espesor de la cobertura de las rocas permeables de la zona de aireación, entre otros factores, que crean condiciones favorables para la formación de considerables recursos de agua subterránea (Llanusa, 1997).

1. Gran heterogeneidad y anisotropía de las rocas con diferentes grados de carsificación.

2. Cuenca abierta con descarga al mar.

3. Compleja composición química del sistema acuífero desarrollado en las calizas y dolomías del mioceno.

4. Extracción intensa e irregular del agua. Se realiza principalmente en el periodo seco.

A lo largo de la costa y con una faja de ancho de aproximadamente uno o tres kilómetros, predomina el suelo cenagoso. El cultivo principal en la región es la caña, se observan también cultivos menores y pastos. Hacia el sur se tiene la vegetación típica costera representada por distintas especies: mangle, patabán y otras.

Orográficamente la región es muy llana, las elevaciones oscilan entre 100 m al norte y 3-5 m al sur. La red fluvial no está bien definida, y es prácticamente imposible delimitar las cuencas hidrológicas. Hacia el sur existe una amplia red de canales, cuya alimentación proviene en lo fundamental de las aguas subterráneas (Llanusa, 1997).

Modelo conceptual propuesto

Toda la modelación matemática del acuífero está basada en un conjunto de abstracciones en síntesis de las características físicas reales del sistema natural que será estudiado con fines de aplicación práctica. Los criterios del modelo conceptual se exponen a continuación:

• El acuífero puede ser simulado con un modelo bidimensional no estacionario en planta, esto quiere decir que sus propiedades hidrogeológicas no varían en dirección vertical.

• El medio filtrante puede simularse como si fuera un medio continuo poroso isotrópico y heterogéneo.

• En el modelo de administración, el acuífero será simulado como confinado.

Discretización de la zona de estudio

En la conformación de los elementos de la malla discretizada se tuvo en cuenta que tanto la geología como las propiedades hidrogeológicas se mantuvieran homogéneas en esos elementos. En la zona paralela a la costa se emplearon elementos que cubrieran el ancho de esta franja cenagosa, no pudiéndose utilizar elementos más pequeños por la carencia de información. La zona de estudio quedó dividida en 69 triángulos, con un total de 158 nodos, como se observa en la figura 1.

Figura 2

Condiciones de frontera

La zona objeto de estudio se encuentra en la vertiente sur de la provincia de La Habana. Por el norte limita con el parteaguas del sur de la Cuenca de la Laguna de Ariguanabo y parte de la divisoria central del país. En esta zona norte, el acuífero se encuentra en contacto con formaciones poco permeables, por lo que fue considerado como un contorno impermeable.

Las fronteras este y oeste del modelo pueden ser asumidas como impermeables debido a que regionalmente existe un flujo del norte hacia el sur en estas zonas que no permite entradas ni salidas de agua subterránea al acuífero por dichas fronteras.

Propiedades hidrogeológicas de la zona de estudio

Se consideran como propiedades hidrogeológicas del acuífero: la conductividad hidráulica, K y la porosidad efectiva o coeficiente de almacenamiento S. Las propiedades hidrogeológicas de los elementos que conforman la zona de estudio son obtenidas de la etapa de calibración dadas en 36 subzonas, columna Material, en el cuadro 1 (Llanusa et al, 1993).

Niveles piezométricos del acuífero: estado inicial

El estado inicial en el acuífero se obtuvo a partir de los valores medidos en los pozos de observación y su cota topográfica. Estos datos se tomaron de una serie de 19 pozos ubicados en la zona acuífera de Güira-Quivicán y pertenecientes a la red informativa de niveles.

Para definir el estado inicial en la modelación se tomaron los datos pertenecientes a la fecha (noviembre del año 1993), publicados en Hernández (2000); estos datos definen una carga media a lo largo de la costa de 0.57 m, producto de la presencia del dique sur.

Alimentación del acuífero debido a la lluvia: recarga

Batista (1995) hace un análisis de la lluvia media hiperanual desde el año 1964 hasta 1993 con la información de 12 pluviómetros. De este análisis determinó la lluvia del 75% de probabilidad (1 113 mm), que prácticamente coincidió con la media obtenida en el año 1975 (1 114 mm). Posteriormente tomó la distribución de lluvia de este año y realizó un balance hídrico diario del suelo de cobertura de la zona de estudio, empleando para esto el modelo matemático EVABAL. Finalmente obtuvo la distribución espacial de la recarga al acuífero debido a la lluvia, así como su cuantificación en volumen. Esta distribución quedó organizada en 12 grupos de triángulos, en los cuales la recarga mensual era diferente. Estos resultados se emplean en el presente trabajo, y se muestran en la figura 3 y en el cuadro 2.

Tratamiento de la explotación al acuífero

Los volúmenes de explotación por nodos que serán empleados en las nuevas simulaciones fueron tomados de Batista (1995) y se muestran en el cuadro 3.

En este cuadro se colocan primeramente las principales extracciones realizadas al acuífero, las cuales corresponden al acueducto y a la agricultura. La suma de estas extracciones asciende a 264.10 hm³/año, lo que representa un 90% de la extracción global, 295.20 hm³/ año. El resto está muy diseminado en otros 64 nodos de la malla; en la última fila del cuadro 3 se colocaron los acumulados de estos valores para este grupo de nodos.

Análisis de la demanda de explotación global solicitada al acuífero

Una forma de examinar si el acuífero puede satisfacer la explotación global solicitada es realizando una simulación con el criterio de que al final del año el volumen de agua en el acuífero sea similar al del inicio de la simulación, concepto equivalente al rendimiento seguro de la escuela norteamericana (Hernández, 1991).

En la figura 4a se muestra un gráfico de volúmenes en el acuífero en cada mes del año simulado con el modelo AQÜIMPE. Es interesante observar que bajo el régimen de explotación y de recarga antes definidas, el acuífero al final del año queda con un volumen ligeramente menor al del inicio de la simulación.

Un análisis más consistente pudiera ser el de repetir este año patrón varios años seguidos, para observar la tendencia de regulación del acuífero. Por ejemplo, se simularían cinco años de baja recarga al acuífero (75% de probabilidad), con una explotación constante (igual al PUA del año 1995). El gráfico de volúmenes de esta variante se muestra en la figura 4b.

En la figura 4b se aprecia una tendencia de estabilidad del acuífero, lo que demuestra que el Plan de Uso del Agua de 295.20 hm³/ año, manejado por Batista (1995), tiene un gran valor práctico, al garantizar el equilibrio hidráulico en el acuífero bajo la lluvia del 75% probabilidad.

En la figura 5 se muestra un mapa de niveles del agua subterránea correspondiente al tiempo cinco de las simulaciones realizadas. En la figura 4a se observa que es el tiempo en el primer año en que el acuífero acumuló menor volumen y por lo tanto se deben esperar en toda la región menores niveles piezométricos. De este mapa es importante resaltar que existe un conjunto de nodos, ubicados hacia la izquierda y en el centro, en donde las cargas hidráulicas están por debajo del nivel cero o nivel medio del mar.

Es por esta razón que en el presente trabajo se aplica un modelo de administración, para proponer una política de explotación que garantice una demanda global aceptable y que conserve los recursos hidráulicos. Además, debe estar distribuida de una forma eficaz en cuanto a garantizar que, en ningún pozo y en ningún tiempo, existan niveles por debajo del nivel medio del mar.

Aplicación del modelo MADA en la optimización de la explotación del acuífero

La primera variante analizada emplea un intervalo mensual, utilizando el mismo estado inicial y la recarga dados anteriormente. Se tendrá en cuenta una extracción global conocida, 295.20 hm³/año.

Serán variables de decisión en el problema los caudales de los 36 nodos que se muestran en el cuadro 3 en cada uno de los tiempos y las cargas piezométricas en todos los nodos y tiempos.

En el cuadro 4 se definen las extracciones mínimas de un grupo de nodos que tienen asociado caudal de extracción para acueductos, así como la máxima capacidad de extracción en ellos.

Por otra parte, el cuadro 5 muestra los valores de extracción mínima y máxima de grupos de nodos vinculados con la explotación con fines agrícolas. Por ejemplo, la suma de los caudales extraídos en los nodos 3, 4 y 5 en cada tiempo está entre 0.01 y 5 hm³/mes.

Los límites de los niveles para todos los nodos en todos los tiempos, exceptuando los que se encuentran en el cuadro 6 coincidentes con pozos de observación, están definidos entre 0 e infinito (∞).

Para cada uno de los intervalos de tiempos (meses) se define un valor de demanda acorde con los datos de extracciones manejados del PUA de 295.20 hm³/año, estos valores organizados por meses, en hm³/mes, son lo siguientes: 28.0, 28.0, 28.0, 34.0, 20.0, 18.0, 18.5, 18.0, 16.5, 17.5, 18.6, 19.0.

Obsérvese que la suma de ellos da un valor de 264.10 hm³/año, el cual es menor que el del PUA y esto tiene sentido, ya que la diferencia es de alrededor del 10% y corresponde con un caudal que no será optimizado, o sea que es conocido e impuesto, última fila del cuadro 3.

La función objetivo utilizada en este estudio es la que responde a la ecuación (11), o sea, implica maximizar la suma de los caudales en los 36 nodos de extracción, en todos los tiempos.

En esta primera optimización, las restricciones de demanda se definieron como menor o igual que. Esto garantiza, por ejemplo, que en el primer tiempo, la suma de todos los caudales en los 36 nodos de optimización sea como máximo igual a 28.0 hm³/mes, y como la función objetivo implica maximizar el caudal, es de esperar que ese sea el resultado. En este caso, es interesante el resultado de la distribución que propone el modelo, que no tiene que coincidir con la forma de explotación real, dada a priori al acuífero, expuesta en el cuadro 3.

Los resultados de los caudales, una vez ejecutado el modelo MADA, se muestran en el cuadro 7 y los volúmenes obtenidos en el acuífero para cada tiempo en la figura 6a.

Es importante recordar que con esta distribución de caudales se cumplen todas las restricciones definidas anteriormente, lográndose que no existan niveles negativos en ninguna zona del acuífero.

Como puede observarse en la figura 6b, los volúmenes acumulados son menores y esto se debe a una fuerte extracción, que ascendió a unos 339.28 hm³/ año, lo que representa un 128.5% del volumen necesario (264.10 hm³/ año) para garantizar el PUA. Esta excesiva explotación lleva al acuífero a terminar en una situación crítica en los niveles al fin del año.

Este resultado obtenido tiene un gran valor práctico, ya que brinda una idea del máximo caudal anual posible a extraer bajo las condiciones de un año seco (75% de probabilidad de lluvia) y garantizando niveles piezométricos por encima del nivel medio del mar.

Finalmente, se realizó una optimización de la explotación, repitiendo el mismo año hidrológico cinco veces, por lo que se tomaron en cuenta en total sesenta intervalos de tiempo. Lo que se busca ahora con este análisis es apreciar el comportamiento del acuífero bajo la acción de la demanda solicitada en un periodo seco prolongado, en este caso de cinco años.

Los resultados obtenidos de la simulación en cuanto al cumplimiento de la demanda exigida cada mes en cada uno de los cinco años se presentan en el cuadro 8. En dicho cuadro se observa que en el primer año se garantiza la demanda planteada, sin embargo en los meses tres y cuatro del segundo año, no se puede entregar la demanda y es lógico que sea precisamente en estos meses que comience la crisis, ya que en los primeros cuatro meses del año no hay recarga al acuífero (ver cuadro 2), sin embargo, se le pide una extracción considerable, sobre todo en el mes cuatro, que asciende a 34.0 hm³ y este efecto puede ser soportado el primer año, pero posteriormente todos los años presentan déficit. Los meses con déficit en cada año pueden observarse en el cuadro 8, en las celdas marcadas.

Con estos resultados se demuestra que el acuífero no puede cumplir con las demandas mensuales planteadas, por lo tanto éstas deberán disminuir, de forma tal que se pueda encontrar una distribución máxima de las demandas que garantice la no ocurrencia de años con déficit. En este caso, la distribución encontrada mediante un ensayo de prueba y error es la siguiente 28.0, 23.0, 20.0, 20.0, 19.0, 18.0, 18.5, 18.0, 16.5, 17.5, 18.6, 19.0.

La suma de estos valores es de 236.10 hm³/ año. La diferencia entre el PUA y este valor es de 28.0 hm³/año. Pudiera pensarse en variar la distribución en función de los objetivos que se quieran garantizar, sin embargo, con esta distribución de demandas por meses que se obtuvo, se garantiza la entrega para las condiciones exigidas.

Análisis de dimensiones y tiempos de ejecución con el modelo MADA

Un aspecto interesante en la última ejecución del modelo MADA está referido a la dimensión de la matriz A_eq , que contiene al modelo AQÜIMPE en el algoritmo de programación lineal. Para este ejemplo, dicha matriz necesitó de 9 480 filas y 11 640 columnas para un total de 110 347 200 elementos. El programa está preparado a tal efecto para trabajar con matrices de tipo "sparse", con las cuales se reduce el espacio de almacenamiento de los datos, ya que se almacenan sólo los valores diferentes de cero de las matrices.

A pesar de las dimensiones del problema solucionado, el método solamente empleó unos dos minutos, con un total de 49 iteraciones en un ordenador personal con microprocesador Intel Core 2 Duo a 2.40 GHZ y 2 GB de memoria RAM, con lo cual se puede concluir que cuando el problema a resolver tiene grandes dimensiones —si bien el método de programación lineal necesita almacenar mucha información— el resultado se obtiene en un tiempo muy corto y con la garantía de que si el problema está bien planteado la solución será la correcta.

Conclusiones

1. Se realiza la concepción teórica general de un modelo de optimización para la administración de la explotación de acuíferos hidráulicamente confinados, llamado MADA. La originalidad del modelo propuesto se basa en introducir las ecuaciones, que emplea el modelo de simulación de acuíferos AQÜIMPE, como restricciones de igualdad, al emplear la programación lineal como técnica de optimización. Además, se tienen en cuenta un conjunto de restricciones a las variables de decisión que enriquecen el valor práctico de este modelo.

2. Queda demostrada la factibilidad del modelo propuesto MADA, al llevarse a cabo una evaluación crítica de la explotación en el modelo regional de la Cuenca Sur de La Habana. En este caso de estudio se demostró que los resultados que brinda el modelo MADA para la distribución de la explotación que garantiza el Plan de Uso del Agua, evalúa las posibilidades de intrusión salina en el acuífero, o sea, se obtiene una distribución de caudales que respeta los valores piezométricos regionales críticos en el acuífero.

3. Se comprueba que el nivel de extracciones correspondiente al Plan de Uso del Agua puede verse afectado por el efecto de varios años secos continuados. En este sentido se realiza una optimización de las extracciones para cinco años secos consecutivos. Los resultados muestran años con meses en los que no se cumple la demanda (meses con déficit). Además se muestra la flexibilidad del modelo para encontrar soluciones a estos problemas, dependiendo de los objetivos que persiga el modelador.

Agradecimientos

Se realiza un especial reconocimiento al doctor José Bienvenido Martínez Rodríguez, creador de la tecnología AQÜIMPE, por sus certeros consejos para acometer la tarea de concebir este modelo. Así como a los especialistas, los doctores Armando Hernández Valdés y David Marón Domínguez del Centro de Investigaciones Hidráulicas (CIH) y del Departamento de Matemáticas de la Facultad de Ingeniería Civil de la CUJAE, respectivamente. También se agradece a la doctora Haydée Llanusa Ruiz por brindar gran parte de la información referente a la Cuenca Sur de La Habana.

BATISTA, C.J. Nuevas experiencia en la modelación de GÜIRA-QUIVICÁN. Tesis de maestría. La Habana: Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 1995.         [ Links ]

CABRERA, E. Simulación de acuíferos con el empleo de herramientas de los sistemas de información geográfica. Tesis de maestría. La Habana: Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 2007.         [ Links ]

DILLA, F. Modelación matemática del flujo subterráneo y las técnicas de optimización aplicadas en la investigación y explotación de acuíferos a escala regional. Capítulo V. Tesis doctoral. La Habana: Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 1993.         [ Links ]

GALLOWAY, D.L., ALLEY, W.M., BARLOW, P.M., REILLY, T.E., and TUCCI, P. Envolving issues and Practices in Managing Ground-Water Resources. Case Studies on the Role of Science. U.S. Department of Interior. U. S. Geological Survey, Circular 1247, 2003, pp. 19-24.         [ Links ]

GARCÍA, J., RODRÍGUEZ, I.J. y VIDAL, J. Aprenda Matlab 7.0 como si estuviera en primero. Madrid: Escuela Técnica Superior de Ingenieros Industriales, Universidad Politécnica de Madrid, diciembre de 2005.         [ Links ]

GÓMEZ, M. Plataforma modular integrada para la gestión del agua subterránea como recurso. Tesis doctoral. La Habana: Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 2009.         [ Links ]

GUTIÉRREZ, V.A. Propuesta de un modelo matemático para la administración óptima de acuíferos. Tesis de grado. La Habana: Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 2008.         [ Links ]

HERNÁNDEZ, A.O. La explotación de acuíferos a escala regional y la modelación matemática como su base tecnológica. Tesis en opción del grado científico de Doctor en Ciencias Técnicas. La Habana: Instituto Superior Politécnico "José Antonio Echeverría", 1991.         [ Links ]

HERNÁNDEZ, C.L. Modelación matemática del acuífero costero cuenca sur de La Habana con información actualizada. Trabajo de diploma. La Habana: Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 2000.         [ Links ]

HSIAO, C.T. and CHANG, L.C. Dynamic Optimal Groundwater Management with Inclusion of Fixed Costs. Journal of Water Resources Planning and Management. January/February, 2002. pp. 57-65.         [ Links ]

LLANUSA, H., MARTÍNEZ, B.J., BATISTA, C.J. y HERNÁNDEZ, A. Tecnología de simulación de los acuíferos costeros Artemisa-Quivican. Informe Técnico. La Habana: Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 1993.         [ Links ]

LLANUSA, H. Modelación matemática de la intrusión salina. Tesis presentada en opción al grado de Doctor en Ciencias Técnicas. La Habana: Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 1997.         [ Links ]

LLANUSA, H., HERNÁNDEZ, A., MARTÍNEZ, B.J. y DILLA, F. Modelo del acuífero costero sur de La Habana: Evaluación del dique. Primer taller de expertos en vulnerabilidad de acuíferos costeros en el caribe insular. La Habana, 2004.         [ Links ]

MARRERO, N. Técnicas de Optimización Aplicadas a la Ingeniería Hidráulica. La Habana: Publicación del Ministerio de Educación Superior (MES), 1985.         [ Links ]

MARTÍNEZ, B.J. Simulación matemática de cuencas subterráneas: flujo impermanente bidimensional. Monografía. La Habana: Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 1989.         [ Links ]

MARTÍNEZ, B.J., LLANUSA, H., HERNÁNDEZ, A., DILLA, F. Manual del usuario del sistema AQÜIMPE. La Habana: Centro de Investigaciones Hidráulicas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 2000.         [ Links ]

MARTÍNEZ, J.B. Comunicación personal. Centro de Investigaciones Hidráulicas CIH, Instituto Superior Politécnico "José Antonio Echeverría", junio de 2006.

MESA, R.H. Solución del problema inverso en modelos de flujo del agua subterránea mediante un algoritmo de convergencia global. Tesis doctoral. La Habana: Departamento de Matemáticas, Facultad de Ingeniería Civil, Instituto Superior Politécnico "José Antonio Echeverría", 2004.         [ Links ]

PÉREZ, M.F. La hidráulica de captaciones y la modelación matemática regional de acuíferos. Tesis Doctoral en Ciencias Técnicas. La Habana: Instituto Superior Politécnico "José Antonio Echeverría", 1995.         [ Links ]

YARAGAL, S.C. and NAGARAJ, M.K. Combined optimization simulation model for groundwater management policy options [en línea], 2007. Disponible en World Wide Web: http://www.ieindia.org/pdf/88/88CV202.pdf.         [ Links ]