Área de estudio

Seis subcuencas hidrográficas se seleccionaron entre 3° 30’ y 5° 8’ S y 79° 10’ y 80° 29’ O, que pertenecen a las cuencas binacionales Catamayo-Chira y la Puyango-Tumbes ubicadas en la frontera entre Ecuador y Perú.

La cuenca Catamayo Chira tiene 17 199.19 km² y 817 968 habitantes. Esta cuenca nace en las cumbres de la divisoria continental de aguas (Ecuador) y desemboca en el Océano Pacifico (Perú), atraviesa escenarios montañosos y costaneros, con climas tropicales, ecosistemas diversos, y varios usos y sistemas administrativos. La geografía de la cuenca es abrupta y con altitudes de 0 y 3700 msnm y tiene 11 zonas de vida, desde el desierto tropical hasta el bosque pluvial montano. La precipitación media anual es 800 mm y varía desde 10 mm en la zona baja, hasta 1000 mm en su cabecera. Esta cuenca está cubierta en 14% por bosque húmedo, 41% por bosque seco, 30% por pastizales, 10% por cultivos y 5% por otros usos.

La cuenca Puyango Tumbes tiene 5 400 km² y comprende las provincias de El Oro y Loja en Ecuador y el departamento de Tumbes en Perú. La cuenca alta del río Puyango Tumbes está formada por terrenos montañosos, con fuertes pendientes y altitudes desde 500 m hasta 3 700 m, y presenta zonas de páramo, bosque natural y cultivos. La cuenca alta se caracteriza por una actividad minera importante, en especial en las cuencas de los ríos Calera y Amarillo. La cuenca baja tiene amplias llanuras y es cultivada con intensidad, con plantaciones de arroz; además hay zonas grandes de bosque seco, pastizales y cultivos agrícolas. La cuenca baja en la zona peruana sufre inundaciones con frecuencia. La precipitación en la cuenca varia de 200 mm a 1150 mm anuales con temperaturas medias de 13 °C a 25 °C La cuenca está formada por tres tributarios principales: el río Calera, el río Amarillo y el río Pindo, y el volumen total anual de agua que produce la cuenca (masa anual promedio) es cerca de 3 400 millones m³.

La disponibilidad de información hidrométrica y características geográficas determinaron la selección de seis subcuencas: las cuencas de los ríos Amarillo y Pindo pertenecientes a la cuenca Puyando-Tumbes, y las cuencas de los ríos Catamayo, Macará, Quiroz y Alamor de la cuenca Catamayo-Chira (Cuadro 1). La ubicación geográfica de las cuencas se muestra en la Figura 1.

Modelo de Témez

El modelo de Témez (^{Témez, 1977}, ^{Estrela Monreal, 1999}) es un modelo hidrológico concentrado que supone que el perfil del suelo está dividido en una zona superior insaturada y otra inferior saturada, y su comportamiento semeja un embalse subterráneo que desagua en la red superficial.

La lluvia (P) se divide en evapotranspiración (ET) y excedente (T). El excedente se descompone en una parte que discurre en la superficie, E (escorrentía superficial), y otra que se infiltra hasta el acuífero (I). La primera evacua por el cauce dentro del período de tiempo considerado, mientras que la segunda permanece en los reservorios subterráneos y se libera después. El excedente se calcula mediante la siguiente expresión:

en la cual:

P _i es la precipitación en el periodo desde el instante i-1 al instante i (mm); P _o es el umbral de escorrentía que define la altura de precipitación, por debajo de la cual no se produce escorrentía (mm); Ti es el excedente en el periodo desde el instante i-1 al instante i (mm); H _max es la capacidad máxima de humedad en el suelo (mm); H _i-t la humedad en el suelo en el instante i-1 (mm); EP _t es la evapotranspiración potencial desde el instante i-1 (mm) y; C es un parámetro del modelo.

La humedad en el suelo, H _i , al final del período resultará:

Una evapotranspiración real Ei (mm) se produce y es igual a:

La ecuación 6 muestra que se puede evapotranspirar toda el agua disponible con la evapotranspiración potencial como el límite superior.

El modelo adopta una cantidad de infiltración (I _i ) como una función del excedente Ti y del parámetro de infiltración máxima I _máx .

will be in mm.

La infiltración aumenta con el excedente, es asintótica para sus valores altos y tiene como límite el valor I _máx .

Esta infiltración, I _i , se convierte en recarga al acuífero, Ri, en tanto que el resto del excedente (E _i = T _i I _i ) será agua superficial. El modelo asume que el tiempo de paso por la zona no saturada es menor al tiempo de paso de la simulación.

El drenaje del acuífero es modelado por una función exponencial del tipo siguiente:

donde Q _i es la descarga en el instante i, a es el coeficiente de la rama de descarga del acuífero, que depende de las condiciones particulares de la cuenca en estudio y t el intervalo de tiempo entre los instantes i-1 e i.

La relación entre el caudal de descarga, Q _i , y el volumen, V _i , almacenado en el acuífero es igual a

La descarga por la infiltración se supone concentrada en la mitad del período, por lo cual la ley de los caudales subterráneos es:

en la cual R _i la recarga al acuífero en el periodo i-1 a i, coincidente con la infiltración I _i .

La aportación subterránea a lo largo del período, A _SUBi , es:

La aportación total, A _Ti , es dada por la escorrentía superficial (T _i I) y la aportación subterránea.

Los parámetros del modelo son cuatro: Hmax, la capacidad máxima de humedad del suelo; Imax, la capacidad máxima de infiltración; C el parámetro de excedente; a es el coeficiente de la rama de descarga. Todos estos parámetros son particulares para cada cuenca y están sujetos a calibración.

Soil and water assessment tool (SWAT; Herramienta de evaluación de agua y suelo)

SWAT es un modelo continuo de simulación usado para pronosticar el impacto de las practicas de manejo del suelo en el agua, sedimentos y nutrientes en una cuenca (^{Di Luzio et al., 2002}). Además, es un modelo continuo semi-distribuido basado en la ecuación de balance hídrico en el perfil del suelo y simula los procesos de precipitación, infiltración escorrentía superficial, evapotranspiración, flujo lateral y percolación. La ecuación de balance hídrico es:

donde SW _t es el contenido de humedad final en el suelo (mm H2O), SW es el contenido de humedad inicial en el suelo en el día i (mm H2O), t es el tiempo (días), Rday es la cantidad de precipitación en el día i (mm H2O), Qsurf es la cantidad de escorrentía superficial en el día i (mm H2O), Ea es el nivel de evapotranspiración en el día i (mm H2O), wseep es la cantidad de agua en la zona de saturación del perfil del suelo en el día i (mm H2O), y Qgw es el flujo de retorno en el día i (mm H2O).

El modelo SWAT realiza una división topográfica de la cuenca en subcuencas en base a un umbral de área, para después subdividir las subcuencas en una o varias unidades de respuesta hidrológicas (HRU) que representan la combinación única de tipo y uso del suelo. La respuesta de cada HRU en términos de agua, sedimentos, nutrientes y pesticidas son determinadas individualmente para luego ser agregadas a nivel de subcuenca y movidas a la salida de la cuenca a través de su red hídrica (^{Bouraoui et al., 2005}).

La escorrentía superficial es estimada con datos diarios de precipitación usando la metodología del número de la curva o la ecuación de Green-Ampt (^{Neitsch et al., 2002}). La evapotranspiración se determina aplicando las metodologías propuestas por Hargreaves, Priestley-Taylor o Penman-Monteith (^{Neitsch et al., 2002}). Un reservorio cinemático que considera variaciones en conductividad hidráulica, pendiente y humedad del suelo, es usado para predecir el flujo lateral en cada una de sus capas. El agua subterránea se divide en dos sistemas de acuíferos: uno no confinado, poco profundo, que contribuye al flujo de retorno; y uno confinado, profundo, desconectado del sistema a menos que se considere agua subterranea (^{Bouraoui et al., 2005}). La tasa de sedimentación es estimada aplicando la ecuación universal de pérdida de suelo modificada (^{Neitsch et al., 2002}), para lo cual se usa la escorrentía superficial, la tasa pico de flujo, la erodabilidad del suelo, la longitud de la pendiente, su inclinación, el factor de cultivo y las prácticas de manejo en la zona. ^{Neitsch et al. (2002)} explicaron las bases teóricas del modelo SWAT.

Implementación, calibración y validación del modelo

Para cada cuenca, la modelación concentrada fue realizada considerando como datos de entrada la precipitación media en la cuenca y la evapotranspiración media. Para fines de calibración se consideró los caudales medios mensuales registrados en la estación de cierre de cada una de las cuencas incluidas en el Cuadro 1. Los datos de precipitación en 24 estaciones ecuatorianas y en 19 estaciones peruanas, así como datos de temperatura en 10 estaciones ecuatorianas y 4 peruanas fueron considerados. La correspondencia de las series de datos fue verificada mediante análisis de curvas de doble masa. Mediante un análisis de regresión entre estaciones que presentaban vecindad geográfica, mismo régimen climático y correspondencia de sus registros, los datos se homogenizaron al período común 1970-2000 en las estaciones de la cuenca Catamayo-Chira, y 1965-1995 en las estaciones de la cuenca Puyango Tumbes. Con el método de Thornthwaite (^{Thornthwaite y Mather, 1957}) se calculó la evapotranspiración potencial en cada una de las 14 estaciones con registros de datos de temperatura.

La precipitación y la evapotranspiración media para cada una de las cuencas analizadas se calculó con base en polígonos de Thiessen (Figura 2).

En cada una de las cuencas seleccionadas se realizó la calibración de los cuatro parámetros del modelo de Témez con base (C, Hma, Imax, a) en aproximaciones sucesivas. Para esta calibración se usaron 70 % de los registros recopilados; los caudales calculados por el modelo se calcularon con los caudales registrados en las estaciones incluidas en el Cuadro 1. Una vez calibrados los parámetros del modelo se realizó su validación generando valores de caudal para un período correspondiente al restante 30 % del registro histórico. La bondad de ajuste entre los valores calculados y observados se midió calculando el coeficiente de correlación (R), la eficiencia de Nash Sutcliffe (EF) (^{Nash y Sutcliffe, 1970}) y el error cuadrático medio (RMSE). El coeficiente de correlación medirá el ajuste de los datos observados y calculado a una línea recta, la eficiencia EF medirá la relación 1 a 1 de los datos observados y calculados, y el RMSE permitirá cuantificar las diferencias medias entre los valores observados y calculados.

La aplicación del modelo de Témez se realizó usando el software CHAC, desarrollado por el CEDEX de España. El modelo Temez se puede encontrar en http://hercules.cedex.es/Chac/.

Con fines de comparación se implementó el modelo SWAT en cada una de las subcuencas. Un DEM SRTM del área de estudio se empleó para la delimitación de cuencas de aportación y el cálculo de los diferentes parámetros morfométricos requeridos por el modelo. La información respecto al tipo de suelo se reunió y los parámetros del suelo se determinaron con pruebas de laboratorio. El uso del suelo se obtuvo de una imagen LANDSAT7 ETM+ del 2 de octubre de 2011, mediante una clasificación supervisada aplicando el criterio del mínimo ángulo espectral (^{Richards y Jia, 2006}; ^{Chuvieco, 2002}). Los parámetros de uso del suelo se tomaron de la base de datos del modelo SWAT, con modificaciones menores. Un mapa de tipo de suelo (^{Valarezo, 2007}) fue validado y adaptado para implementar el modelo SWAT y esta validación se realizó mediante perforaciones y pozos en cada una de las unidades edafológicas identificadas. La prueba de infiltración se realizó in situ. Muestras inalteradas de suelo fueron obtenidas y las características generales de los horizontes de los suelos fueron registradas. Con todos los datos y los resultados de los análisis de laboratorio, los parámetros del tipo de suelo fueron determinados para cada unidad edafológica identificada.

La combinación del uso y tipo de suelo permitió determinar unidades de respuesta hidrológica (UHR). La información de tipo de suelo se utilizó para caracterizar el acuífero no confinado. Además, datos de precipitación y temperatura diaria, así como de precipitación, temperatura, radiación solar, velocidad del viento, humedad relativa mensual para cada una de las estaciones fue usada para calcular la evapotranspiración empleando el método de Penmman-Monteith. Los datos faltantes de precipitación fueron estimados mediante un análisis de correlación ortogonal entre estaciones vecinas geográficamente y con similares condiciones climáticas. En el caso de la temperatura, una ecuación de regresión entre elevación-temperatura fue determinada. Un proceso similar se aplicó para las otras variables.

Antes de la calibración, un análisis de sensibilidad fue realizado mediante el cambio sistemático de los valores de los parámetros del modelo y observando el efecto dichos cambios en los resultados obtenidos, seleccionando un número reducido de parámetros sensibles de importancia más alta en el cálculo y optimizándolos en la fase de calibración. El análisis de sensibilidad se realizó con el módulo que el modelo SWAT tiene para ello, y combina el método de muestreo del hipercubo Latino con la simulación de un factor a la vez (OAT) (van Griensven et al., 2006). El muestreo mediante el hipercubo Latino, y en contraposición al método tradicional de Monte Carlo, realiza un muestreo estratificado entre el rango de los valores posibles de cada parámetro; además, la simulación OAT asegura que los cambios en el modelo se pueden atribuir a los cambios realizados en la variable de entrada en cada simulación (van Griensven et al., 2006).

El Modelo SWAT fue calibrado manualmente ajustando los valores de los parámetros más sensibles cambiando su valor inicial no más de 12% usando el 70% de los registros recopilados y optimizando R, EF y RMSE. Los parámetros del modelo SWAT relativos a la pérdida de suelo no fueron considerados en el proceso de calibración porque el modelo Temez no reproduce pérdida de suelo; así, la comparación entre los modelos solo se realizó con datos de caudal. La validación se realizó con el 30% restante de los registros. Después se efectuó una comparación gráfica de los caudales observados y calculados. En ^{Oñate-Valdivieso y Bosque (2014)} hay información adicional de la implementación del modelo SWAT en las cuencas en estudio.

Calibración de parámetros del modelo Témez

Los parámetros calibrados del modelo Témez para cada cuenca en estudio se presentan en el Cuadro 2.

La humedad máxima (H _máx ) permite definir el umbral de escorrentía (P _o ), es decir la altura de precipitación por debajo de la cual no se produce escorrentía y así definir el excedente (T _i ) compuesto por la escorrentía superficial y la infiltración hacia los acuíferos. En el Cuadro 2 se muestra que los valores de H _máx superan los 150 mm en la mayoría de los casos, con un valor máximo de 400 mm. Los valores que alcanza este parámetro evidencian una alta capacidad de almacenamiento de los suelos de la zona de estudio, los cuales están sometidos a un considerable déficit hídrico la mayor parte del año. Esta situación puede encontrarse en la subcuenca del río Alamor; en cambio, lo contrario se presenta en la subcuenca del río Amarillo que es árida en su parte baja pero presenta niveles mayores de precipitación durante el año, por lo cual los suelos presentan mejores condiciones de humedad y una menor capacidad máxima de infiltración.

El parámetro de excedencia C se usa para definir el umbral de escorrentía P _o actuando como un factor de ponderación de la diferencia entre la humedad máxima del suelo y su humedad en un instante determinado. Los valores del parámetro C incluidos (Cuadro 2) son menores a la unidad, controlando posibles subestimaciones del excedente producidas al considerar factores de humedad relativamente altos.

La infiltración real I es directamente proporcional a la infiltración máxima I _máx , y es una fracción de ésta. Al observar los valores de este parámetro (Cuadro 2) se observan importantes niveles de infiltración máxima, lo cual concuerda con los valores de H _máx . Al ser una zona árida es de esperarse que existan grandes niveles de infiltración.

El coeficiente α alcanza valores medios entre 0.01 y 0.08. Este parámetro se relaciona con las aportaciones de agua subterránea al flujo base de cada una de las cuencas en análisis. Los valores mayores se presentan en las cuencas de mayor superficie que por su extensión abarcan zonas en las que hay aporte de aguas subterráneas de diferente magnitud. En las partes altas, generalmente más húmedas y en mejores condiciones de conservación vegetal, el aporte de las aguas subterráneas es mayor, a diferencia de las partes bajas donde ocurre lo contrario. El parámetro calibrado en las cuencas de mayor tamaño es representativo de dichas condiciones.

Evaluación del modelo Témez

Los parámetros usados para medir la bondad de ajuste entre los caudales observados y calculados para el período de validación se incluyen en el Cuadro 3. La relación gráfica entre los caudales observados y calculados se presenta en la Figura 3.

En las cuencas de menor tamaño (Amarillo, Pindo y Alamor) el modelo reproduce los caudales de la temporada seca de mejor manera que en la temporada lluviosa. Al aumentar el tamaño de la cuenca se observa una sobreestimación de los caudales llegando casi a desaparecer la correspondencia entre los caudales observados y calculados. La sobreestimación de los caudales se debe a que el caudal se calcula como el producto de la aportación (mm/mes) por el área de la cuenca, considerando así que los niveles de precipitación y evapotranspiración son uniformes en toda la cuenca. Esto no se cumple necesariamente en cuencas de tamaño medio o grande.

La cuenca del río Amarillo presenta un coeficiente de correlación de 0.64 y una eficiencia EF de 0.63 (Cuadro 3), valores que confirman la existencia de un ajuste aceptable entre los valores observados y calculados. El error cuadrático medio es 7.32 m³ s^-1 el cual, comparado con el caudal promedio registrado (14.32 m³ s^-1), es excesivo. En la Figura 3 se muestra que el modelo no reproduce con fidelidad los caudales picos y los caudales de la temporada seca son sobreestimados.

En la cuenca del río Pindo, aproximadamente el doble del tamaño de la del río Amarillo, hay un coeficiente de correlación entre los datos observados y calculados de 0.63 que es apenas aceptable, con una EF de 0.49 inferior a la observada en la cuenca del río Amarillo. El error cuadrático medio es 14.47 m³ s^-1, que comparado con el caudal medio en el río Pindo (23.43 m³ s^-1) resulta excesivo. El modelo sobreestima los caudales y se observa en la Figura 3 que el mejor ajuste se produce en los caudales de la temporada seca, y los caudales picos muestran el peor ajuste.

En la cuenca del río Alamor (Cuadro 3) se observa una correlación relativamente buena entre los valores calculados y observados (R=0.72), aunque hay un valor bajo de la EF (0.35), con un error cuadrático medio de 9.11 m³ s^-1 que supera al caudal medio observado que es 7.3 m³ s^-1. Podría esperarse que con un coeficiente alto de correlación el RMSE debería ser menor, pero al observar la Figura 3 se comprueba que la tendencia es bien reproducida por el modelo, de ahí el valor satisfactorio del coeficiente de correlación. Además los caudales pico registrados son sobreestimados por el modelo (lo cual influye en el RMSE). La EF tiene un valor un poco más realista, con 0.35 que no es óptimo (EF=1), pero da una idea de que el modelo tiene un desempeño hasta cierto punto aceptable.

La cuenca del río Catamayo (Cuadro 3) presenta los peores resultados y no existe correlación entre los caudales observados y calculados. Al obtenerse un valor negativo de la eficiencia de Nash-Sutcliffe se considera que la media de los valores observados es mejor predictor que el modelo en sí.

La cuenca del río Quiroz presenta valores aceptables del coeficiente de correlación (R=0.69) y de la eficiencia de Nash-Sutcliffe (EF=0.62), aunque el excesivo error cuadrático medio, que duplica el caudal medio observado, evidencia la sobreestimación de los caudales calculados.

La cuenca del río Macará presenta el mejor coeficiente de correlación (R=0.82), pero el valor bajo de la eficiencia de Nash-Sutcliffe (EF=0.27) implica que el modelo no puede reproducir el comportamiento hidrológico de la cuenca. Además, el error cuadrático medio casi duplica el caudal medio observado por sobreestimar los caudales en general.

La zona de estudio es principalmente montañosa, con una distribución irregular de la precipitación y en promedio solo hay una estación pluviométrica cada 270 km², y se presume que la red de estaciones meteorológicas en la zona no posee la densidad suficiente para reflejar con precisión la variación espacial de la precipitación. Además, los registros históricos de precipitación, temperatura y sobre todo de caudal, presentaron vacíos notables de información que, hasta cierto punto, evidencian problemas en el manejo de las estaciones meteorológicas e hidrológicas en la zona.

Al igual que la precipitación, la pendiente, el tipo de suelo y su uso presentan una significativa variación espacial. Por lo tanto, considerar un valor único, concentrado, representativo de cada condición puede ser válido en cuencas muy uniformes y de dimensiones reducidas.

Comparación con el modelo SWAT

La obtención de datos meteorológicos en la zona de estudio fue difícil porque el número de estaciones climáticas es reducido, solamente 14 en el área de estudio. La información sobre velocidad del viento, radiación solar y humedad relativa fue aún más escasa.

El análisis de sensibilidad reportó los 35 parámetros más sensibles los cuales fueron clasificados en una escala de 35 niveles en la que 1 corresponde al parámetro más sensitivo sensible. El Cuadro 4 muestra los 10 primeros parámetros reportados por dicho análisis y su nivel de influencia en el modelo.

En el Cuadro 4 se observa que el modelo fue muy sensible al número de la curva CN2, porque este parámetro permite separar la precipitación efectiva de la precipitación total, y por ende tiene un efecto directo en el cálculo de la escorrentía. Además, el contenido de humedad del suelo SOLAWC condiciona el escurrimiento superficial porque causa mayor o menor infiltración. La constante de recesión del flujo base ALPHA_BF regula las aportaciones de las aguas subterráneas al caudal. El albedo del suelo húmedo SOL_ALB y el coeficiente de compensación de la evaporación del suelo ESCO regularon las pérdidas por evapotranspiración. Una influencia importante presentaron la pendiente media SLOPE y el de retardo de la escorrentía superficial SURLAG, afectando directamente el tiempo de concentración de cada subcuenca y variando la ocurrencia temporal de los caudales picos. Los parámetros de pérdida de suelo no fueron considerados.

Un análisis gráfico de la relación entre los caudales observados y calculados con el modelo Témez y el modelo SWAT se presenta en la Figura 4. Allí se observa que el modelo SWAT presenta una mejor representación de las tendencias y un ajuste mayor con los valores observados, incluso si se consideran los caudales de crecida de la temporada húmeda, así como los de estiaje en la temporada seca.

El análisis comparativo del desempeño del modelo Témez con el modelo SWAT se resume en el Cuadro 5. Allí se observa que en todos los casos el modelo SWAT presenta un mejor coeficiente de correlación y en la mayoría de los casos alcanza valores superiores al 80%. La eficiencia de Nash-Sutcliffe es significativamente al calcular con los valores obtenidos con el modelo SWAT. De manera similar, el error cuadrático medio fue significativamente menor en los valores calculados con el modelo SWAT.

El SWAT es un modelo semidistribuido y se podía esperar un mejor desempeño en todas las cuencas, pero aún en las cuencas más pequeñas el desempeño del modelo Témez es menor al del modelo SWAT, lo cual podría aceptarse dada la menor cantidad de información requerida para implementar.

Si se considera que en muchos lugares la carencia de información es una realidad, el modelo Témez puede ser una opción viable en cuencas relativamente pequeñas. Además, el modelo Témez se podría usar como base de un modelo semidistribuido aplicable a cuencas más grandes.