En las últimas décadas, uno de los desastres naturales que más pérdidas causa son las inundaciones (Ripple et al., 2020; Wallemacq et al., 2015), siendo en España el segundo desastre natural que más muertes provocó (215) entre 2000 y 2020 según datos de Protección Civil (Protección Civil, 2021). No sólo son pérdidas humanas lo que acarrea la presencia de estos eventos, sino también pérdidas naturales, económicas y paisajísticas. Sólo en la cuenca del Ebro, las inundaciones ocurridas entre 1980 y 2018 generaron alrededor de 1000 millones de euros en pérdidas, lo que representa el 77% de las pérdidas económicas registradas en Aragón por desastres naturales en ese período (Heraldo de Aragón, 2018b). Por esta razón, es necesario desarrollar y optimizar modelos predictivos que proporcionen información acerca de la evolución de estos eventos, ayudando a cuantificar los daños potenciales (Hu et al., 2023; Knijff et al., 2010; Thielen et al., 2009). Entre estos, los modelos de resolución numérica de las ecuaciones que gobiernan el flujo en lámina libre constituyen una herramienta útil, que ha sido desarrollada y perfeccionada en las últimas décadas (Kalyanapu et al., 2011; Özgen-Xian et al., 2021; Vacondio et al., 2016). A su vez, es imprescindible disponer de medidas reales y precisas de la evolución del flujo para poder calibrar dichos modelos.

Dada la complejidad del estudio de este tipo de flujos, se incorporan frecuentemente aproximaciones que facilitan la resolución de las ecuaciones que los describen. Las distintas hipótesis consideradas dan lugar a sistemas de ecuaciones que conforman los distintos modelos matemáticos. Aunque el fenómeno de una inundación asociada a la crecida en un río es, en esencia, un problema tridimensional (3D), es habitual su estudio a gran escala mediante el promedio de las ecuaciones en la vertical, obteniendo un modelo bidimensional (2D) (Hu et al., 2019; Noh et al., 2018; Vreugdenhil, 1994); e incluso, mediante el promedio en la sección transversal para reducir el fenómeno a una aproximación unidimensional (1D) (Horrit y Bates, 2002; Mark et al., 2004; Samal et al., 2025). Además, en regiones donde la velocidad del flujo es aproximadamente nula, como en un embalse, es frecuente el uso de modelos agregados o 0D, que permiten describir con relativa precisión y alta eficiencia la evolución de ciertas variables hidráulicas (Dazzi et al., 2024; García-Alén et al., 2023). Recientemente, diversos estudios han demostrado la utilidad de combinar modelos 0D con modelos 1D o 2D para representar sistemas hidrológicos complejos, como cuencas con embalses (Echeverribar et al., 2021; Vallés et al., 2025). Este acoplamiento permite aprovechar la eficiencia computacional de los modelos agregados en zonas de flujo cuasiestático, mientras se mantiene una descripción más detallada del flujo en los tramos fluviales. Todos estos modelos permiten representar un mismo fenómeno, diferenciándose en la precisión de sus resultados y en su ámbito de validez.

Por tanto, el estudio de la dinámica fluvial en grandes sistemas hidrográficos requiere la implementación de modelos computacionales de altas prestaciones que permitan una representación precisa y fiable del comportamiento del flujo bajo distintos episodios de avenida. En este contexto, el presente trabajo analiza el comportamiento del tramo del río Ebro comprendido entre Logroño y Mequinenza (Figura 1), incluyendo sus principales afluentes: Ega, Arga y Aragón. Esta región se caracteriza por sufrir frecuentes inundaciones debidas a las crecidas del río Ebro y sus afluentes, dada la presencia de extensas llanuras de inundación, salvo en la región aguas abajo del dominio, donde la presencia de la presa de Mequinenza y el acusado relieve generan un flujo encauzado y cercano al reposo. Dada la amplia extensión del dominio de estudio (más de 700 km de cauce), se ha optado por dividir el tramo en diferentes subdominios, cada uno modelizado con distintas aproximaciones numéricas para optimizar el tiempo de cálculo. La selección del tipo de modelo en cada tramo (0D, 1D, 2D o sus combinaciones) ha estado condicionada por la morfología del cauce, la extensión de la lámina de agua durante episodios de avenida y la necesidad de equilibrar precisión de los resultados con eficiencia computacional. Para evaluar el desempeño de los modelos, se ha realizado una comparación con mediciones procedentes de las estaciones de aforo gestionadas por la Confederación Hidrográfica del Ebro (CHE), las cuales registran niveles de agua a intervalos de 15 minutos (https://www.saihebro.com/homepage/estado-cuenca-ebro). Los resultados obtenidos evidencian la idoneidad de los modelos 2D en zonas con amplias llanuras de inundación, donde la representación del flujo lateral es crucial para capturar la evolución hidrodinámica del sistema. En contraste, en tramos encauzados y de menor anchura, el modelo 1D ha demostrado ser más eficiente y preciso. Asimismo, en el embalse de Mequinenza, la combinación de modelos 2D y 0D ha permitido reducir significativamente el tiempo de cálculo sin comprometer la precisión de los resultados. De este modo, se concluye que este trabajo contribuye al desarrollo de estrategias de modelización del flujo hidrodinámico en ríos, proporcionando una herramienta útil para la gestión de los recursos hídricos y la evaluación de eventos hidrológicos extremos. Los modelos computacionales utilizados en este estudio, resueltos mediante esquemas numéricos explícitos de primer orden basados en el método de volúmenes finitos, han sido desarrollados por el Grupo de Hidráulica Computacional de la Universidad de Zaragoza, y han sido validados extensamente en configuraciones bidimensionales (Echeverribar et al., 2019b; García-Navarro et al., 2019), unidimensionales (Echeverribar et al., 2021; Murillo y García-Navarro, 2010), y en regiones susceptibles de modelización 0D (Echeverribar et al., 2021; Vallés et al., 2025).

A continuación, se presentan las ecuaciones que conforman los modelos matemáticos utilizados para describir la evolución del flujo (Figura 2). Estas ecuaciones deben resolverse mediante métodos numéricos, ya que carecen de soluciones analíticas. Dichos métodos no se describen en el presente trabajo, pero sí en las siguientes referencias: Murillo y García-Navarro (2010) y García-Navarro et al. (2019).

El principal interés en analizar esta región radica en las frecuentes inundaciones que ha sufrido en las últimas décadas (Heraldo de Aragón, 2018b). Dado que el dominio abarca más de 700 kilómetros de cauce y aproximadamente 1600 km² de extensión, este se divide en 5 tramos (Figura 1b) con el objetivo de analizar de forma más precisa cada subdominio, así como aplicar diferentes aproximaciones numéricas que permitan optimizar la eficiencia computacional:

Todos los tramos mencionados se modelizan mediante un enfoque 2D y 1D, a excepción del tramo TM3, el cual se modeliza mediante modelos acoplados 2D-0D y 1D-0D debido a la presencia del embalse en su región final.

A continuación, se muestran los resultados obtenidos para el evento de inundación ocurrido en abril de 2018 en los tramos analizados, dado que se trata de una de las avenidas que más pérdidas ha provocado en la región de Aragón en las últimas décadas (Heraldo de Aragón, 2018a). Este análisis se realiza de manera individual en cada tramo, evaluando la precisión y la eficiencia de cada uno, con el objetivo de obtener un modelo global previamente validado. Para cada tramo, se muestra el hidrograma correspondiente a la condición de contorno de entrada durante dicho evento, así como la evolución temporal simulada y observada del caudal en las distintas estaciones de aforo ubicadas a lo largo del cauce (Figuras 10, 11, 12, 13, 14 y 15). Estos hidrogramas de entrada (Figuras 10a, 11a, 12a, 13a, 14a y 15a) corresponden a datos medidos y proporcionados por la CHE. Dado que cada tramo se analiza de forma individualizada, estos hidrogramas no son acumulados ni resultan de la simulación, sino que representan los caudales reales registrados en el río principal del tramo y en sus afluentes durante el evento. El resto de las condiciones iniciales y de contorno se recogen en la Tabla 2. Además, las láminas de inundación simuladas y observadas para los tramos principales TM1, TM2 y TM3 pueden consultarse de forma interactiva mediante la herramienta BeMaps, desarrollada por Geospatiumlab S.L. en colaboración con GEOT (Grupo de Estudios de Ordenación del Territorio de la Universidad de Zaragoza), en el enlace: https://bemaps.unizar.es/user-map/684be1f3c8c0ee487a5317f9/view.

A partir de los resultados mostrados en la sección anterior, se presenta un análisis de su precisión, utilizando parámetros cuantitativos, así como una discusión sobre la eficiencia computacional de los modelos. Para ello, se emplean distintos parámetros que permiten evaluar la precisión de los modelos tanto en el cauce como en la llanura de inundación.

Los valores obtenidos para las métricas de evaluación se presentan en la Tabla 3 para el modelo 1D y en la Tabla 4 para el modelo 2D. En la Tabla 3 no se incluyen las métricas relativas a la llanura de inundación, ya que este modelo, al estar promediado en la sección transversal, no permite representar el flujo lateral ni, por tanto, la inundación en esta región del dominio. Asimismo, no se incluyen dichas métricas para los tramos TM1.1, TM1.2 y TM1.3 en la Tabla 4 por falta de datos observacionales. Como se aprecia en ambas tablas, los errores MAE y PPD son, en general, menores en los tramos principales (TM1, TM2, y TM3) cuando se utiliza el modelo 2D. Por ejemplo, el parámetro MAE disminuye un 9%, un 26% y un 20% en los tramos TM1, TM2 y TM3, respectivamente, cuando se usa el modelo 2D en vez del modelo 1D. Esto se debe a que en estos tramos el flujo tiende a desbordarse hacia la llanura de inundación, donde el modelo 2D reproduce con mayor realismo los gradientes laterales de velocidad y la redistribución del caudal. Sin embargo, en los afluentes (TM1.1, TM1.2 y TM1.3), donde el flujo permanece mayoritariamente encauzado y la llanura de inundación no desempeña un papel relevante en la evolución de la avenida, el modelo 1D proporciona resultados más precisos, ya que la componente lateral del flujo tiene menor relevancia. A pesar de ello, el parámetro PTD muestra valores más elevados con el modelo 1D en comparación con el modelo 2D. Esta diferencia se debe a que el modelo 1D, al no representar la llanura de inundación, no captura el efecto de laminación que ocurre cuando el flujo se extiende sobre esta zona, mientras que el modelo 2D sí lo contempla, resultando en predicciones más ajustadas temporalmente (PTD más bajos). Por otro lado, las métricas de evaluación de la llanura de inundación reportadas en la Tabla 4 indican una elevada precisión del modelo 2D en la representación de esta región del dominio. En particular, se obtienen valores altos de HR (superiores a 0.9) y valores reducidos de FAR (inferiores a 0.2), dando lugar a valores de CSI cercanos o superiores a 0.8 en todos los tramos. Esta precisión es especialmente destacable en el tramo TM3, lo que demuestra que la utilización de un modelo 0D para representar el embalse no reduce la precisión de los resultados. Finalmente, con el objetivo de evaluar la eficiencia computacional de los modelos utilizados, se presenta la Tabla 5, donde se recogen los tiempos computacionales para cada tramo y modelo, así como el valor de la variable FTRT (Faster Than Real Time), definido como:

Los valores de FTRT (Tabla 5) evidencian una elevada eficiencia computacional en todos los casos, especialmente en los modelos 1D, donde las simulaciones se ejecutan entre tres y cinco órdenes de magnitud más rápido que el tiempo físico del evento. En el caso de los modelos 2D, los valores de FTRT oscilan entre 67 y 180, siendo superiores en los afluentes debido al menor número de celdas mojadas en el dominio respecto a los tramos principales.

Tras la realización y análisis de las simulaciones correspondientes al evento de avenida de 2018, se ha comprobado que los modelos desarrollados para cada uno de los tramos proporcionan resultados satisfactorios, tanto en el cauce como en las llanuras de inundación. En particular, los modelos construidos para los tramos principales TM1, TM2 y TM3 han demostrado una adecuada capacidad para reproducir tanto la altura como el momento de llegada de los picos de caudal en las estaciones de aforo de referencia. Las métricas obtenidas para la llanura de inundación corroboran también una alta precisión en la representación espacial del evento. En cuanto a los afluentes Ega, Arga y Aragón, la evaluación se ha limitado a la evolución del flujo en el cauce, dado que no se dispone de datos observacionales sobre la extensión de la inundación en estas regiones. Aun así, los resultados en estos tramos muestran una buena concordancia entre las simulaciones y las observaciones en las estaciones de aforo. Por otro lado, se ha mostrado que una elevada precisión no implica necesariamente un alto coste computacional, como se refleja en los valores elevados del parámetro FTRT, especialmente en los modelos 1D. En los tramos principales (TM1, TM2 y TM3), el uso de modelos 1D conlleva una pérdida notable de precisión, así como la imposibilidad de representar la dinámica en la llanura de inundación. Sin embargo, en los afluentes (TM1.1, TM1.2 y TM1.3), los resultados obtenidos con los modelos 1D son similares o incluso más precisos que los del modelo 2D, pese a su menor coste computacional. Esto se debe a que el flujo en estos afluentes está mayoritariamente encauzado, por lo que su comportamiento se ajusta bien a una modelización unidimensional, basada en una aproximación a lo largo de la sección transversal.

Algunas de las discrepancias observadas en la evolución del flujo en el cauce en tramos como Zaragoza-Mequinenza (Figura 15) podrían explicarse por alteraciones topográficas durante el evento, como la rotura de motas. Asimismo, errores o simplificaciones en la digitalización del cauce (véase sección Modelo digital de terreno), junto con la imposibilidad de recoger los cambios en la topografía, constituyen una limitación en los modelos, al igual que la falta de información sobre procesos hidrológicos como la infiltración o la precipitación, que pueden provocar variaciones significativas del volumen de agua dentro del dominio. Como línea futura de trabajo, se propone incorporar al modelo la capacidad de simular este tipo de fenómenos. Esto permitiría no solo mejorar la precisión de las simulaciones, sino también dotar al modelo de un mayor realismo y utilidad como herramienta predictiva frente a eventos extremos de inundación. A su vez, se prevé incluir procesos hidrológicos, como la infiltración, exfiltración o precipitación, que podrían contribuir a una mayor precisión en la predicción del flujo. Finalmente, se plantea desarrollar una simulación global que integre todos los tramos del sistema, empleando los hidrogramas calculados en cada tramo como condiciones de contorno para el siguiente. Este enfoque permitiría evaluar de forma más realista la propagación del flujo y validar el comportamiento conjunto del modelo.