En este artículo realizaremos un análisis de la "Riada de Santa María" de septiembre de 2019, la más fuerte y catastrófica en 140 años en la Región de Murcia. Para ello haremos uso de un programa de modelización hidrológica, llamado HEC-HMS, que nos permitirá crear un modelo hidrológico que realizará simulaciones de los hidrogramas producidos (y por ende, la virulencia de la riada) por la cuenca del Campo de Cartagena, vertiente al Mar Menor, dónde se dieron los peores resultados del episodio.
Figura 1. Composición RGB extraída del satélite SENTINEL 2 mediante QGIS donde se muestra el impacto que tuvo el episodio. Fuente: Elaboración propia.
Para la creación del modelo hidrológico utilizaremos la Rambla del Albujón, red de drenaje más importante del Campo de Cartagena. Su estudio nos permitirá entender y comprender con mayor claridad los eventos del pasado septiembre de 2019. Además, la creación del modelo anteriormente mencionado podrá utilizarse en el futuro para estimar avenidas, así como para comprender mediante simulaciones las riadas pasadas. Sin más, comenzamos:
HEC-HMS
HEC-HMS es un programa diseñado para calcular el hidrograma producido por una cuenca en un determinado episodio pluviométrico. Esto lo consigue a partir de unos datos de entrada relacionados con las características de la cuenca y de las precipitaciones. HEC-HMS simplifica el sistema a partir de la definición de cuatro fases (Figura 2).
En la primera fase (A) se obtiene la precipitación neta (Pe), parte de la precipitación caída que se convierte en escorrentía superficial (el resto, “pérdidas", pasa a infiltrarse al acuífero o se almacena en el suelo y cobertura vegetal). En la fase (B) la Pe se transforma en escorrentía directa (B) y en (C) se obtiene el flujo base procedente del acuífero. En la fase (D) se aproxima el transito del hidrograma por los cauces. Para la obtención del modelo mediante las fases anteriormente mencionadas, HEC-HMS cuenta con una serie de componentes en los cuales han sido introducidos los datos.
Figura 2. Fases en las que HEC-HMS divide el sistema hidrológico. Fuente: Sánchez San Román (2017).
Proyecto:
El primer paso para modelizar un evento en HEC-HMS es la creación de un nuevo proyecto. En este estudio se han establecido como parámetros básicos los siguientes:
Para las abstracciones no se considera el componente de cobertura vegetal (Canopy), que simula la interceptación de lluvia y la evaporación entre eventos, ni tampoco el almacenamiento en depresiones en la superficie del terreno (Surface).
Para las pérdidas y la transformación de la precipitación se ha usado el método de número de curva del Soil Conservation Service (SCS).
No se considera flujo base ya que el caudal de la Rambla del Albujón es efímero.
Como se ha creado un modelo semidistribuido, existe traslación del hidrograma por las subcuencas, así como pérdidas y ganancias durante el tramo del cauce; para ello se ha usado el Método de Muskingum para aproximar el tránsito por los cauces.
Las entradas de precipitación han sido determinadas por los hietogramas obtenidos de las 11 estaciones pluviométricas y su interpolación a partir del IDP.
Componente de cuenca (subbasin):
La Rambla del Albujón es una cuenca de gran extensión perteneciente al Campo de Cartagena. Para que la modelación de esta sea correcta, es imprescindible dividirla en subcuencas. Se han obtenido 24, determinadas por las distintas características requeridas por HEC-HMS:
Área (Km2).
Número de Curva (CN).
Porcentaje de impermeabilidad (%).
Tiempo de retardo (Tlag) (min).
La obtención del área se realizó mediante el procesamiento de la capa vectorial de subcuencas en un sistema de información geográfica QGIS. Al seleccionar como método de pérdidas y transformación el Método del SCS, el procesamiento del resto de parámetros se ha realizado de la siguiente manera:
A) Método de pérdidas
El método de pérdidas del SCS permite obtener qué parte de la precipitación se convierte en escorrentía directa (Pe) y qué parte en abstracciones, mediante la obtención de únicamente dos parámetros: el número de curva (CN) y el porcentaje de impermeabilidad. CN es un valor adimensional que sirve para obtener el “umbral de escorrentía” (P0) (valor a partir del cual se produce escorrentía). Varía entre 0 (condiciones de permeabilidad muy altas) y 100 (zonas totalmente impermeables). El método se basa en las hipótesis de que la precipitación comienza a producir escorrentía directa (Pe) cuando la precipitación total caída hasta ese momento supera un umbral o abstracción inicial (P0). Si en un momento del transcurso de la precipitación la capacidad del suelo está al 30% de su capacidad máxima, hasta ese mismo instante habrá generado escorrentía directa el 30% de la precipitación (Sánchez San Román, 2017). Se ha comprobado empíricamente que P0 es aproximadamente el 20% de la abstracción máxima, por lo que puede obtenerse conociendo CN a partir de la siguiente ecuación:
𝐶𝑁 = 25400 254 + 𝑃0/0,2 (1)
Por medio de esta ecuación, y la obtención de una capa de valores P0 mediante algebra de mapas, resultante de la capa de los grupos hidrológicos del suelo, tipo de uso del suelo y la pendiente, se ha obtenido una capa de valores en QGIS que representa el CN medio por subcuenca. Para obtener una mayor aproximación al CN real de la subcuenca ha sido necesario corregir dicho valor según las condiciones precedentes del suelo, que en nuestro caso fueron secas, por lo que se utilizó la siguiente ecuación:
𝐶𝑁 = 4,2 ∙ 𝐶𝑁 𝑟𝑒𝑎𝑙 10 − 0,058 ∙ 𝐶𝑁 (2)
Al igual que el NC, el porcentaje de impermeabilidad se ha obtenido tras la creación de una capa de zonas impermeables de las subcuencas de la Rambla del Albujón obtenida de CLC2018.
B) Método de transformación
El método del SCS para la transformación de la precipitación permite averiguar qué parte de la precipitación neta se transforma en escorrentía directa. Para ello es necesario obtener el tiempo de retardo (Tlag), definido como el tiempo de desfase entre el centroide de la precipitación neta y el caudal punta del hidrograma. Se ha demostrado empíricamente que en cuencas de pequeño tamaño y características relativamente homogéneas el Tlag es aproximadamente el 60% del tiempo de concentración (Tc). Este último se define como el tiempo mínimo necesario para que todos los puntos de la cuenca estén aportando escorrentía simultáneamente al caudal recibido en la salida, y se ha comprobado experimentalmente que es independiente de la magnitud de la precipitación y depende únicamente de las características de la cuenca. (Témez, 1987). Para calcularlo, Témez (1987) planteó la siguiente fórmula:
𝐿 0,76 𝑇𝑐 = 0,3 ∙ (𝑆 0,25) (3)
Donde:
Tc se mide en horas
L = longitud del cauce principal hasta la divisoria de aguas (Km)
S = pendiente del cauce principal en tantos por uno.
Los datos necesarios como la longitud del cauce (L) y la pendiente media (S) han sido procesados por medio de QGIS. Finalmente, mediante la fórmula de Témez y la mencionada relación entre el tiempo de concentración y el tiempo de retardo (𝑇𝑙𝑎𝑔 = 0,6 ∙ 𝑇𝑐) se ha obtenido para su introducción en HEC-HMS el Tc para cada subcuenca (ver ANEXO I. Tablas).
C) Método de tránsito o translación
HEC-HMS aproxima el tránsito del hidrograma a partir del elemento cauce (Reach). Dicho elemento permite determinar el tiempo y la magnitud del caudal en un punto concreto del curso de agua. Para el tránsito del hidrograma se ha usado el Método de Muskingum. Este método modela el almacenamiento volumétrico (S) del flujo en un tramo de un cauce por medio de la unión de dos sumandos determinados por el almacenamiento en prisma, que sería proporcional al caudal de salida (O), y el almacenamiento en cuña, cuyo volumen final vendría dado por la diferencia del caudal de entrada y de salida (I – O). Esta relación viene determinada por dos parámetros; el coeficiente de proporcionalidad (K) y el factor de ponderación (X). K es el coeficiente de proporcionalidad medido en horas y representa el tiempo de tránsito de la onda de flujo a través del tramo del cauce. Se puede estimar de la misma manera que Tlag, asumiendo que será el 60% de Tc, de la subcuenca por la que pasa el hidrograma. X es un factor de ponderación que depende de la forma del almacenamiento en cuña y se encuentra entre los valores de 0 y 0,5. Normalmente se usan valores entre 0,2 (para tramos con poca pendiente) y 0,3 (para tramos con pendiente elevada). Para su incorporación a los elementos “reach” de HEC-HMS se ha optado por dar el valor 0,2, dado que será recalculado en la etapa de calibración.
Componente modelo meteorológico
Este componente es el encargado de calcular para cada subcuenca la precipitación asociada a partir de los datos de entrada. Como ya se ha descrito, parte de los hietogramas introducidos se han obtenido mediante el método IDP. IDP es un método de interpolación mecánica a partir de medias ponderadas basado en el principio de la autocorrelación espacial (Tobler, 1970). El valor de dicha ponderación se estima como la inversa de la distancia que los separa del punto a estimar. En el caso de HEC-HMS se generan series obtenidas mediante IDP para los centroides de cada subcuenca.
Componente especificaciones de control
En este componente hemos establecido la fecha y hora de inicio y final de la simulación además de la duración de los intervalos de tiempo, llamada resolución temporal: 08/09/2019 a las 0:00 y la fecha fin el 17/09/2019 a las 23:00, con una resolución temporal de 10 minutos.
Calibración del modelo y bondad de ajuste de la simulación
Esta fase del trabajo se basa en la ejecución del modelo cambiando los parámetros estimados en las fases anteriores a partir de un método de búsqueda de aquellos que minimizan el error de la simulación respecto a los datos observados de aforo. Para ello utiliza una función objetivo a minimizar (F0); ésta suele ser un estadístico de error. En este trabajo se han utilizado las siguientes (definidas en detalle en Feldman, 2000):
Porcentaje de error en el pico de descarga (PEP), que se calcula como el ratio entre la diferencia del pico simulado y el observado dividido por el pico observado. Es una F0 utilizada para priorizar los valores de pico sobre la forma del hidrograma, de gran utilidad en estudios de riesgo de inundación.
Porcentaje de error en volumen de descarga (PEV), similar al anterior, pero utilizando el volumen simulado y observado del episodio. Suele utilizarse para priorizar el volumen total de agua descargada, de gran utilidad en estudio de recursos hídricos.
RMSE ponderado por el pico (PwRMSE), función objetivo de compromiso entre elvolumen y los valores máximos o de pico, al minimizar la distancia media entre el caudal observado y simulado, pero dando un mayor peso a los flujos mayores al promedio del hidrograma.
Como método de búsqueda se ha utilizado el de Nelder y Mead (1965), método numérico basado en simplex para minimizar F0 en un espacio multidimensional a partir de mínimos locales.
Para evaluar la bondad de ajuste se han utilizado los siguientes estadísticos (desarrollado en detalle en Gomariz-Castillo, 2016)
Coeficiente de determinación r2, que determina la proporción de variación de los resultados que puede explicarse por el modelo. Varía de 0 (sin concordancia) a 1.
Índice de eficiencia de Nash y Sutcliffe (1970): Similar a r2 , varía entre 0 (sin concordancia) y 1 (concordancia perfecta); es más realista que el anterior al tener en cuenta desviaciones constantes.
Error cuadrático medio (RMSE): Estima el error cometido por el modelo en m3 /s. Varía entre 0 y +infinito. 23
Porcentaje de sesgo (PBIAS): Su valor óptimo es 0 y valores cercanos a éste indican precisión de ajuste. Valores positivos indican sesgo debido a una subestimación en la serie estimada, mientras que valores negativos indican sobreestimación.
RESULTADOS
En este apartado se presentan los resultados obtenidos tras la calibración del modelo implementado. Aunque inicialmente se partió de 6 estaciones de aforo del SAIH para su calibración, debido a diversos problemas detectados se han utilizado solo los aforos incluidos en las “estaciones intermedias”. En este trabajo se han identificado tres grandes grupos de estaciones de aforo en función de su localización:
Estaciones de cabecera 06A06A01 y 06A05A01. En ellas se ha observado cómo el modelo sobreestima casi 10 veces el hidrograma de caudal; sin embargo, el coeficiente de determinación r2 es elevado aún sin calibrar.
Estaciones intermedias 06A04A01 y 06A03A01. Son las finalmente utilizadas en el trabajo para obtener el modelo final.
Estaciones próximas a la desembocadura 06A01A01 y 06A02A01. El problema detectado en ellas es la falta de datos: disponen de datos de aforo en el inicio y final del hidrograma de escorrentía, pero en la parte central, coincidente con el pico del episodio, no dispone de información. Este problema puede deberse al desbordamiento de la rambla. En los siguientes subapartados se describen estos grupos de estaciones de aforo y un apartado
ESTACIONES DE CABECERA:
La superficie aportante a los puntos de aforo 06A05A01 y 06A06A01 es de 24.102,43 ha y 3.111,65 ha, representando el 34,79 % y el 4,49 % de toda la cuenca respectivamente. El problema en este caso es la elevada sobreestimación del hidrograma de caudal sin calibrar. Sin embargo, r2 es muy elevado y la sobreestimación observada no se propaga aguas abajo (en los puntos de aforo 06A04A01 y 06A03A01), con valores de sesgo y diferencias de caudal pico muy inferiores. Este comportamiento puede deberse a una serie de errores cuyas causas pueden ser:
Especificación errónea de los parámetros iniciales, principalmente NC, que controla la proporción de agua precipitada que se infiltra.
Errores en la toma de datos de los pluviómetros, que dan como resultado volúmenes mayores a los reales.
Datos observados erróneos.
Especificación errónea de la definición de las subcuencas vertientes.
En lo referente a los datos erróneos de los aforos, no podemos detectar si éstos son erróneos y hemos decidido utilizarlos como información ya contrastada por el órgano competente. La cuarta causa creemos que es la que afecta en nuestro estudio. Es un error bastante común en modelización hidrológica de zonas donde las vertientes de las subcuencas no pueden definirse bien, debido a la baja pendiente de la zona, o elementos geográficos no recogidos por el MDT25. Para comprobar este problema es necesario volver a generar los cauces y subcuencas a partir del MDT a una mayor resolución (por ejemplo, el MDT5). En este estudio utilizamos las subcuencas y cauces oficiales del CEDEX, realizadas a partir del MDT25, y es posible que en zonas como la del Campo del Mar Menor, puedan ocurrir este tipo de errores.
A continuación, vamos a analizar los estadísticos de bondad de ajuste de las estaciones de cabecera, dónde son especialmente interesantes los resultados mostrados en la Tabla 1: r2 en la simulación sin calibrar para ambas estaciones es bastante elevado, indicando un buen ajuste en la forma de los hidrogramas. Sin embargo, nse es muy inferior a 0, como resultado del elevado porcentaje de sesgo (PBIAS), que indica la elevada sobreestimación del modelo no calibrado; este resultado de nse indica que el rendimiento del modelo es peor para representar al hidrograma observado que el propio valor medio de los datos observados. Sin embargo, en los resultados obtenidos en las estaciones intermedias en la simulación sin calibrar (y por tanto partiendo de los parámetros no calibrados de los elementos aguas arriba de estas estaciones) no indican un comportamiento similar, lo que puede estar indicando algún problema de los indicados en los datos observados. Debemos resaltar que la calibración se ha centrado en el ajuste mediante la F0 Rems como solución de compromiso; aun así, se obtienen valores de nse y de r2 elevados, sobre todo teniendo en cuenta que estamos modelizando eventos con paso de tiempo horarios.
Tabla 1. Estadísticos de bondad de ajuste en las estaciones de aforo. Q. no calibrado (m3 /s) se refiere a la comparación de los datos aforados respecto a la simulación con los parámetros no calibrados; Q. calibrado (m3 /s) se refiere a la comparación de los datos aforados respecto a la simulación con los parámetros calibrados utilizando los hidrogramas de los aforos 06A05A01 y 06A06A01 (estaciones de cabecera).
Fuente: Elaboración propia
El ajuste se obtiene a costa de reducir notablemente el valor de CN de forma sistemática, de forma que se aumenta la infiltración (pérdidas) y ser reduce la precipitación efectiva. Se puede observar cómo se reducen significativamente los valores de NC para tratar de aumentar la infiltración y reducir la escorrentía superficial. En el modelo se puede observar que muy poca precipitación es transformada en escorrentía superficial, es decir, hay muy poca precipitación efectiva. Dicho problema mejoró tras la calibración, pero aun así los valores de infiltración resultantes fueron bastante altos. Esto se puede verificar obteniendo el volumen de precipitación que se convierte en escorrentía mediante la diferencia entre el exceso de volumen de precipitación (excess volume) y el volumen de la precipitación total (precipitation volume).
Figura 3. Hidrogramas de caudal (observado, simulado y simulado calibrado) generados en la cuenca vertiente hasta (a) el punto de aforo 06A05A01 y (b) el punto de aforo 06A06A01.
Observamos que en el aforo 06A05A01 la simulación calibrada se ajusta bien en el segundo pico, pero denota un pequeño desfase en la hora y el caudal del primero, siendo los datos de la simulación calibrada del primer pico de 8,26 m3 /s a las 08:40 y del caudal observado de 10,85 m3 /s a las 9:00. En cuanto al volumen de descarga podemos concluir que, debido al pequeño desfase horario, el modelo data de un volumen de 2,68 mm a la calibración frente a los 1,60 mm de los datos observados, y, por tanto, según los valores de los estadísticos de validación r2=0.84 y nse=0,703, es un buen ajuste. En el aforo 06A06A01 la simulación del modelo se ajusta notablemente a los datos aforados. Esto se observa en el coeficiente de determinación r2 que obtiene un valor de 0,755 y el nse con valor de 0.718, siendo valores bastantes altos que justifican la buena representatividad del ajuste a pesar del pequeño desfase horario. Por tanto, el pico de crecida tiene un valor de 1.7 m3 /s a las 7:20 en la simulación calibrada con un volumen de 2,33 mm y los datos de aforo muestran un pico de 1,6 m3 /s a las 6:00 con un volumen de 1,65 mm. Concluimos pues que el modelo se ajusta notablemente en los aforos 06A06A01 y 06A05A01, ya que realiza un buen ajuste de los picos observados y calibrados, finalidad conforme a los objetivos perseguidos en estudios de inundación.
ESTACIONES INTERMEDIAS
Al aforo 06A04A01 vierte el 49,68 % de toda la cuenca del Albujón (34.423,80 ha), y hasta el aforo 06A03A01 el 51, 55 % (35.715,24 ha). Curiosamente, el caudal pico máximo identificado en 06A04A01 es mayor al observado en 06A03A01, debido probablemente a la laminación del hidrograma a su paso por el cauce de la subcuenca 15, de 5.978 m y una pendiente muy baja de 0,027 m/m.
Como paso previo a la calibración final se probaron tres funciones objetivo (F0) siguiendo el esquema de calibración conjunto para los puntos de aforo 06A04A01 y 06A03A01. Tanto en las tablas como en las figuras se puede ver cómo el modelo tiende a fallar de forma generalizada en el ajuste del segundo pico del hidrograma. Ello es debido probablemente a dos aspectos fundamentales: la complejidad de la cuenca estudiada en lo referente a su hidrología y al método empleado para estimar las pérdidas por infiltración: el método de pérdidas del SCS. Éste se comporta correctamente a la hora de aproximar el pico máximo del evento al utilizar la F0 PEP; sin embargo, para conseguirlo sobreestima el segundo pico. Al utilizar la F0 PEV para aproximar el volumen generado subestima el pico máximo, con objeto de garantizar la mayor aproximación al volumen total, obviando para ello dicho pico.
Tabla 2. Ajuste de las calibraciones realizadas en función de la F0 seleccionada. Resumen en el punto de aforo A06A04A01: Pico observado=138,4 m3 /s; Volumen observado=6,51 mm.
Fuente: Elaboración propia
Figura 4. Hidrogramas de caudal, observado, simulado y simulados calibrado para las tres F0 utilizadas, generados en la cuenca vertiente hasta el punto de aforo 06A04A01.
Los estadísticos de validación para las funciones objetivo usadas en el aforo 06A04A01 nos muestran que solamente las F0 PwRMSE y PEV son válidas para realizar el ajuste ya que ambas tienen un coeficiente de determinación y un nse relativamente satisfactorio. El porcentaje de sesgo PBIAS arroja valores negativos para todas las F0, lo que determina la sobreestimación del modelo, aun así, los valores más cercanos a 0 se encuentran en las dos funciones mencionadas anteriormente, al igual que los valores del RMSE, los cuales se observan relativamente más altos en la F0 PEP. Por tanto, se ha determinado como función más representativa para este aforo la F0 PwRMSE.
Se puede observar cómo las funciones PwRMSE y PEV son las que más se ajustan al hidrograma. En este aforo el modelo ya calibrado nos muestra que el pico de descarga es de 68,36 m3 /s a las 9:30 con un volumen de 10,30 mm, mientras que los datos aforados muestran un pico de 138,41 m3 /s a las 07:00 con volumen de 6,51 mm. La generación de escorrentía directa viene dada en parte por la subcuenca 6, conectada directamente con el aforo, que proporciona un volumen de escorrentía de 35.12 mm y unas pérdidas por infiltración de 129.87 mm. Es decir, se genera un 21% de escorrentía del total de la precipitación caída.
Tabla 3. Ajuste de las calibraciones realizadas en función de la F0 seleccionada. Resumen en el punto de aforo A06A03A01: Pico observado=100,05 m3 /s; Volumen observado=7,01 mm.
Fuente: Elaboración propia
Figura 5. Hidrogramas de caudal, observado, simulado y simulados calibrado para las tres F0 utilizadas, generados en la cuenca vertiente hasta el punto de aforo 06A03A01
En el aforo 06A03A01 ocurre lo mismo que en el anterior; de las funciones objetivo usadas solo se ajustan al hidrograma observado las F0 PwRMSE y PEV, siendo la primera la que mejor lo representa. Esto se observa en los estadísticos de validación donde, como en el anterior aforo, el error cuadrático medio, el coeficiente de determinación y nse son mayores en la F0 PEP. Además, la función objetivo PwRMSE presenta un mayor r 2 (0,831) por lo que se usará esta como la más representativa.
Observando el hidrograma podemos comprobar que el modelo con la F0 PwRMSE ajusta satisfactoriamente el primer pico del hidrograma, determina una descarga del pico de 75.49 m3 /s a las 7:30 con un volumen de 11.29 mm frente a los 100.50 m3 /s a las 8:00 con volumen de 7.05 mm de los datos aforados. Las subcuencas 5 y 15, conectadas directamente con el aforo, generaron un 24% y un 21% respectivamente de escorrentía respecto a la precipitación caída en cada una de ellas, todo debido al gran volumen de precipitación que recibieron las subcuencas, llegando a tomar valores de 169.36 mm y 165.39 mm. Deducimos que la diferencia entre los picos es debida a problemas en la estimación de los parámetros necesarios para el método del SCS. También podemos interpretar claramente que es en estos aforos intermedios cuando se produjeron las grandes crecidas, causadas por las grandes precipitaciones, que provocaron los innumerables daños en la cuenca del Campo de Cartagena y que, como veremos más adelante, colapsaron los sistemas de medición de las estaciones de aforo más cercanas a la desembocadura. Como se puede observar, en ambos aforos la versión calibrada mediante la F0 PwRMSE resulta ser una solución de compromiso aceptable entre caudales pico, el volumen total y la forma del propio hidrograma, caracterizada nse. Se puede observar cómo se sobreestima el segundo pico, consecuencia de las propias especificaciones del modelo de pérdida SCS; este modelo tiene problemas con episodios con más de un pico, debido a que considera una vez que se satura el suelo en las pérdidas que no se produce infiltración en el resto de dicho episodio. Sobre el segundo pico en el hidrograma simulado hemos determinado que es un problema derivado del método de pérdidas utilizado. En este estudio hemos tratado de simular un proceso tan complejo como el del episodio de septiembre de 2019 con el método de pérdidas del SCS. Su ventaja principal es el reducido número de parámetros (NC y Tlag); sin embargo es una aproximación demasiado simple que, aunque reproduce bien picos del hidrograma, tiene problemas en episodios con más de un pico: una vez estimadas las abstracciones iniciales en el primer pico e iniciada la escorrentía, el modelo asume en toda la simulación que el suelo se ha saturado y en el segundo pico no se produce abstracción inicial, por lo que para un hietograma con un menor volumen de precipitación se producen picos en proporción mayores.
ESTACIONES PRÓXIMAS A LA DESEMBOCADURA
Al aforo 06A02A01 vierte el 68,36 % de toda la cuenca del Albujón (47365,71 ha), y hasta el aforo 06A01A01 el 78, 05 % (54074,33 ha). Una vez evaluadas las diferentes calibraciones, hemos decidido utilizar los parámetros calibrados para la F0 PwRMSE, como solución de compromiso. Es preciso puntualizar que, como hemos mencionado anteriormente, debido a la situación geográfica de los aforos, próximos a la desembocadura, estos han recibido una gran cantidad de m3 de agua y sedimentos que provocaron el desbordamiento del cauce y con ello el colapso de los sistemas de medición. Por ello, los datos de las estaciones de aforo de los siguientes hidrogramas (Figura 6 y Figura 7) están incompletos y no se ha podido calcular la diferencia del caudal pico y el volumen respecto al hidrograma observado. El resto de estadísticos, mostrados en las Tabla 4 y Tabla 5, se han obtenido a partir de los pares de puntos con observaciones.
Tabla 4. Ajuste entre las series simuladas sin calibración y con calibración. Resumen en el punto de aforo A06A02A01
Fuente: Elaboración propia
Los estadísticos de validación para el aforo 06A02A01 muestran que el caudal simulado tras la calibración tiene un valor de r 2 de 0.901, muy cercano a 1, por lo que podríamos decir que el ajuste es satisfactoriamente representativo. Sin embargo, nse es relativamente bajo (-24.964) y el PBIAS es muy bajo (-391.26). A pesar de ello, parece ser que en las zonas donde se da la existencia de datos, el modelo se ajusta notablemente y llega a ser representativo. En ocasiones, esta carencia de datos se puede suplir con una estimación, pero hemos preferido no sobrecargar el modelo para evitar errores.
Figura 6. Hidrogramas de caudal, observado, simulado sin calibrar y tras la calibración previa, generados en la cuenca vertiente hasta el punto de aforo 06A02A01.
Como podemos observar en el hidrograma, se aprecia una disminución del volumen de agua tras la calibración, siendo el valor para el caudal simulado no calibrado de 26.55 mm y el del calibrado de 19,84. Además, se alcanzan picos de caudales muy altos, pudiendo alcanzar hasta los 167,32 m3 /s en el caudal simulado, ya que el caudal pico antes de la calibración tenía un valor de 163.1 m3 /s. Por tanto, determinamos una reducción del volumen y un aumento de los caudales pico tras la calibración. Estos grandes caudales son el resultado de la precipitación caída en toda la cuenca.
Tabla 5. Ajuste entre las series simuladas sin calibración y con calibración. Resumen en el punto de aforo A06A01A01
Fuente: Elaboración propia
Tanto el nse (-24.964) como el porcentaje de sesgo PBIAS (-391.26) arrojan valores negativos, alejándose este último bastante del 0. El r 2 , sin embargo, nos ofrece un valor bastante alto (0.901) con lo que podríamos determinar una buena representatividad.
Figura 7. Hidrogramas de caudal, observado, simulado sin calibrar y tras la calibración previa, generados en la cuenca vertiente hasta el punto de aforo 06A01A01.
Al igual que en el anterior, en el aforo 06A02A01 disminuye el volumen del caudal simulado y aumenta el pico del caudal tras la calibración. El modelo data el volumen del caudal simulado no calibrado en 29.83 mm mientras que tras la calibración disminuyó a 23.95 mm. El primer pico llego a alcanzar los 184.5 m3 /s y tras la calibración aumentó a 209.1 m3 /s. Como se comentaba en la introducción del apartado de resultados, solucionar el problema de la falta información en los hidrogramas observados requiere técnicas de relleno de datos. Al disponer de dos estaciones intermedias hemos preferido no forzar la calibración a datos estimados. No obstante, resulta interesante analizar los resultados de la simulación respecto a la parte del hidrograma observado existente:
Las zonas con datos observados sirven para ver el ajuste del hidrograma
Siempre que el ajuste sea similar, es posible utilizar el hidrograma modelizado como estimación consistente del hidrograma de caudal.
HIDROGRAMA FINAL
En la Figura 8 podemos observar el hidrograma final de nuestro modelo, es el resultante del episodio de precipitaciones objeto de estudio, ya que se encuentra en la desembocadura de la Rambla del Albujón, vertiendo directamente al Mar Menor.
Figura 8. Hidrogramas de caudal, simulado sin calibrar y tras la calibración, generados en la desembocadura de la cuenca.
Si analizamos detenidamente los hidrogramas de los aforos conforme a la simulación calibrada del modelo, podemos determinar que el episodio de precipitaciones fué mucho más intenso en los dos primeros días (12 y 13) de lo que determina la simulación sin calibrar, y además, las zonas más cercanas a la costa fueron las más castigadas por las lluvias. Esto produjo que la capacidad de retención de agua del suelo en estas zonas se saturara, provocando la no generación de precipitación efectiva debido en gran parte a la textura francoarenosa-limosa que caracteriza a los xerosoles, tipo mayoritario de suelo en el Campo de Cartagena. Cuando el suelo comenzó a saturarse durante los días 12 y 13, unido a la intensidad de las precipitaciones, se llegó al nivel de máxima recarga, es decir, toda la cuenca empezó a contribuir en la generación del caudal mostrado en el hidrograma, por lo que se generó una gran cantidad de escorrentía que provocó el primer pico del hidrograma el día 13 de septiembre a las 12:20, alcanzándose un caudal de 424,6 m3 /s en la desembocadura y un volumen de 36.82 mm. Tras esto, la intensidad de las precipitaciones disminuyó y con ello los flujos de caudales, manteniéndose los aportes debido únicamente a las reservas de agua subterránea. Es lo que se conoce como la llamada “curva de agotamiento”. Todo parece indicar que, tras la breve disminución de la intensidad, a partir del día 14 las precipitaciones aumentaron, en menor medida que al inicio, pero provocando el segundo pico debido a la baja tasa de infiltración del suelo, donde se llegaron a alcanzar los 107.83 m3 /s el día 14 a las 14:15. Por otro lado, la simulación sin calibrar subestimaba el primer pico (383,1 m3 /s) y sobrestimaba el segundo (125,6 m3 /s), además, el modelo no calibrado databa en 39.37 mm el volumen del hidrograma en la desembocadura siendo este valor mayor que el volumen del modelo calibrado (36.82). Por otra parte, las pérdidas por infiltración durante el final del día 12 y principio del 13 llegaron a alcanzar los 55.94 mm a las 2:50 y 56.50 mm a las 7:20. Debido a las fuertes precipitaciones en todas las zonas de la cuenca, la cantidad de agua que llegó a los pueblos costeros provocó grandes inundaciones y daños.
En la Figura 9 podemos observar cómo prácticamente toda la zona litoral del noroeste está anegada por el agua, la cual arrastró grandes cantidades de sedimentos, que posteriormente, y unidos a una serie de factores externos, provocaron una mortandad masiva de seres vivos en la laguna. Existen autores que analizan la respuesta del Mar Menor a este tipo de situaciones, cada vez más recurrentes, y concluyen que en un futuro es posible que la laguna quede colmatada, en parte, debido a eventos de este calibre (Erena Manuel et al, 2020).
Figura 9. Composición RGB donde se muestra el NDWI calculado mediante algebra de mapas en QGIS para la visualización de las zonas inundadas tras la riada. Fuente: Elaboración propia.
CONCLUSIONES
Como se puede ver, el modelo tiende a fallar de forma generalizada en el ajuste del segundo pico del hidrograma. Ello es debido probablemente a la propia hidrología de la cuenca (caracterizada por su elevada fragmentación espacial en coberturas de uso) y al comportamiento del propio episodio de precipitación; la elevada complejidad hace que el método empleado para estimar las pérdidas por infiltración (método de pérdidas del SCS) no funcione todo lo bien que debiera. Se hace evidente su debilidad a la hora de simular episodios complejos con más de un pico, siendo conveniente probar con modelos más complejos en caso de necesitar estimar volúmenes más realistas o más de un pico de tormenta. En este sentido, en los episodios simples de un solo pico el uso del método del SCS es idóneo, mientras que en episodios más complejos puede ser conveniente probar otros métodos como el Soil Moisture Accounting (SMA); el problema de este método, utilizado para modelización orientada a recursos hídricos para estimar la demanda de agua (modelización continua en el tiempo) es que cada subcuenca requiere de 14 parámetros, frente a los dos del SCS.
Por otro lado, se ha comprobado que el uso de diferentes F0 cambia sustancialmente los resultados obtenidos en el ajuste final. Se aconseja por tanto el uso de una F0 determinada en función del objetivo perseguido en un estudio concreto:
PEP: obtiene muy buenos resultados para el ajuste del pico máximo (con valores en la diferencia entre el observado y simulado próximos a 0, siendo apropiada para estudios por riesgo de inundación.
PEV: obtiene buenos resultados en el volumen final simulado, debido a que trata de minimizar los errores en los residuos de todo el conjunto del hidrograma. Es recomendable para estudios relacionados con los recursos hídricos, como la disponibilidad de agua o el estudio de la infiltración al acuífero.
PwRMSE: Se observan en los resultados que es una F0 de compromiso entre ambos aspectos, al tratar de minimizar el error cuadrático medio de todo el hidrograma, pero ponderando (dando más peso) a las diferencias en los mayores caudales. Por ello obtiene un buen rendimiento en el ajuste del hidrograma al mismo tiempo que trata de compensar los caudales pico.
Además, este tipo de estudios pueden servir para inferir las partes del hidrograma no observados, como podemos ver para 06A01A01 y 06A02A01 o la extrapolación de la calibración a subcuencas próximas no aforadas, como es el resto de la cuenca hasta la desembocadura. En este sentido pueden desarrollarse futuras líneas de trabajo o investigación en las que el modelo puede ser eficientemente utilizado:
Validaciones del modelo en episodios anteriores para los que no existiesen datos de aforo. Desde la puesta en funcionamiento en 2016 de las estaciones de aforo del SAIH destacan multitud de eventos por su peligrosidad y susceptibilidad de ser utilizados para validar el modelo. Estos eventos son: diciembre de 2016, noviembre de 2018, septiembre de 2019, diciembre de 2019, enero de 2020.
Modelización en el resto de la Cuenca del Mar Menor. En este sentido, el trabajo presentado puede permitir extrapolar y regionalizar los parámetros calibrados al resto de ramblas (11), sin datos observados de aforo.