Descripción del Modelo

El software en código abierto Soil and Water Assessment Tool (SWAT) es un modelo lluvia-escorrentía, de parámetros semidistribuidos, capaz de simular numerosos procesos físicos a escala temporal continua (anual, mensual, diaria y subdiaria). Su objetivo principal es predecir la producción de agua y sedimentos en cuencas hidrográficas, así como el impacto de prácticas de gestión agrícola sobre la calidad del agua (nutrientes y pesticidas). Presenta una precisión razonable en cuencas de gran tamaño, con variedad en el relieve, en tipos y usos del suelo. Su elevada resolución espacial permite que sea implementado tanto a escala continental como de cuenca hidrológica (^{Arnold et al., 1998}, ²⁰¹²). No obstante, el uso de SWAT a gran escala requiere simplificaciones importantes. ^{Gassman et al. (2007)} y ^{Abbaspour et al. (2015)} muestran múltiples aplicaciones de SWAT a escala continental y de cuenca en Europa, África, Asia y EEUU. La componente hidrológica de SWAT permite calcular los elementos del balance hídrico y, en consecuencia, los recursos hídricos (agua azul, verde, etc.) a nivel incluso de subcuenca. La fase terrestre del ciclo hidrológico es simulada con base en la siguiente ecuación del balance hídrico (Ecuación 1).

Donde SW _t y SW ₀ son el contenido final e inicial de agua en el suelo, R _day es la precipitación, Q _surf es la escorrentía superficial, E _a la evapotranspiración, w _seep el agua que percola en la zona vadosa desde el perfil del suelo, y Q _gw es el flujo de retorno. Las variables se expresan en mm de H₂O para el día, i, y el tiempo, t, en días (^{Neitsch et al., 2011}).

En SWAT, cada cuenca hidrológica es dividida en subcuencas a partir de las coberturas disponibles para el relieve del terreno y los cursos fluviales. A su vez, cada subcuenca se subdivide en polígonos, irregulares y de tamaños variables, denominados Unidades de Respuesta Hidrológica (HRU por sus siglas en inglés). Las HRU resultan de la superposición de un mismo tipo de suelo, uso de suelo y pendiente en una subcuenca dada. Suponen el principal avance de SWAT con respecto a otros modelos hidrológicos similares (^{Winchell et al., 2013}). Dicha ventaja permite que, en SWAT, la simulación hidrológica en la fase terrestre no se realice únicamente a nivel de subcuenca sino también a escala de HRU (^{Arnold et al., 2012}).

Al margen de la componente hidrológica, el desarrollo actual del modelo permite que sea utilizado en: la toma de decisiones frente a alternativas de gestión de recursos hídricos; evaluación de contaminación por fuentes no puntuales incluso en grandes cuencas fluviales; predecir el impacto agrícola sobre sedimentos, compuestos químicos y recursos hídricos de cuencas con diferentes condiciones de gestión, tipos y usos de suelo durante largos períodos de tiempo (^{Neitsch et al., 2011}; ^{Abbaspour et al., 2015}).

Datos de entrada en SWAT

El modelo fue implementado como una extensión del software ArcGIS 10.2.2 denominada ArcSWAT 2012.10.2.18 (^{Arnold et al., 2012}; ^{Winchell et al., 2013}), en adelante SWAT. La delimitación de la cuenca alta del río Taibilla, y de las subcuencas de éste, se realizó a partir del Modelo Digital de Elevación del Terreno (MDE), con resolución de 25 m x 25 m, propiedad del Centro Nacional de Información Geográfica (^{CNIG, 2017}). Junto con el MDE, se agregó a SWAT la cobertura de la red hidrográfica disponible en la Confederación Hidrográfica del Segura (^{CHS, 2017}). La salida de la cuenca se ubicó en la estación de aforo 7017 localizada en el embalse del Taibilla. A partir de lo anterior, SWAT delimitó la cuenca alta del Taibilla con 315 km² distribuidos en 25 subcuencas (Figura 3).

Figura 3 Coberturas del área de estudio en SWAT: relieve (a); usos de suelo (b); edafología (c); clases de pendiente (d). 

En la elaboración de las HRU se empleó el Mapa Nacional de Usos del Suelo del Proyecto Corine Land Cover (CLC2006) disponible en ^{CNIG (2017)} a escala 1:100000. Los códigos de uso de CLC2006 fueron reclasificados a los de SWAT utilizando la metodología descrita por ^{Bobal et al. (2010)} y ^{El-Sadek e İrvem (2014)}. La edafología de la cuenca se extrajo de la base de datos Harmonized World Soil Database (^{HWSD, 2012}, ²⁰¹⁷), a escala 1:1000000, que recoge las unidades de suelo (códigos) establecidas en ^{FAO (1990)} y ^{WRB (2014)}. A partir de la metodología de ^{Gomariz-Castillo (2015)}, los códigos de suelo de ^{HWSD (2017)} fueron incluidos en la base de datos de SWAT. La pendiente fue dividida en cinco clases (0%-10%, 10%-20%, 20%-30%, 30%-40% y >40%). Dichos umbrales se fijaron ya que la cuenca presenta tramos altos, de cabecera o de alta montaña, con pronunciadas pendientes en las laderas vertientes al cauce del río (^{AYESA, 2010a}; ^{Quiñonero-Rubio et al., 2013}). Los terrenos con pendientes suaves quedan relegados casi por completo al cauce del río Taibilla y sus arroyos (Figura 3).

A partir de las coberturas mostradas (Figura 3) se obtuvieron 769 HRU. Debido al gran número de HRU generadas, SWAT se comporta como un modelo de parámetros distribuidos. El tamaño de cada HRU resulta por tanto menor que la superficie de las celdas en las que se divide el territorio, situada habitualmente entre 12 km² y 5 km² (^{Cabezas, 2015}). Por ello, los parámetros hidrológicos que caracterizan a la cuenca varían espacialmente en cada punto (HRU < celda).

Tras elaborar las HRU, se incorporaron en SWAT las series diarias de precipitación, temperatura máxima y mínima de las estaciones meteorológicas localizadas en las inmediaciones y en el interior de la cuenca (periodo 01/01/1944 - 30/04/2016). Los citados registros fueron facilitados por AEMET, y posteriormente tratados en el software R3.3.2 (R, ²⁰¹³) para incluir valores ausentes, años bisiestos, etc. A partir de estas series, se seleccionó el método de Hargreaves (^{Hargreaves et al., 1985}) para calcular la Evapotranspiración Potencial (ETP). Por otro lado, SWAT asigna al centroide de cada subcuenca, calculado previamente, los registros de la estación más cercana (^{Winchell et al., 2013}). La distribución de las ocho estaciones escogidas en la cuenca resulta homogénea (Figura 1), permitiendo con ello que el modelo realice una asignación diferenciada y adecuada, de las variables precipitación y temperatura, a los distintos grupos de subcuencas. Gracias a la distribución que presentan dichas estaciones, no se incluyeron en SWAT las denominadas bandas de elevación con las que aplicar un gradiente a las series de precipitación y temperatura en función de la altitud (^{Winchell et al., 2013}).

Calibración, validación y análisis de sensibilidad de SWAT

Los resultados (simulaciones) de los modelos hidrológicos deben ser corregidos con datos medidos. Se restringe así la incertidumbre de las simulaciones pudiendo representar, aunque simplificado, el comportamiento de una cuenca hidrográfica (^{Gupta et al., 2009}; ^{Abbaspour et al., 2015}). Este procedimiento se realiza con las técnicas de calibración, validación y análisis de sensibilidad (^{Moriasi et al., 2007}). La calibración evalúa los parámetros del modelo comparando sus predicciones, para un conjunto de condiciones de contorno asumidas, con datos observados para las mismas condiciones. En la validación se lanza el modelo utilizando parámetros determinados en calibración, de modo que éste sea capaz de realizar simulaciones próximas a los datos reales de la serie temporal destinada a validar el modelo (^{Krause et al., 2005}; ^{Moriasi et al., 2007}). En SWAT, la calibración y validación, pueden realizarse con medidas de caudal, nitratos, fosfatos, etc. (^{Abbaspour, 2015}; ^{Abbaspour et al., 2015}). En el presente estudio, se escoge la variable caudal a la salida de las subcuencas para calibrar y validar SWAT. El análisis de sensibilidad determina el cambio en las salidas del modelo con respecto a las entradas (parámetros) de éste. Es un proceso necesario para identificar parámetros clave y la precisión de los parámetros requerida en calibración (^{Ma et al., 2000}; ^{Moriasi et al., 2007}).

Descripción del algoritmo SUFI-2 en SWAT-CUP

La optimización (calibración) y el análisis de sensibilidad de SWAT se realizaron con el algoritmo SUFI-2 (Sequential Uncertainty Fitting version-2) incluido en el software SWAT-CUP2012 (^{Abbaspour, 2015}). En SUFI-2 la calibración se realiza con una serie de iteraciones que incluyen numerosas simulaciones. Cada nueva iteración es alimentada con los resultados de la anterior. Esto persigue aproximar (optimizar) la variable simulada, caudal a la salida de las subcuencas, a la observada. Los resultados de las iteraciones son un conjunto de valores (rangos) asignados a los parámetros que representan los procesos hidrológicos, las características físicas y la dinámica de cada cuenca hidrográfica. Cada nueva iteración presenta intervalos (rangos) de los parámetros recursivamente más próximos a su valor real. Ello pretende acotar la incertidumbre existente en los rangos iniciales de los parámetros, dado que a menudo no se dispone de mediciones de estos (^{Abbaspour, 2015}; ^{Abbaspour et al., 2007}). Así, a partir de mediciones de caudales, introducidas en SUFI-2, éste proporciona iteración tras iteración mayor exactitud en los rangos de los parámetros de cada área de estudio. Dicho procedimiento se denomina Modelización Hidrológica Inversa (^{Beven y Binley, 1992}; ^{Abbaspour et al., 2004}, ²⁰⁰⁷). En el estudio se escogen 14 parámetros hidrológicos condicionados por los métodos de cambio R-relative y V-replaces (^{Abbaspour, 2015}). Así mismo, en cada iteración, SUFI-2 dispone de estimadores (estadísticos) para evaluar la bondad del ajuste entre el caudal simulado (Qsim) y el observado (Qobs: Tabla 1, ecuaciones 2 a 6). En consecuencia, se escoge la simulación más próxima al caudal observado, es decir, la que presenta mejores resultados en dichos estimadores. Por ello, se consideran óptimos los valores de los parámetros obtenidos en dicha simulación. La validación por su parte se implementó en R3.3.2 para contrastar la serie de Qsim con Qobs mediante los citados estimadores.

Por otro lado, para iniciar la calibración, SUFI-2 requiere escoger un estadístico como función objetivo. La selección de esta función depende de los objetivos del estudio, de las variables a calibrar y de los parámetros hidrológicos a optimizar. Asimismo, la elección de la función afecta a los resultados de la calibración (parámetros hidrológicos, variable simulada-optimizada, etc. ^{Abbaspour et al., 2004}). En el estudio se han realizado cinco iteraciones (de 2000 simulaciones) con NSE y otras cinco con KGE. NSE es el estadístico más usado en modelación hidrológica como función objetivo dado que incrementa el rendimiento (ajuste) entre simulaciones y observaciones. Sin embargo, en ocasiones presenta errores por divergencia entre la serie simulada y la observada (^{Nash y Sutcliffe, 1970}). KGE posee un comportamiento similar a NSE, aunque proporciona un mejor ajuste del modelo en caso de disponer de mediciones con numerosos picos (caudales punta). Corrige así el sesgo generado por el modelo entre series, observadas y simuladas, con caudales punta (^{Gupta et al., 2009}; ^{Gomariz-Castillo, 2015}).

El análisis de sensibilidad, implementado en SUFI-2 tras la calibración, indica la relación entre la función objetivo y los parámetros hidrológicos seleccionados, cuyos valores cambian según presenten mayor o menor sensibilidad a la función. Los parámetros con mayor sensibilidad (significativa) poseen un P-valor ≤ 0.05. Un número elevado de parámetros con sensibilidad significativa indica una correcta elección de la función objetivo y de los parámetros. No obstante, a medida que las iteraciones son ejecutadas el número de parámetros significativos disminuye. Por ello, los parámetros significativos en las últimas iteraciones serán aquellos cuyo cambio afectará en mayor medida a los resultados del modelo (^{Abbaspour, 2015}; ^{Abbaspour et al., 2015}). Información detallada respecto a los procedimientos de calibración, validación y análisis de sensibilidad en el algoritmo SUFI-2 puede ser encontrada en ^{Abbaspour (2015)}.

Datos de entrada en SUFI-2

La cuenca alta del río Taibilla posee cuatro estaciones de aforo (Figura 1). En la calibración y validación de SWAT únicamente se consideraron los registros de caudales de la estación 7017 (periodo 01/01/1996-31/12/2012) facilitados por el CEDEX. Las restantes estaciones presentan series obsoletas (mediciones de 1930 a 1950) y con numerosos valores ausentes. La estación 7017, localizada en el embalse del Taibilla (salida de la cuenca), registra las reservas y salidas, a nivel mensual, del embalse. Dado que los embalses constituyen elementos de regulación de caudales, otorgando un régimen hídrico alterado, se realizó una operación de balance hídrico para restituir el régimen natural de los caudales en la salida de la cuenca. A partir de las reservas y salidas se obtuvieron las entradas medias mensuales en el embalse (Ecuación 7).

Donde E _n son las entradas al embalse en el mes actual (m³/s), R _n las reservas del embalse en el mes actual (m³/s), R _n-1 las reservas en el mes anterior (m³/s), S _n las salidas en el mes actual (m³/s), y Ev _n las pérdidas por evaporación en el mes actual (m³/s).

Las pérdidas por evaporación (m³/s) se calcularon a partir del área inundada del embalse (763009 m² variable según el año hidrológico) y de los registros de evaporación en el embalse (mm/mes) disponibles en ^{CHS (2015h)}. Dichas pérdidas afectan considerablemente al balance debido a la climatología del área de estudio, con clima semiárido, registrando en el embalse una evaporación media situada entre 52.5 mm/mes y 272.4 mm/mes (^{CHS, 2015h}). Las pérdidas del embalse por filtración se estiman despreciables debido a la formación del vaso del embalse compuesta por materiales impermeables (arcillas y margas) cuya acción se une al sedimento acumulado en el lecho del embalse (^{CHS, 2015c}).

Las entradas mensuales al embalse, obtenidas por balance (periodo 01/01/1996-31/12/2012), suponen la aportación total que éste recibe (superficial y subterránea). Los dos registros ausentes identificados en dicha serie fueron contrastados y completados con el aporte superficial del embalse durante el periodo 01/01/1996-31/12/2010 (^{AYESA, 2010b}) y con el caudal registrado en el Manantial Molino de las Fuentes (07/03/2007-07/04/2016) facilitado por la CHS.

Escenarios regionalizados de cambio climático en SWAT

Tras calibrar, validar y analizar la sensibilidad de SWAT, éste puede ser lanzado en escenarios pasados y futuros ya que los valores optimizados en los parámetros de la cuenca no varían con el tiempo (^{Winchell et al., 2013}). En coherencia con los objetivos del estudio, se incorporan en SWAT los escenarios de emisión RCP4.5 y RCP8.5 (periodo 2010-2050). Los modelos de circulación global (GCMs) seleccionados se denominan: bcc-csm1-1-m, MPI.ESM.MR y bcc.csm1.1. Se escogen estos modelos puesto que se trata de los GCMs más próximos al promedio anual de todos los modelos en la DHS (serie 2010-2050). Asimismo, con el fin de adecuar las predicciones de precipitación, temperatura máxima y mínima de los tres GCMs al área de estudio, se acoplan en SWAT las proyecciones de éstos a partir de la regionalización estadística realizada por ^{Petisco de Lara (2008a}, ^2008b) con el método de análogos (Tabla 2).

A partir de la regionalización por análogos y por regresión lineal, AEMET ha incorporado las series futuras, a nivel diario, de las citadas variables para el periodo 2006-2100, utilizando como referencia la serie 1961-2000, en prácticamente la totalidad de su red de estaciones meteorológicas en el territorio nacional (^{AEMET, 2017}). La red de estaciones de AEMET en la DHS, que incluye proyecciones de cambio climático, consta de 127 estaciones de precipitación y 39 de temperatura, de las cuales se consideraron en el estudio las indicadas en la Figura 1 dada su proximidad con la cuenca. Los ficheros a escala diaria de estas estaciones, para las variables, GCMs, y escenarios mencionados, fueron adaptados al formato requerido por SWAT (.txt). La operación se realizó automáticamente con un bucle ejecutado en el software bash4.3 en el que se incorporaron como valores ausentes (NA) los del 29 de febrero, caso de los años bisiestos, no incluidos en las proyecciones de AEMET. Asimismo, la tendencia lineal de las series anuales de precipitación y temperatura, resultado del promedio de los tres GCMs utilizados, se evaluó mediante una regresión lineal analizando la significancia de dicha tendencia con el contraste no paramétrico de Mann-Kendall (^{Kendall, 1976}; ^{Sneyers, 1992}).

Calibración, validación y análisis de sensibilidad con NSE y KGE

En la calibración, validación y análisis de sensibilidad de SWAT se contrastaron las mediciones del caudal, a nivel mensual, en la salida de la cuenca con las simulaciones de éste. Los estimadores de la bondad del ajuste, aplicados en la serie de validación y calibración, indican que el caudal simulado adapta su comportamiento a los datos observados. El periodo con caudales medidos (01/01/1996-31/12/2012) es dividido para destinar seis años a validación (01/01/1996-31/12/2001) y 11 años a calibración (01/01/2002-31/12/2012). De este modo, se escoge para calibrar SWAT la serie de mediciones más actual, que presenta menores lagunas y errores de medición. En el caudal simulado en ambas series se destina el primer año como periodo de calentamiento (NYSKIP). Por tanto, este año es eliminado en cada simulación evitando así parcialmente el sesgo inicial del modelo en la simulación (^{Arnold et al., 2012}, Figura 4).

Figura 4 Ajuste de SWAT en validación y calibración con NYSKIP=1 y NSE como función objetivo. 

El ajuste con NSE como función objetivo muestra, en calibración, que de las 2000 simulaciones lanzadas (distribución de probabilidad o 95PPU: ^{Abbaspour, 2015}) la mejor simulación obtenida con

SUFI-2 se adapta satisfactoriamente al caudal observado (Figura 4). Los resultados de la bondad del ajuste mejoran sensiblemente en la serie de calibración con respecto a la de validación, siendo clasificados en ambos casos como buenos y muy buenos según ^{Moriasi et al. (2007)}. Asimismo, en ambas series SWAT infraestima el caudal observado (PBIAS positivo). Sin embargo, el modelo ajusta correctamente el caudal base simulado al observado (ausencia de caudal base nulo), es decir, consigue reproducir los aportes de origen subterráneo recibidos en el embalse del Taibilla durante el periodo de estudio. Además, los caudales punta son simulados de un modo aceptable en dicho periodo pese a la subestimación de los meses 02 y 03 del año 2010. Esta infraestimación indica probablemente, por parte del modelo, una respuesta lenta ante precipitaciones intensas (posible conversión deficitaria lluvia-escorrentía).

La eficiencia de Kling-Gupta (KGE) fue escogida como función objetivo ya que a menudo proporciona un ajuste robusto frente a observaciones con numerosos caudales punta (Figura 5).

Figura 5 Ajuste de SWAT en validación y calibración con NYSKIP=1 y KGE como función objetivo. 

Debido a la naturaleza de KGE como función objetivo, el caudal simulado (Figura 5) presenta mayor número de picos (caudales punta) que el simulado con NSE (Figura 4). Esta fluctuación del simulado con KGE provoca que, en la serie de validación, SWAT infraestime el caudal observado de un modo perceptible (PBIAS=6.18), mientras que en calibración lo sobreestima tímidamente (PBIAS=-0.55). Por ello, la representación del caudal base con KGE empeora respecto a la simulación realizada con los parámetros hidrológicos empleando NSE como función objetivo. No obstante, la mencionada fluctuación del caudal simulado favorece que en los meses 02 y 03 del 2010 el simulado se aproxime al observado. Como consecuencia de lo enunciado, en validación y calibración, KGE muestra una bondad del ajuste inferior a la obtenida con NSE (Tabla 3). Estadísticos como NSE, KGE y R ² se alejan de rangos satisfactorios según la clasificación de ^{Moriasi et al. (2007)} y ^{Gupta et al. (2009)}.

De este modo, la modelación hidrológica en la cuenca, para el periodo presente y futuro, se realiza a partir de los resultados proporcionados con NSE como función objetivo. Dichos resultados se trasladan a los parámetros hidrológicos a través de la columna valor fijado en la Tabla 4. Se corrigen así los valores originales de los 14 parámetros seleccionados para representar el comportamiento hidrológico de la cuenca alta del río Taibilla (modelización hidrológica inversa). Estos parámetros, junto con sus rangos de valores iniciales en la Iteración uno, fueron escogidos a partir de experiencias de modelización con objetivos similares a los enunciados en este estudio (^{Abbaspour, 2015}; ^{Abbaspour et al., 2015}; ^{Gomariz-Castillo, 2015}; ^{Senent-Aparicio et al., 2017}).

Los parámetros de la Tabla 4 poseen diferencias notables en el análisis de sensibilidad (iteración 5). ALPHA_BF.gw, SOL_K (…).sol y CN2.mgt se muestran significativamente sensibles a KGE como función objetivo (P-valor≤0.05). Por ello, un cambio en dichos parámetros afecta en gran medida a los resultados del modelo. El análisis con NSE presenta un número mayor de parámetros con sensibilidad significativa que en el caso anterior (P-valor≤0.05), favoreciendo así la elección de NSE como función objetivo. ALPHA_BF.gw, CN2.mgt, SOL_BD(…).sol y SOL_AWC(…).sol son los que poseen mayor sensibilidad. Dado que los parámetros ALPHA_BF.gw y CN2.mgt también mostraron sensibilidad significativa con KGE como función objetivo, se considera que ambos resultan esenciales en el comportamiento y dinámica hidrológica de la cuenca alta del Taibilla. Análogamente, los parámetros de las capas de suelo (base de datos .sol) presentan un comportamiento relevante con NSE y KGE como funciones objetivo.

Análisis de las proyecciones de cambio climático

Los GCMs seleccionados, en la zona de estudio con la regionalización estadística de análogos, se denominan: bcc-csm1-1-m, MPI.ESM.MR y bcc.csm1.1. En la Figura 6 se muestran, para las variables precipitación y temperatura media, los promedios anuales de las simulaciones resultantes de estos GCMs a escala diaria en la cuenca alta del Taibilla. Dichas simulaciones se obtuvieron a nivel de unidad de respuesta hidrológica (HRU) por lo que fueron agregadas espacialmente a escala de cuenca.

Figura 6 Evolución de las proyecciones de cambio climático en la cuenca (periodo 2010-2050). 

Durante las próximas décadas, en los escenarios de emisión RCP4.5 y RCP8.5, el promedio de la precipitación para los GCMs implementados presenta una variabilidad interanual notable (Figura 6). El análisis mediante el contraste no paramétrico de Mann-Kendall (^{Kendall, 1976}; ^{Sneyers, 1992}) no detecta tendencias significativas. El escenario RCP8.5 presenta una tendencia moderadamente decreciente con una reducción de la precipitación media de 8.4 mm/década, mientras que el escenario RCP4.5 muestra una tendencia prácticamente constante (+0.3 mm/década). El promedio de la temperatura en cambio registra un claro incremento en el periodo de estudio (2010-2050). Los dos escenarios considerados prevén un aumento de la temperatura estadísticamente significativo (P-valor<0.001) en dicho periodo (Figura 6). El incremento previsto es de 0.31 y 0.25 °C/década para los escenarios RCP8.5 y RCP4.5 respectivamente. Estas modificaciones resultan coherentes con los cambios proyectados, a lo largo del siglo XXI, en ^{CEDEX (2012)} y ^{Senent-Aparicio et al. (2017)} en la cabecera de la Demarcación Hidrográfica del río Segura (DHS). Del mismo modo, las tendencias identificadas (Figura 6) concuerdan con las variaciones, para precipitación y temperatura, estimadas a escala global en ^{Allen et al. (2000)} e ^{IPCC (2013)} debido al aumento de las emisiones y concentración atmosférica de gases de efecto invernadero (GEI).

Así mismo, la evolución de los diferentes modelos, en las variables y escenarios de emisión analizados (Figura 6), indica que el que el GCM bcc.csm1.1 muestra mayor variabilidad interanual (mayor número de picos) que los restantes con respecto a la serie promedio. A pesar de ello, las últimas décadas del periodo de estudio, especialmente la 2040-2050, presentan una convergencia y estabilización de las proyecciones de los modelos utilizados, aspecto ya señalado en ^{IPCC (2013)}. No obstante, este informe señala que las predicciones de temperaturas y precipitaciones futuras poseen una gran incertidumbre asociada al GCM utilizado. La variable precipitación muestra, además, elevada incertidumbre en función de la técnica de regionalización implementada (^{IPCC, 2013}).

En virtud de las proyecciones de cambio climático analizadas, para las variables precipitación y temperatura, se estima que el clima actual de la cuenca alta del Taibilla, clasificado en ^{Boix-Fayos et al. (2007)} como mediterráneo subhúmedo en transición al semiárido, derivará hacia una climatología de tipo semiárida-árida. Dichos cambios se deben principalmente a los incrementos, con tendencia significativa (P-valor<0.001), identificados en la temperatura (Figura 6).

Respuesta hidrológica en los periodos de estudio

En este apartado, se analiza la respuesta hidrológica de la cuenca en el periodo actual (2000-2010) junto con el corto plazo (2021-2030) y el mediano plazo (2041-2050) como consecuencia de las proyecciones de cambio climático acopladas al modelo SWAT tras evaluar la bondad del mismo (principalmente mediante el fichero Output.std obtenido en cada simulación). Los promedios anuales de los elementos del balance hídrico mostrados (Tabla 5) proceden de la agregación temporal y espacial realizada, a escala de cuenca, con el software SWAT-CHECK (^{White et al., 2012}) durante los periodos mencionados. En la Tabla 5 se incluyen, para el corto y el mediano plazo, los resultados procedentes del promedio de los tres GCMs considerados en el estudio. Los valores medios de precipitación, escorrentía superficial y recarga de acuíferos se relacionan en mm/año, mientras que la temperatura media en °C. Así mismo, se indican las modificaciones que dichas variables presentan en los periodos y escenarios futuros con respecto a sus respectivos valores en el periodo presente (2000-2010).

Los resultados obtenidos (Tabla 5) aportan valores y tendencias verosímiles en la cuenca alta del Taibilla dada la climatología actual y futura de la misma. Así mismo, los cambios registrados en las variables del balance hídrico (Tabla 5) concuerdan con las modificaciones mostradas en ^{CEDEX (2012)}, ^{CHS (2015a)} y ^{Senent-Aparicio et al. (2017)} para la cabecera de la DHS durante el periodo 2000-2050. Ello indica, junto con los resultados de la calibración y validación del modelo, una correcta elección de la función objetivo, la eficiencia de NSE, así como de los 14 parámetros hidrológicos seleccionados para representar el comportamiento hídrico de la cuenca.

La precipitación media anual en la cuenca alcanza, actualmente, los 585 mm según los registros facilitados por AEMET mientras que SWAT estima en el periodo 2000-2010 596.46 mm/año (Tabla 5) debido a que el modelo completa los valores ausentes existentes en las series de AEMET. Los periodos 2021-2030 y 2041-2050, para el escenario RCP4.5, muestran reducciones de la pluviometría media de en torno al 3.8% con respecto a los 596.46 mm/año del 2000-2010. En cambio, el descenso del escenario RCP8.5, en ambos periodos, se sitúa en el 12.7%. A pesar de ello, en dichos escenarios, las modificaciones estimadas en la precipitación para el promedio del periodo 2010-2050 mostraron tendencias no significativas.

A diferencia de la pluviometría, la temperatura media en la cuenca registra para las próximas décadas incrementos significativos (P-valor<0.001) en los escenarios de emisión considerados. En el RCP4.5, con respecto al periodo actual, la temperatura aumentará 1.47 °C en el corto y el mediano plazo, mientras que en el RCP8.5 ascenderá 1.72 °C (Tabla 5). Dichos aumentos se aproximan al umbral de los 2 °C fijados como punto de no retorno a partir del cual la influencia del cambio climático sobre los recursos hídricos y ecosistemas asociados tendrá efectos irreversibles (^{Green et al., 2011}; ^{IPCC, 2013}; ^{Klove et al., 2014}). Así mismo, los citados cambios afectarán a las tasas de evapotranspiración real (ETR) y potencial (ETP) estimadas en el periodo presente con SWAT en 402.33 mm/año y 947.97 mm/año respectivamente. En coherencia con los resultados de ^{CEDEX (2012)} para el sureste mediterráneo, durante las próximas décadas, se prevén en la cuenca reducciones de la ETR, debido a la menor disponibilidad de agua de lluvia, e incrementos en la ETP por el aumento previsto de la temperatura media.

Las modificaciones estimadas para precipitación y temperatura provocarán una disminución de los recursos hídricos superficiales y subterráneos en la cuenca, cuya disponibilidad neta se verá además condicionada por las elevadas tasas de ETR propias del sureste mediterráneo. El 64% de la lluvia caída genera escorrentía superficial (número de curva promedio de 64.06) debido a la humedad, tipos y usos del suelo (zona próxima al régimen natural). En el periodo 2000-2010, la escorrentía superficial se sitúa en los 40.04 mm/año (Tabla 5). En el escenario RCP4.5 para el corto y el mediano plazo desciende un 10.2%, mientras que en el RCP8.5 disminuye un 18.6%.

La recarga de acuíferos presenta en el periodo actual un valor promedio de 103.01 mm/año (Tabla 5). En el escenario RCP4.5, durante las series 2021-2030 y 2041-2050, se reduce un 4.4% y en el RCP8.5 un 16.3%. Estos cambios provocarán reducciones en las tasas de retorno de acuíferos y por ende del caudal base o ecológico de los cursos superficiales afectando negativamente a los ecosistemas asociados (^{Klove et al., 2014}). No obstante, los descensos estimados en la recarga de acuíferos resultan menores que los obtenidos para la escorrentía superficial durante las próximas décadas. Además, los elevados tiempos de tránsito y residencia del flujo hídrico en el subsuelo, variables desde días a miles de años, impiden identificar de un modo directo los efectos del cambio climático sobre las aguas subterráneas (^{Green et al., 2011}). De hecho, en la cuenca, el escenario RCP4.5 prevé un incremento del 1.1% en las tasas de recarga durante el periodo 2041-2050. Del mismo modo, determinados estudios realizados en la DHS proyectan aumentos puntuales de la recarga de acuíferos durante los próximos decenios (^{Pulido-Velázquez et al., 2015}). Ello indica que actualmente existe una gran incertidumbre asociada a los escenarios de emisión, GCMs y técnicas de regionalización estadística empleadas en las proyecciones de cambio climático.