nueva página del texto (beta)
Español (pdf)
Artículo en XML
Referencias del artículo
Traducción automática
Enviar artículo por email
Citado por SciELO 
Similares en
SciELO 
Permalink
Introducción
El escurrimiento proveniente de corrientes superficiales y manantiales ―en ausencia de precipitaciones― es alimentado por las reservas subterráneas de los acuíferos, las cuales se han acumulado en dichas estructuras geológicas en el transcurso de una recarga determinada (Pérez-Franco, 1995).
Este proceso puede ser simulado a partir de modelos matemáticos (agregados o distribuidos), que establecen relaciones entre la precipitación y el escurrimiento. Este tipo de modelos se ha venido utilizando para evaluar los recursos hídricos renovables de un sistema subterráneo, dependiendo de la complejidad y finalidad de los modelos de gestión adoptados. Uno de los elementos que soporta la elección de los modelos de referencia es la disponibilidad de datos. No siempre se cuenta con datos piezométricos y de explotación que permitan utilizar modelos distribuidos, por lo que en varias ocasiones, a pesar de la variabilidad hidráulica del acuífero en estudio, debe recurrirse razonablemente a modelos agregados, los cuales permiten, utilizando pocos parámetros, globalizar el comportamiento del acuífero a partir de datos en ríos o manantiales relacionados con las recargas y/o extracciones con cierta precisión.
De acuerdo con Arenillas et al. (s.f.), estos modelos agregados de precipitación-escurrimiento (MPE) forman parte de un grupo específico de herramientas computacionales para la simulación hidrológica, cuyo principal objetivo es la generación (o reconstrucción) de series temporales de escurrimiento, partiendo de una serie de datos de precipitaciones y extracciones por bombeo, así como temperatura existente en la zona de estudio.
En dichos modelos intervienen parámetros físicos que son característicos del sistema que se pretende simular y, en sentido general, son desconocidos, resultando de suma importancia abordar la formulación de un problema inverso que conduzca a una estimación óptima de los mismos, pues intervienen en la evaluación de los recursos hídricos del sistema en cuestión.
La solución de un problema inverso puede ser hecha a través de algoritmos de optimización, con el propósito de minimizar (o maximizar) cierta función objetivo. Básicamente existen dos clases de métodos de optimización: determinísticos y estocásticos. Los primeros dependen del gradiente de la función objetivo, haciendo uso de las particularidades del funcional; en cambio, los métodos estocásticos trabajan con la selección de configuraciones aleatorias; por ejemplo, técnicas computacionales inspiradas en la naturaleza, con el objetivo de barrer (estocásticamente) y el espacio de búsqueda completo, al procurar el óptimo global (Silva-Neto & Becceneri, 2009).
En el presente artículo se abordarán diferentes modelos agregados para el estudio de descargas de acuíferos. Para establecer una comparación entre ellos, el autor inicialmente unificará el enfoque teórico que permitirá obtener las respectivas expresiones que estimarán los volúmenes de descarga en acuíferos kársticos. En la sección 2, al inicio de este estudio, será descrito el modelo conceptual MEDA, propuesto por Iglesias (1984), el cual aplica la fórmula propuesta por Maillet (citado por Pérez Franco, 1995) para el estudio del agotamiento del escurrimiento proveniente de un acuífero. Otras formulaciones, como las desarrolladas por Tisson et al., así como Forkasiewicz y Paloc (citados por Pérez-Franco, 1995), también serán objeto de análisis, con el propósito de formular modelos que permitan estimar las descargas en cuestión. En la sección 3, el autor presentará y desarrollará un modelo no lineal. La sección 4 describe el problema inverso donde es formulada la función objetivo y los algoritmos de búsqueda global utilizados, mientras que en la sección 5 se describe el caso de estudio. Posteriormente se muestran los resultados obtenidos y su discusión. Por último se tienen las conclusiones derivadas de los objetivos perseguidos en esta contribución.
Modelos agregados para cuantificar caudales de descarga de acuíferos
El caudal de determinada fuente, o la descarga de un acuífero en sentido general, no se correlaciona con la precipitación con carácter mensual. Es decir, las descargas dependen del volumen almacenado en el acuífero y éste guarda correlación con las precipitaciones precedentes en varios meses e incluso en años anteriores. En cambio, si se consideran la precipitación de cierto periodo y el volumen de descarga acontecido a lo largo del tiempo como consecuencia de esa precipitación es posible establecer una clara interrelación. El método para determinar el volumen de descarga consecuente a una precipitación dada ha sido abordado en detalle por Castillo (1981), citado por Arenillas et al. (s.f.). Con base en series de precipitación (P) y volúmenes observados (V), se han obtenido ecuaciones de ajuste de tipo lineal y potencial, con coeficientes de correlación (r) del orden de 0.9 en la forma:
siendo m y n = parámetros de ajuste. En función del tipo de ecuación adoptada, los parámetros m y n presentan dimensiones. Esta observación es pertinente, pues en Arenillas et al. (s.f.), ambos parámetros son literalmente presentados como adimensionales, lo cual constituye, sin duda, un error en dicha publicación. En efecto, para el enfoque lineal (m) = L 2 y (n) = L 3, mientras que para el potencial (m) = L 3 - n y (n) = adimensional, al tratarse del exponente.
Formulación de Iglesias
Tal y como fue enunciado con anterioridad, Iglesias (1984) introduce el modelo propuesto por Maillet (citado por Pérez Franco, 1995) para la curva de agotamiento de la descarga de un acuífero, el cual se expresa como:
siendo Q t = descarga del acuífero en el instante de tiempo t; Q o = descarga al inicio del periodo analizado; 𝛼 = coeficiente de agotamiento, característico del acuífero. Acorde con Pérez-Franco (1995), este parámetro 𝛼 = 1/𝜏, siendo 𝜏 el tiempo necesario para que Q t = Q o /e ≅ 0.3678 Q o .
Tal y como se muestra en la Figura 1, en
el i-ésimo mes (mes i), el caudal de
descarga se denota por Q
i
. En caso de no existir recarga como consecuencia de precipitaciones,
entonces el caudal del mes i + 1 será denotado como
La precipitación de referencia P i , acorde con Iglesias (1984), habrá de incrementar el volumen almacenado en el acuífero en una cantidad equivalente al área comprendida entre las curvas de agotamiento representadas en la Figura 1 (Arenillas et al., s.f.).
En efecto, aplicando la ecuación (2), el área bajo las curvas de agotamiento identificada como el volumen de descarga V puede ser expresada por la diferencia de las siguientes integrales:
siendo Δt i = intervalo entre el mes i + 1 y mes i. Como se pretende determinar Q i+1 , cualesquiera de las expresiones (1) puede ser utilizada. Arenillas et al. (s.f.) recomiendan el empleo de la ecuación (1b) correspondiente a la estructura potencial de forma generalizada, ya que la ecuación (1a) tiene un rango de validez restringido, y puede llegar a dar resultados anómalos e incluso negativos. Para evitar dicho inconveniente y dar continuidad al desarrollo teórico aquí presentado, el autor será consecuente con dicha recomendación, sin ser excluyente con el posible empleo de la ecuación (1a).
Sustituyendo la ecuación (1b) en (3) y despejando Q i+1 es posible obtener la siguiente expresión:
La ecuación (4) es la base conceptual del Modelo para el Estudio de Descargas de Acuíferos (MEDA) del Instituto Tecnológico GeoMinero de España, y con el mismo es posible realizar las estimaciones pertinentes de los caudales de descarga de cada mes, sobre la base del conocimiento de precipitación, y caudales del mes precedente y los parámetros m, n y 𝛼. Posteriores análisis del modelo han aportado mejoras notables en su formulación conceptual. Entre éstas se puede mencionar la introducción del concepto de precipitación útil (P u ), definida según Arenillas et al. (s.f.), como la diferencia entre la precipitación (P) y la evapotranspiración real (E tr ), la cual forma parte de la recarga en el acuífero y está estrechamente relacionada con la temperatura sobre la superficie del suelo (T) y situación geográfica, entre otros parámetros del clima. Arenillas et al. (s.f.), al analizar valores de evapotranspiración potencial y real en España, presentan la siguiente ecuación de cálculo de la precipitación útil para el i-ésimo mes:
donde el parámetro b = exponente, que puede variar entre 1.3 y 1.6. Expresiones similares a la anterior pueden ser obtenidas para condiciones locales, características de regiones y climas concretos. No obstante, de acuerdo con Arenillas et al. (s.f.), la propuesta abordada no resta generalidad al modelo MEDA y, por lo tanto, la ecuación (4) podrá ser rescrita como:
De igual forma, es posible incluir el volumen de los bombeos extraídos mensualmente en el acuífero (W). Esta inclusión, tal y como reportan Arenillas et al. (s.f.), presenta la limitación, en ocasiones insalvable, de que éstos no produzcan afectaciones de carácter dinámico al punto de descarga. Por tanto, de ser considerado el volumen bombeado en el i-ésimo mes (W i ), la ecuación (6) podrá ser replanteada como:
La ecuación (7), para las consideraciones establecidas hasta su conceptualización, representa la expresión más general y utilizable del modelo MEDA (de tipo MPE), la cual, en ausencia de datos de temperatura y/o bombeos mensuales, presenta como casos particulares las ecuaciones (4) y (6), respectivamente.
Formulación de Tisson, Werner y Dunsquit
Pérez-Franco (1995) ha desarrollado la siguiente expresión para el agotamiento de la descarga de un acuífero:
De forma análoga a la propuesta de Maillet, empleada por Iglesias (1984) en el desarrollo de MEDA, el parámetro α se identifica como el coeficiente de agotamiento del acuífero. Como ya es conocido 𝛼 = 1/𝜏; no obstante esta vez 𝜏 es el tiempo necesario para que Q t = Q o /4. Este razonamiento sugiere que para el mismo valor de 𝛼, las expresiones (2) y (8) arriban a caudales de descarga diferentes, siendo mayores aquellos calculados a través de la ecuación (2). Para obtener el volumen de descarga V, con el empleo de la ecuación (8) se procede de forma análoga al cálculo de las siguientes integrales:
Tomando en cuenta las mejoras introducidas por Iglesias (1984) (citado por Arenillas et al., s.f.) en la cuantificación de V, es posible obtener un modelo para la estimación de las descargas de un acuífero (ver ecuación (10)), donde es de destacar la alteración que sufre el segundo término de la ecuación (7), que caracteriza el agotamiento del acuífero a partir del i-ésimo mes, como consecuencia de introducir la formulación (8):
Formulación de Forkasiewicz y Paloc
Esta fórmula fue propuesta sobre la base de observaciones en una fuente kárstica y presenta la siguiente estructura (Pérez-Franco, 1995):
donde β = coeficiente de decrecimiento del caudal de descarga. Las unidades
de este parámetro son
La ecuación anterior no sólo sugiere que β depende de la capacidad de agotamiento del acuífero, expresada a través de 𝛼, sino también del caudal de descarga Q o y además del tiempo. Puede obtenerse un resultado similar al emplear las ecuaciones (11) y (8). Sin embargo, su desarrollo no reportaría ningún elemento cualitativo adicional al obtenido por la vía aquí expuesta. Precisamente, en esta multidependencia radica la complejidad de la formulación analizada. Aunque en Pérez-Franco (1995) no se presentan elementos adicionales que permitan esclarecer aún más la naturaleza de dicho parámetro, es de suponer que sería útil para reconstruir series históricas de descargas de acuíferos kársticos en determinado periodo, y por tal motivo ha sido introducida en el presente artículo. El volumen de descarga V en este caso podrá expresarse como:
El modelo que resulta para la estimación de descargas de un acuífero a partir de la ecuación (11) presenta una estructura diferente a las presentadas en las ecuaciones (7) y (10). Este modelo (ver ecuación (14)) aporta además la presencia de un nuevo parámetro (β), que caracteriza el decrecimiento de la descarga de una fuente existente en un acuífero, que como ya se ha explicado, es de naturaleza compleja y su determinación debe ser rigurosa para no incurrir en errores no deseados que afecten la calidad de las estimaciones efectuadas en la descarga.
Propuesta de modelo no lineal para determinar la descarga para un acuífero
Tal y como se ha podido apreciar en las secciones anteriores, existen varios modelos para describir la dinámica del agotamiento de la descarga de un acuífero. Uno de los más usados es el propuesto por Maillet (citado por Pérez-Franco, 1995). Desde el punto de vista de su estructura (ver ecuación (2)), dicho modelo es solución de la siguiente ecuación diferencial ordinaria:
sujeta a la siguiente condición inicial Q(0) = Q o . Independiente de los niveles de precisión que puede brindar este tipo de formulación, una pequeña modificación en la ecuación (15) puede ser introducida, considerando que la variación de la descarga en el acuífero es una función no lineal. Ésta puede ser escrita como:
siendo 𝜂 = exponente y 𝜅 = parámetro que caracteriza el agotamiento del acuífero, y cuyas unidades son L 3(1-𝜂) T 𝜂-2 . La integración de la ecuación (16) conduce a la siguiente solución analítica:
Como ya es usual, el volumen de descarga V podrá ser determinado según:
Para la actual formulación, resulta de interés para el autor expresar las integrales de la ecuación (18) empleando procedimientos apropiados para el cálculo de límites, dado que se está en presencia de integrales impropias, por lo tanto:
En la ecuación (19), para 0 < η < 1, los límites del primer y tercer términos convergen a cero y el volumen de descarga resulta de la diferencia entre el segundo y cuarto término en dicha ecuación. Esta particularidad es aprovechada para establecer la formulación del modelo propuesto en esta sección, de forma análoga, a como se procedió para obtener las ecuaciones (7), (10) y (14). En efecto, el modelo no lineal para establecer las descargas de acuíferos, válido para 0 < η < 1, se expresa por la ecuación (20) siguiente:
Formulación del problema inverso
Acorde con Silva-Neto y Moura-Neto (2005), así como con Silva-Neto y Becceneri (2009), los problemas inversos pueden ser clasificados tomando en cuenta la dimensión del espacio donde el modelo es definido, así como la dimensión del espacio donde los objetos que son estimados pertenecen. Esta clasificación se resume en la Tabla 1 y es útil para brindar una idea a priori de la dificultad del problema en consideración.
Tabla 1 Clasificación de problemas inversos.
| Modelo matemático | Estimación | |
|---|---|---|
| Dimensión finita (constantes) | Dimensión infinita (función) | |
| Dimensión finita * | Tipo I | - |
| Dimensión infinita ** | Tipo II | Tipo III |
* por ejemplo: sistema de ecuaciones algebraicas; ** por ejemplo: sistema de ecuaciones diferenciales.
En los modelos dados por las ecuaciones (7), (10), (14) y (20) es posible estimar de modo simultáneo los parámetros que son desconocidos. Como la cantidad de mediciones debe ser superior (o al menos igual) a la cantidad de incógnitas a ser estimadas (parámetros), su identificación se basa en la resolución del problema inverso, tratado como la minimización de cierta función objetivo en un espacio de búsqueda de dimensión finita.
Función objetivo
La función objetivo usada en este estudio es multicriterio, y ha sido adaptada de Barco, Wong y Stenstrom (2008):
sujeta
Tabla 2 Parámetros de los modelos.
| Modelo | Vector de parámetros |
|---|---|
| Ecuación (7) | θ = (m, n, b, 𝛼)T |
| Ecuación (10) | |
| Ecuación (14) | θ = (m, n, b, β)T |
| Ecuación (20) | θ = (m, n, b, 𝜅, 𝜂) T |
En esta investigación se empleó el algoritmo de búsqueda global Cuckoo (Cuckoo search, CS), desarrollado por Yang y Deb (Yang & Deb, 2009; Yang & Deb, 2010). Tal y como reportan Civicioglu y Besdok (2011), la introducción de este algoritmo puede proveer resultados más robustos que algoritmos ampliamente utilizados, como el algoritmo de optimización por enjambre de partículas (Particle Swarm Optimization, PSO) desarrollado por Clerc y Kennedy (2002), así como el basado en una colonia artificial de abejas (Artificial Bee Colony Algorithm, ABC) desarrollado por Karaboga y Basturk (2007a, 2007b). Además, se ha optado por emplear un enfoque híbrido, acoplando CS con algún método determinístico basado en el gradiente. En este caso, será introducida la herramienta FMINCON del asistente matemático MatLab, la cual consta de varios algoritmos y permite, como bien es conocido por los utilizadores de este asistente, encontrar el mínimo de una función no lineal, multidependiente sujeta a restricciones. En esta contribución será empleado el algoritmo punto interior (interior-point), el cual le incorpora rapidez y precisión al problema inverso abordado, al tratarse de minimizar la función objetivo (ecuación (21)), sujeta a restricciones de los parámetros.
El algoritmo CS
El algoritmo CS fue desarrollado por Yang y Deb (Yang & Deb, 2009; Yang & Deb, 2010) basado en la aplicación metaheurística proveniente de la naturaleza e inspirado en el comportamiento de pájaros fascinantes, no sólo por los sonidos que pueden hacer sino por su agresiva estrategia de reproducción. Acorde con dichos investigadores, algunas especies de aves, como Ani y Guira cuckoos, colocan sus huevos en nidos comunes, donde ellos pueden llegar a remover los huevos de otras especies de pájaros para incrementar la probabilidad de producir sus propios huevos.
Un número considerable de especies de cuckoos (o cucos) presentan un comportamiento parásito al momento de anidar, colocando sus huevos en el nido de otras especies para que éstas se encarguen de cuidar y alimentar a sus polluelos. Pueden llegar a establecerse conflictos directos con los cucos intrusos y si una especie descubre que los huevos en su nido no son propios, simplemente aparta estos huevos o abandona ese nido, y construye otro en un lugar diferente. Acorde con Yang y Deb (2009), algunas especies de cucos, como el New World brood-parasitic Tapera han evolucionado de tal manera que las hembras han podido asemejar en color y patrones de los huevos de un número selecto de especies de pájaros. Esto reduce la probabilidad de que sus propios huevos sean abandonados e incrementa su descendencia. El algoritmo CS emplea las siguientes representaciones:
Cada huevo en un nido representa una solución y el huevo de cuco representa una solución nueva. El propósito es usar las nuevas y mejores soluciones (cucos) para reemplazar aquellas que no son tan buenas en los nidos.
En la forma más simple, cada nido tiene un huevo; sin embargo, el algoritmo puede ser extendido a casos más complicados, en los que cada nido tiene múltiples huevos, representando un conjunto de soluciones.
El algoritmo CS está basado en tres reglas idealizadas:
Cada cuco coloca un huevo a la vez y lo hace en un nido escogido de forma aleatoria.
Los mejores nidos, con alta calidad de huevos, serán transferidos a la próxima generación.
El número de enjambres de nidos es predeterminado (N d ) y el huevo colocado por el cuco es descubierto por las especies de pájaros con una probabilidad p a = ∈ (0, 1).
Además, Yang y Deb (2009) descubrieron que el estilo de búsqueda del camino aleatorio es mejor ejecutado por los vuelos Lévy. El pseudocódigo puede ser resumido como se ilustra en la Figura 2.
La importante ventaja de este algoritmo es su simplicidad, así como también su facilidad para ser implementado computacionalmente. Es válido resaltar que el algoritmo propuesto por Yang y Deb (2009) ha sido modificado por diferentes investigadores en este campo. El autor del presente artículo ha utilizado el código de CS programado en MatLab por Yang, en la Universidad de Cambridge, entre noviembre de 2008 y junio de 2009, con base en la propuesta original del algoritmo. Como ya se mencionó, se ha acoplado a la función FMINCON de MatLab, de manera que, en adelante, para hacer referencia a la estrategia adoptada esta combinación será denominada CS-FMIN. En la Figura 3 se muestra el diagrama de flujo para la estimación de parámetros en modelos de descarga de acuíferos. Para la aplicación de CS-FMIN, inicialmente será utilizada una versión débil de CS, con el objetivo de generar una aproximación inicial para el algoritmo IP. Tal y como muestra la Figura 4, en primer lugar, se establece un procedimiento iterativo para estimar los parámetros del modelo de descarga analizado, pero con el objetivo de mejorar los resultados es que se introduce la estimación definitiva a través de FMINCON.
Figura 3 Diagrama de flujo para la estimación de los parámetros en modelos de descargas de acuíferos.
Figura 4 Ubicación geográfica y cuenca superficial asociada al área del acuífero Ciego-Morón. Adaptado de Vidal-Olivera y González-Abreu (2013).
Materiales y métodos
Estudio de caso
A continuación se expone el análisis de las descargas del acuífero Ciego-Morón, tomando como referencia principal el estudio de Martínez-Rodríguez, Hernández-Valdés, Llanusa-Ruiz y Dilla-Salvador (1986), durante el periodo comprendido entre noviembre de 1983 a diciembre de 1985 (un total de 24 meses). Dichos autores dirigieron su investigación hacia el conocimiento del comportamiento hidráulico del mencionado acuífero, el cual juega un papel fundamental en el desarrollo socioeconómico de la provincia Ciego de Ávila, Cuba. El acuífero Ciego-Morón se establece en las calizas y dolomías del Mioceno, con espesores variables, que pueden alcanzar 300 m, con alta presencia de karst, abarcando más de 80% del territorio analizado (Martínez-Rodríguez et al., 1986).
Entre los paquetes de programas para realizar el estudio del comportamiento de este acuífero se reporta el modelo MEDA (Iglesias, 1984; Arenillas et al., s.f.), con el objetivo de completar la serie histórica de descarga del acuífero en el periodo de calibración seleccionado (de noviembre de 1983 a diciembre de 1985). Al no darse detalles del proceso de calibración por parte de Martínez-Rodríguez et al. (1986) es de suponer que ésta fue realizada de forma manual (mediante ensayo y error), pues no existe evidencia que indique si la versión utilizada por estos autores es aquella que hace uso del algoritmo de Levenberg-Marquardt para la calibración automática de sus parámetros, implementada por Erloza (1985), citado por Arenillas et al. (s.f.).
Definición del área a modelar y datos para la estimación de la descarga
En el área en estudio (ver Figura 4), la descarga del acuífero Ciego-Morón se produce hacia el norte. Los caudales medidos se reportan en la estación de aforo de la cuenca superficial La Yana (área de 1 503 km2) en el periodo comprendido en la sección de referencia (Figura 5b).
Figura 5 a) Registro de precipitaciones pluviómetro 124. Precipitaciones entre diciembre de 1983 y diciembre de 1985, intervalo estudiado por Martínez Rodríguez et al. (1986); b) descargas en la estación de aforo de La Yana, de diciembre de 1983 a diciembre de 1985.
Se seleccionó el periodo de precipitaciones entre enero de 1929 a diciembre de 2010 (Figura 5a), a partir de los datos de un pluviómetro ubicado en la zona central del acuífero e identificado como P-124, coincidiendo con el utilizado por Martínez-Rodríguez et al. (1986), pero para una serie histórica entre enero de 1967 y diciembre de 1986. De acuerdo con dicho autor, su selección se justificó al resultar imposible hallar la isoyeta media mensual, la cual hubiera sido de mucha utilidad para establecer la correlación entre la precipitación mensual y los caudales de descarga observados en la estación de aforo La Yana.
Es válido destacar que en dicho estudio se introdujo el volumen de descarga de acuerdo con una relación del tipo ecuación (1a); además, no se reportaron datos de bombeo y la precipitación útil se toma como el acumulado mensual, es decir, sin considerar la evapotranspiración real, estableciendo de esta manera los datos de entrada en los modelos abordados en esta investigación.
Experimento numérico
Un experimento numérico fue conducido con cada uno de los modelos en estudio, donde la dimensión del problema quedó previamente establecida a partir del número de parámetros, al ser estimados en cada uno de ellos y sus posibles variantes. En el algoritmo CS-FMIN, la tolerancia se estableció por etapas. Para determinarla se condujeron varios experimentos previos, tomando en cuenta que se está empleando una versión débil de CS para generar la aproximación inicial a utilizar en FMIN. La tolerancia en CS está definida como el mínimo valor que puede obtener la función objetivo en estudio, dada en este caso por la ecuación (21), tal y como se muestra en la Tabla 3.
Tabla 3 Experimento numérico.
| Modelo | Tolerancia CS | Nd | pa |
|---|---|---|---|
| Ecuación (7) | 3.2* | 5 | 0.075 |
| Ecuación (10) | |||
| Ecuación (14) |
16* 5.5** |
5 25 |
0.075 0.25 |
| Ecuación (20) |
3.2* 5.5** |
5 25 |
0.075 0.25 |
(*) Considerando parámetros constantes.
(**) Considerando variabilidad temporal del parámetro de agotamiento (p. ej., ecuación (14): β y ecuación (20): 𝜅.
En las ejecuciones preliminares de CS-FMIN se detectó que cuando se emplea el modelo basado en la curva de agotamiento propuesta por Forkasiewicz y Paloc, según la ecuación (14) y considerando el parámetro β como constante, la tolerancia de CS debe ser elevada para arribar a una aproximación inicial que permita la ejecución de FMIN. Esto no sucede en el caso del modelo de descarga no lineal propuesto por el autor, según la ecuación (20), donde la tolerancia se mantuvo en los mismos valores de los modelos dados por las ecuaciones (7) y (10), lo cual es favorable para esta formulación.
Cuando los parámetros β y 𝜅 son considerados variables, tomando en cuenta que la dimensión del problema de estudio crece al incrementarse el número de parámetros a identificar, el autor prefirió ser conservador y admitir una tolerancia ligeramente superior (Tol = 5.5), para agilizar el proceso de búsqueda de la aproximación inicial. Por el contrario, la tolerancia en la verificación de las restricciones en FMIN se mantuvo constante e igual a 1e-06 para todos los casos analizados.
De igual forma, el número de nidos disponibles N d y la probabilidad p a de que el huevo depositado por el cuco fuera descubierto por otras especies de aves se seleccionó a partir de las ejecuciones preliminares a las que se hacía referencia con anterioridad. En efecto, pudo detectarse que resulta ventajoso que N d sea del orden del número de parámetros a ser estimados; por esta razón se escogieron los valores de 5 y 25. Cuando el número de parámetros se elevó de manera considerable, como consecuencia de asignar la variabilidad temporal a β y 𝜅, también fue necesario incrementar la probabilidad p a , resultando los valores de 0.075 y 0.25, respectivamente (ver Tabla 3).
Se estableció un total de 20 réplicas para cada modelo en estudio, resultando un total de 120 ejecuciones de CS-FMIN, las cuales fueron llevadas a cabo en una computadora con 12 Gb de memoria RAM y procesador Intel (R) Core (TM) i7-4500U CPU@ 2.40 GHz. En la Tabla 4 se muestra el dominio que representa el espacio de búsqueda para cada uno de los parámetros a identificar en cada uno de los modelos, los cuales también resultaron de las pruebas preliminares ya mencionadas.
Tabla 4 Experimento numérico para los modelos en estudio.
| Modelo | Dimensión | Dominio |
|---|---|---|
| Ecuación (7) | 3 |
|
| Ecuación (10) | ||
| Ecuación (14) | 3 |
|
| 26 |
|
|
| Ecuación (20) | 4 |
|
| 27 |
|
(*); (**): ÍdemTabla 3, en este particular 24 representa el número de meses.
Resultados y discusión
En la Tabla 5 se muestran los parámetros que, en las réplicas conducidas por el autor, presentaron el mejor desempeño en la evaluación de la función objetivo dada por la ecuación (21). En la Figura 6 se expone el comportamiento de los parámetros β y 𝜅, donde puede apreciarse que para este último (𝜅) existe una alta correlación entre su valor medio mensual y la distribución que minimiza la función objetivo en su mejor ejecución. En la Figura 7 se muestra la comparación de caudales de descarga simulados y observados para los mejores desempeños de cada modelo. En la Tabla 6 se presenta el procesamiento estadístico (media μθ, desviación estándar σθ y coeficiente de variación Cvθ) de las 20 réplicas ejecutadas, donde se detecta que el parámetro m exhibe las mayores variaciones con relación a su valor medio.
Tabla 5 Parámetros con mejor desempeño en la evaluación de la F.O.
| Modelo | m | n | 𝛼 (e-07) | β (e-06) | 𝜅 (e-07) | 𝜂 | F.O (ecuación 21) |
|---|---|---|---|---|---|---|---|
| Iglesias (1984) | 68.908 | 2.137 | 2.86 | - | - | - | 3.051 |
| Tisson et al. (Pérez-Franco, 1995) | 63.341 | 2.247 | 1.74 | - | - | - | 2.871 |
| Forkasiewicz y Paloc * (Pérez-Franco, 1995) | 79.919 | 1.787 | - | 5.671 | - | - | 15.962 |
| Forkasiewicz y Paloc ** (Pérez-Franco, 1995) | 76.050 | 1.975 | - | Ver Figura 6a | - | - | 1.148 |
| Autor * | 74.814 | 2.126 | - | - | 2.626 | 0.809 | 2.829 |
| Autor** | 67.208 | 2.144 | - | - | Ver Figura 6b | 0.788 | 0.989 |
Figura 6 Variabilidad temporal de los parámetros de decrecimiento/agotamiento: a) Modelo de Forkasiewicz y Paloc; b) modelo propuesto por el autor.
Tabla 6 Resultados estadísticos para los experimentos numéricos.
| Parámetros | Modelos de descarga de acuíferos | |||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Iglesias (1984) | Tisson et al (Pérez-Franco, 1995) | Forkasiewicz-Paloc * (Pérez-Franco, 1995) | Forkasiewicz-Paloc ** (Pérez-Franco, 1995) | Autor * | Autor ** | |||||||||||||
| μθ | σθ | Cvθ | μθ | σθ | Cvθ | μθ | σθ | Cvθ | μθ | σθ | Cvθ | μθ | σθ | Cvθ | μθ | σθ | Cvθ | |
| m | 48.633 | 242.113 | 0.320 | 42.736 | 283.959 | 0.394 | 68.020 | 19.275 | 0.065 | 64.947 | 101.228 | 0.155 | 57.183 | 184.134 | 0.237 | 51.423 | 253.964 | 0.310 |
| n | 2.212 | 0.004 | 0.029 | 2.337 | 0.007 | 0.037 | 1.817 | 1.879 E-04 | 0.008 | 1.975 | 0.001 | 0.020 | 2.180 | 0.002 | 0.021 | 2.200 | 0.006 | 0.034 |
| α | 2.853 E-07 | 4.884 E-18 | 0.008 | 1.727 E-07 | 4.052 E-18 | 0.012 | - | - | - | - | - | - | - | - | - | - | - | - |
| β | - | - | - | - | - | - | 5.751 E-06 | 1.019 E-13 | 0.056 | ver figura 7d | ver figura 7d | ver figura 7d | - | - | - | - | - | - |
| 𝜅 | - | - | - | - | - | - | - | - | - | - | - | - | 2.576 E-07 | 6.166 E-17 | 0.030 | ver figura 7f | ver figura 7f | ver figura 7f |
| 𝜂 | - | - | - | - | - | - | - | - | - | - | - | - | 0.803 | 0.001 | 0.037 | 0.787 | 6.707 E-05 | 1.041 E-02 |
Figura 7 Comparación de los caudales simulados y observados: a) Iglesias (1984); b) Tisson et al. (Pérez-Franco, 1995); c) Forkasiewicz y Paloc* (Pérez-Franco, 1995), con parámetro de decrecimiento β constante; d) Forkasiewicz y Paloc** (Pérez-Franco, 1995), con parámetro de decrecimiento β variable; e) modelo no lineal propuesto por el autor*, con parámetro de agotamiento 𝜅 constante, f) modelo no lineal propuesto por el autor**, con parámetro de agotamiento 𝜅 variable.
Tal y como se muestra en la Figura 8, el modelo basado en la ecuación de agotamiento dada por Forkasiewicz y Paloc presenta variaciones notables del parámetro de decrecimiento de caudal β, lo que redunda en una pobre estimación de éste en el experimento conducido, a pesar de haber mostrado el mejor desempeño en una de las réplicas de las 20 ejecutadas. El modelo no lineal presentado por el autor fue muy estable en cualquiera de sus alternativas (𝜅 constante o variable) y competitivo, en comparación con el modelo de Iglesias (1984) y Tisson et al. (Pérez-Franco, 1995).
Figura 8 Varianza de los parámetros de decrecimiento/agotamiento: a) en el modelo de Forkasiewicz y Paloc (Pérez-Franco, 1995); b) en el modelo propuesto.
A pesar de disponer de una sola estación pluviométrica en este estudio para validar cada uno de los modelos tratados, es posible llevar a cabo simulaciones con cada uno de ellos, con el objetivo de restituir los datos históricos dentro del periodo enero de 1929 a diciembre de 2010 (serie de 82 años), siendo válido destacar que el modelo de Forkasiewicz y Paloc (Pérez-Franco, 1995) se ha excluido en cualquiera de sus variantes. En la Figura 8 puede apreciarse una adecuada convergencia entre ellos para el régimen de recarga por precipitación a que fue sometido el acuífero Ciego-Morón en el periodo de referencia. Cabe señalar que aunque se observan descargas máximas mensuales que pueden llegar a superar los 10 m3/s, la tendencia de éstas es decreciente.
Conclusiones
En el presente artículo diferentes modelos agregados para estimar las descargas de acuíferos se implementaron en el asistente matemático MatLab y acoplados al algoritmo CS-FMIN, demostrándose la robustez de la estrategia abordada para reproducir series de caudales de descarga en acuíferos kársticos. En cada uno de los modelos se identificaron los parámetros relacionados con el volumen de precipitación drenado por el acuífero y aquellos que caracterizan el agotamiento.
Fue desarrollado por el autor un modelo de descarga no lineal y su desempeño demostró ser satisfactorio en los experimentos numéricos conducidos, al ser contrastado con las observaciones en la estación La Yana, y con las formulaciones de Iglesias (1984) y Tisson et al. (Pérez-Franco, 1995), así como de Forkasiewicz y Paloc (Pérez-Franco, 1995), siendo este último desarrollado por sus autores para medios kársticos, como el caso aquí abordado. Sin embargo, la demostrada variabilidad temporal del parámetro de decrecimiento del caudal limita esta formulación sólo al completamiento de series de caudales de descarga dentro de intervalos de tiempo específicos y no restituir datos fuera de éstos. En este sentido, la propuesta no lineal mostró gran adaptabilidad, al considerar parámetros constantes o variables en el tiempo, así como las correlaciones a las series observadas, que fueron elevadas.
