Iber v3. Manual de referencia e interfaz de usuario de las nuevas implementaciones

Resumen

Iber es un modelo numérico bidimensional de simulación de flujo turbulento en lámina libre en régimen no-permanente, y de procesos medioambientales en hidráulica fluvial. Desde el lanzamiento de la primera versión, que incluía un motor de cálculo hidrodinámico completamente acoplado con procesos de transporte de sedimento y turbulencia, Iber ha ido evolucionando hasta convertirse en una herramienta de modelización de flujo en lámina libre de procesos ambientales de elevada complejidad. En este documento se presentan los desarrollos realizados para la versión 3, donde los avances son principalmente en cuatro líneas de investigación: un nuevo módulo de drenaje urbano, un avance significativo en las capacidades del módulo de hidrología, un nuevo módulo de erosión de suelos, y un nuevo módulo para el cálculo del transporte de sedimentos considerando material no uniforme (mezclas). Asimismo, todos los trabajos han ido acompañados por una tarea transversal de mejora de la interfaz, tanto para los módulos ya existentes, como la creación de nuevas ventanas y menús para los nuevos módulos con el objetivo de mejorar todo el flujo de trabajo.

Palabras clave: drenaje urbano, hidrología, erosión de suelos, granulometría no-uniforme, Iber

Abstract

Iber is a two-dimensional hydraulic model for the simulation of free surface flow in rivers and estuaries, and the simulation of environmental processes in fluvial hydraulics. Since the release of the first version of Iber, which included a hydrodynamic calculation engine fully coupled with sediment transport processes and turbulence, it has evolved to become a free surface flow modelling tool for highly complex environmental processes. This document presents the developments made for version 3, where the advances are applied mainly in four current research lines: a new urban drainage module, a significant advance in the capabilities of the hydrological process module, a new soil erosion module, and a new module for calculating sediment transport considering non-uniform material (mixtures). Likewise, all the work has been accompanied by a cross-cutting task of improving the interface, both for existing modules and the creation of new windows and menus for new modules aiming to improve the whole workflow.

Keywords: urban drainage, hydrology, soil erosion, non-uniform granulometry, Iber

1 Introducción

Iber es un modelo numérico de simulación de flujo turbulento en lámina libre en régimen no-permanente, y de procesos medioambientales en hidráulica fluvial [1,2]. El modelo fue desarrollado en primera instancia en coordinación por los grupos universitarios Flumen, ahora Instituto Flumen, de la Universitat Politècnica de Catalunya, el grupo GEAMA (Grupo de Enxeñería da Agua e do Medio Ambiente) de la Universidade da Coruña, el Centre Internacional de Mètodes Numèrics en Enginyeria (CIMNE), y el Laboratorio de Hidráulica del Centro de Estudios Hidrográficos del CEDEX.

El modelo Iber se presentó oficialmente (Iber 1.0) en junio de 2010. Entonces Iber contaba con tres módulos de cálculo principales: un módulo hidrodinámico, un módulo de turbulencia y un módulo de transporte de sedimentos (considerando sedimento de granulometría uniforme). Desde entonces, Iber ha incorporado nuevos módulos adicionales como el módulo de rotura de presas y balsas [3], el módulo de calidad de aguas [4,5], el módulo de hábitat hidráulico [6], y mejorando también las capacidades de los módulos ya existentes [7], así como la amigabilidad de la interfaz, y el tiempo de cálculo.

En este documento se presentan los desarrollos realizados para la versión 3, donde los avances son principalmente en cuatro líneas de contenido: un nuevo Módulo de Drenaje Urbano, un avance significativo en las capacidades del Módulo de Hidrología Distribuido, un nuevo Módulo de Erosión de Suelos, y un nuevo módulo para el cálculo del Transporte de Sedimentos considerando material no uniforme (mezclas). Asimismo, todos los trabajos han ido acompañados por una tarea transversal de mejora de la interfaz, tanto para los módulos ya existentes, como la creación de nuevas ventanas y menús para los nuevos módulos. Por ello, antes de la descripción de cada uno de los módulos, se describen los desarrollos de interfaz generales.

Los apartados de este documento, correspondientes a cada uno de los nuevos módulos, además del apartado de descripción general de la interfaz, se estructuran en dos partes diferenciadas: fundamentos teóricos e interfaz del usuario del propio módulo. En la presentación de las bases teóricas se describe, en primer lugar, los fenómenos objeto de los nuevos módulos, tanto su descripción física como las ecuaciones de gobierno que permiten describirlos, y, en segundo lugar, se relatan los esquemas numéricos para la resolución de las ecuaciones implementados en Iber. El apartado de interfaz es una descripción de los nuevos menús y ventanas desarrollados en Iber para cada uno de los nuevos módulos; así, que por un lado se presentan los desarrollos de mejora de la interfaz propiamente dicha, pero por el otro se dispone ya de una guía para el uso de las nuevas herramientas.

1.1 Breve descripción de los nuevos módulos

Iber ha pasado de ser un modelo hidrodinámico y de procesos morfodinámicos, a un modelo numérico que, estando basado en la resolución de las ecuaciones de aguas profundas, permite simular múltiples procesos ambientales de base hidrodinámica. La incorporación de nuevos módulos, así como la mejora de los ya existentes, se ha llevado a cabo de forma totalmente acoplada, siendo el nexo de unión entre todos ellos el módulo hidrodinámico. La Fig. 1 representa, de forma esquemática, las interrelaciones entre los distintos módulos que incorpora Iber.

Fig. 1. Interrelaciones entre los distintos módulos principales de cálculo de Iber v3.

1.1.1 Módulo de drenaje urbano

El nuevo Módulo de Drenaje Urbano se ha construido con la finalidad de ser complementario del modelo Iber ya existente, de cálculo bidimensional en superficie, para poder incorporar los procesos asociados al drenaje de agua en zona urbana a través de la red de alcantarillado.

El elemento principal de dicho módulo es el cálculo, en una dimensión, de la propagación del flujo a través de los conductos cerrados que conforman la red de alcantarillado, ya sea en lámina libre, o con flujo a presión. La consideración de flujo a presión se ha implementado a través de la metodología de la ranura de Preissmann.

Un proceso del drenaje urbano es el intercambio de agua entre la superficie y la red, ya sea solo de captación, a través de las rejas de alcantarillado, o de salida de la red a superficie en episodios de grandes lluvias. Ambos fenómenos han sido incorporados en Iber a través de dos metodologías distintas: considerando las rejas como orificios, o mediante la metodología de evaluación de la eficiencia hidráulica de las rejas desarrollada en los laboratorios del Instituto Flumen, en la UPC.

Para la simulación de las uniones de colectores, que también pueden funcionar en lámina libre o a presión, se ha optado por considerarlos como elementos puntuales de superficie, y calcularlos con un esquema bidimensional, semejante al utilizado por Iber en superficie.

En la red de alcantarillado se ha implementado una opción para poder considerar elementos especiales de la propia red (por ejemplo, depósitos de retención, uniones complejas de colectores, etc.), pudiéndose simular estos en un esquema bidimensional. Asimismo, este módulo se puede emplear para realizar la simulación únicamente en superficie, tanto en superficie como en la red, o tan solo en la red de drenaje.

1.1.2 Nuevos desarrollos en el módulo de hidrología

El Módulo de Hidrología Distribuido de Iber permite calcular la transformación lluvia-escorrentía en toda una cuenca, haciendo posible su utilización como modelo hidrológico distribuido en cuencas de pequeño y mediano tamaño. El módulo está basado en las ecuaciones hidrodinámicas de aguas poco profundas bidimensionales, las cuales se utilizan para resolver el movimiento de la escorrentía superficial sobre el terreno, al que se le han añadido una serie de nuevos procesos, como son: el aporte de agua de lluvia, la infiltración de agua en el terreno, el flujo subsuperficial, la evapotranspiración y la percolación.

El módulo hidrológico es completamente bidimensional y distribuido en el espacio, y trabaja sobre la misma malla de cálculo que el módulo hidrodinámico. La base del módulo hidrológico es el módulo hidrodinámico, al que se le han añadido los procesos de aporte de agua por precipitación y la infiltración del agua en el terreno. De esta forma, el exceso de precipitación se mueve por el terreno y por los cauces de acuerdo con las ecuaciones de aguas someras bidimensionales. El agua infiltrada pasa a formar parte del subsuelo, en donde se puede realizar un balance de masa de agua distribuido en el espacio a nivel de elemento de la malla de cálculo. En dicho balance de agua en el terreno se consideran los procesos de infiltración, evapotranspiración, percolación y flujo en el subsuelo. Los tres primeros implican un movimiento de agua en la columna vertical, mientras que el flujo en el subsuelo se considera bidimensional y paralelo al terreno. Todos estos procesos pueden activarse o desactivarse por el usuario.

1.1.3 Módulo de erosión en laderas por escorrentía

El nuevo Módulo de Erosión de Suelos de Iber permite determinar las zonas en las que se produciría erosión del suelo en una cuenca y cuantificar la misma, así como evaluar la evolución temporal y espacial de las concentraciones de sólidos en suspensión a lo largo de los cauces de la red fluvial.

El módulo puede trabajar con distintos niveles de complejidad de la estructura del suelo, mediante la activación o no de los diferentes procesos de erosión considerados. De esta manera, se pueden considerar una única capa de suelo erosionado (no cohesivo), o dos capas de suelo, una correspondiente a la matriz original del suelo (cohesiva) y otra correspondiente al suelo erosionado (no cohesivo). Las partículas de ambas capas pueden ponerse en movimiento por el impacto de la lluvia y por la fricción que genera el agua con el fondo. Una vez en movimiento, los procesos de transporte de sólidos considerados son por carga de fondo y en suspensión, si bien en la mayor parte de los casos el transporte en suspensión de partículas finas será predominante en las laderas de las cuencas. Asimismo, se consideran diferentes clases o tipos de sedimento, por lo que permite trabajar con granulometrías distribuidas si el usuario lo considera necesario.

1.1.4 Módulo de transporte de sedimento no uniforme (mezclas)

Desde el principio Iber incorporó un módulo de transporte de sedimentos por arrastre de fondo, y un módulo de transporte en suspensión, pero considerando la hipótesis, o simplificación, de material uniforme. En esta versión de Iber se ha implementado y validado un nuevo Módulo de Transporte de Sedimentos con granulemtría no uniforme, o granulometrías extendidas. Ello proporciona no solo mayor precisión en los resultados, sino también permite la simulación numérica de procesos que antes no se podían contemplar, como son: el acorazamiento del lecho, erosiones diferenciales por efecto de las distintas granulometrías, o la clasificación del material del fondo por zonas según la capacidad de arrastre.

Los desarrollos se basan en la teoría de la capa activa, y se ha realizado su implementación numérica mediante el método de los volúmenes finitos. De esta manera, el fondo de un cauce se puede dividir en distintos estratos, cada uno de ellos con una granulometría distinta representada por una serie de clases y su fracción. La ecuación de estimación del caudal sólido se hace para cada una de las clases, considerando el fenómeno de escondimiento para la interacción entre los distintos tamaños de grano.

2 Interfaz de Iber 3.0

2.1 Generalidades

2.1.1 Apariencia

La versión actual de Iber, al igual que versiones anteriores, se ha embebido dentro de la interfaz de pre- y post-proceso GiD [8-10]. Para ello se ha empleado la versión 15 de GiD, que incorpora múltiples ventajas y mejoras con respecto a versiones anteriores (más información en www.gidhome.com).

Fig. 2. Visión general de la interfaz de Iber v3.

La apariencia general del programa es, a priori, muy similar a la versión precedente (Fig. 2). El principal cambio en la interfaz, concretamente en la parte de pre-proceso, ha sido la simplificación de los menús y el agrupamiento de los diferentes módulos de cálculo en Plug-ins como se verá en el apartado 2.1.2.

En versiones anteriores de Iber los menús y submenús, así como las opciones generales de cada uno de los módulos de cálculo, se mostraban por defecto, aunque el/la usuario/a no los emplease durante la generación y cálculo de un modelo. Esto podía entorpecer el flujo de trabajo y, además, generar confusión a la hora de definir los parámetros.

Puesto que Iber v3 incluye nuevos módulos de cálculo, se ha creído conveniente agrupar todas y cada una de las opciones, menús y parámetros de cada módulo de cálculo. Esto ha redundado en una simplificación significativa de toda la interfaz, mejorando el flujo de trabajo para que el/la usuario/a pueda generar un modelo de forma eficaz, eficiente y bajo una experiencia lo más amigable posible, sin interferencias ni distracciones.

Así, por defecto, tan solo aparecen los menús, opciones y parámetros para el cálculo de la hidrodinámica, que es el módulo básico para el resto de procesos que incorpora Iber. Como se puede ver en la Fig. 3a, el menú ‘Datos’ se ha simplificado enormemente mostrando, únicamente, aquellos submenús relacionados con los parámetros generales del modelo (‘Datos del Problema’), la ‘Hidrodinámica’ (condiciones de contorno, iniciales, internas, etc.) y aquellos parámetros que guardan relación con los esfuerzos tensionales (rugosidad del fondo, turbulencia y viento). En la Fig. 3b se muestra la ventana ‘Datos del Problema’ que, al igual que el menú ‘Datos’, por defecto tan solo muestra opciones generales del módulo hidrodinámico.

Fig. 3. (a) Disposición por defecto del menú ‘Datos’. (b) Ventana de ‘Datos del Problema’ por defecto.

Asimismo, ha habido pequeños cambios en la interfaz que también modifican ligeramente el flujo de trabajo (como se verá en el apartado 2.1.2). Es el caso de la herramienta ‘Brecha’, antes ubicada en el menú ‘Datos’, que ha sido trasladada al menú ‘Herramientas Iber’.

2.1.2 Flujo de trabajo

El flujo de trabajo de la versión 3 de Iber sigue siendo el mismo que en versiones anteriores, pero ahora si un/a usuario/a desea realizar cálculos con alguno de los módulos adicionales a la hidrodinámica, lo tendrá que activar. Para ello, deberá emplear la ventana de Plug-ins (Fig. 4a), que es accesible desde Herramientas Iber >> Plug-ins de la barra de herramientas de Iber (Fig. 4b).

Fig. 4. (a) Ventana de selección de plug-ins. (b) Activación de los plug-ins mediante el atajo en la barra de herramientas de Iber.

Esta ventana permite a usuarios/as seleccionar, en primera instancia, el método de paralelización cálculo mediante tarjetas gráficas (GPU) y, en segundo lugar, el módulo de cálculo deseado. La elección de un método de paralelización del cálculo mediante GPUs es opcional y, por el momento, está condicionado al tipo o tipos de módulos de cálculo que se desee emplear.

Fig. 5. Ejemplo de activación del módulo de ‘Transporte de Sedimentos’: (a) cambios en la ventana de activación de Plug-ins; (b) incorporación de los menús y submenús en ‘Datos’; (c) actualización de la ventana ‘Datos del Problema’ con las nuevas opciones (c).

Una vez seleccionado el módulo de cálculo, se aplicarán los cambios a la interfaz para mostrar las funcionalidades de dicho módulo. A modo de ejemplo se selecciona el módulo de ‘Transporte de Sedimentos’. De este modo, automáticamente se desactivarán los módulos que no son compatibles con el que se ha activado (Fig. 5a) (esto de describe en el apartado 2.2). Además, el menú ‘Datos’ (Fig. 5b) y la ventana ‘Datos del Problema’ (Fig. 5c) ahora incorporarán las opciones, parámetros y funcionalidades propias de dicho módulo. El flujo de trabajo a partir de aquí es idéntico al de versiones anteriores.

2.2 Compatibilidad entre módulos de cálculo

Iber es un modelo en continuo desarrollo en el que, a lo largo de su trayectoria, ha ido incorporando nuevos módulos de cálculo y procesos que amplían y agilizan la computación de procesos hidrodinámicos, así como otros íntimamente relacionados con el cálculo hidrodinámico en lámina libre.

Si bien todos los módulos de cálculo dependen de la hidrodinámica y, por tanto, están completamente integrados con ella, no todos son compatibles entre sí. Esto sucede, por ejemplo, con el módulo de ‘Hábitat’, el cual no es compatible con el módulo de ‘Drenaje Urbano’ ni con el de ‘Erosión de suelos’ porque el hábitat fluvial no guarda relación directa con el drenaje urbano y los procesos de erosión de suelos en cuencas. Otro caso sería el de ‘Transporte de Sedimentos’, que tampoco es compatible con los dos módulos anteriores. El motivo por el cual no es compatible con ‘Erosión de suelos’ es porque este módulo ya incorpora su propio submódulo de transporte de sedimentos; mientras que no es compatible con ‘Drenaje Urbano’ porque todavía no se ha integrado con este módulo.

Asimismo, desde la versión 2.5, Iber incorpora la posibilidad de realizar el cálculo hidrodinámico y de procesos hidrológicos mediante dos posibilidades: código secuencial parcialmente paralelizado mediante CPU; y cálculo paralelo mediante GPU (tarjetas gráficas) con arquitectura CUDA. El cálculo con CPU incorpora todos y cada uno de los últimos módulos de cálculo y funcionalidades, mientras que el cálculo con GPU incorpora algunos de los módulos de cálculo y funcionalidades del cálculo en CPU. Es por ello que cuando se activa en Plug-ins la opción de paralelización ‘Plus’ o ‘R-Iber’, este último recientemente incorporado fruto del presente desarrollo, no es posible seleccionar diversos módulos de cálculo ni tampoco utilizar otras funcionalidades que sí están disponibles en la versión secuencial.

2.3 Nuevas herramientas y funcionalidades mejoradas

Otros cambios relacionados con el flujo de trabajo, pero que a la vez se trata de nuevas funcionalidades y herramientas, es la consideración de flujo a presión en tramos cubiertos y la eliminación de la herramienta ‘Brecha’ del menú ‘Datos’. De este modo, para ‘Cubiertas’, el/la usuario/a deberá proceder como en versiones anteriores, imponiendo la cubierta donde considere, pero ahora deberá activar su cálculo a través de su correspondiente pestaña en ‘Datos del Problema’ (Fig. 3b). Aquí tendrá, además, la posibilidad de calcular cubiertas mediante la metodología tradicional de la ranura de Preissmann 2D [11,12] o el nuevo método basado en la teoría de cubiertas elásticas [13]. En lo que respecta a la herramienta ‘Brecha’ (Fig. 6a), así como pasa con ‘Alcantarillas’ (Fig. 6b), el/la usuario/a no tiene que activarla, tan solo tiene que definirla como se ha hecho hasta ahora en versiones anteriores de Iber y decidir si desea o no calcularla.

Fig. 6. Integración de la opción de cálculo en: (a) ‘Brechas’; (b) ‘Alcantarillas’.

Asimismo, hay un ligero cambio en cómo el/la usuario/a solicita los resultados hidrodinámicos, así como aquellos que dependen de la hidrodinámica. Éstos se siguen indicando desde la pestaña ‘Resultados’ de la ventana ‘Datos del Problema’, sin embargo, en la versión 3 de Iber, se han ordenado y agrupado para una mayor agilidad y comodidad tal como muestra la Fig. 7a. En la opción ‘Hidrodinámica’ se encuentran resultados propios de la resolución de las ecuaciones de flujo en lámina libre y que, por defecto, están casi todos activados. Los términos de ‘Peligrosidad’ se han separado según los criterios del RD 9/2008 [14] y del ACA 2003 [15], agilizando así el cálculo, y, además, la pestaña ‘Peligrosidad personalizada’ se ha integrado en este apartado. Asimismo, se ha creado una nueva opción para extraer resultados, llamada ‘Sondas’, que permite extraer valores hidrodinámicos (calado, velocidad, etc.) de puntos concretos del modelo a intervalos de tiempo diferentes a los de la escritura de resultados (Fig. 7b). Esta opción es especialmente útil para modelos cuyo tiempo de cálculo es muy extenso y en los que no es necesario escribir todos los resultados con pasos de tiempo reducidos.

Fig. 7. Organización de la pestaña ‘Resultados’ y nuevas opciones: (a) vista general; (b) nueva herramienta ‘sondas’.

En cuanto al post-proceso, la principal novedad está en una mejora sustancial del tiempo de respuesta y de cálculo para con la generación de rasters de resultados (Herramientas Iber >> Resultados a Raster). Además, en la actual versión existe la posibilidad de cambiar el tipo de análisis, seleccionar en qué instante de tiempo se desea extraer el resultado o resultados, y se han mejorado las opciones de exportación para que sean más rápidas y eficientes. De este modo, un/a usuario/a puede seleccionar uno o varios resultados de un análisis concreto, exportarlos en el paso de tiempo que está activo, en todos o en uno concreto que no se esté visualizando, con un tipo u otro de interpolación y un tamaño de celda determinado. Los resultados se guardarán dentro de la carpeta del modelo preservando la clasificación por medio del tipo de análisis y paso de tiempo. Por ejemplo, si se seleccionan ‘Velocidad’ y ‘Cota de agua’ (Fig. 8a) y se aplica, en la carpeta del modelo se crearán los n rasters de velocidad y cota de agua asociados a los n pasos de tiempo analizados (Fig. 8b).

Fig. 8. (a) Ejemplo de la herramienta ‘Exportar rasters de resultados’. (b) conjunto de rasters de resultados exportados.

Finalmente, para el módulo de Transporte de Sedimentos para mezclas se ha diseñado una nueva funcionalidad que permite conocer la distribución granulométrica de la muestra. Para ello, es necesario ir a Herramientas Iber >> Transporte de Sedimentos >> Granulometría. La Fig. 9 muestra, a modo de ejemplo, la granulometría en un instante de tiempo concreto en los elementos 23, 28 y 33. Como se observa, la gráfica está en escala logarítmica.

Fig. 9. Visualización de la distribución granulométrica de una muestra.

3 Módulo de drenaje urbano

3.1 Fundamentos teóricos

El flujo de agua en zona urbana presenta una complejidad considerable respecto al flujo en zona no urbana. Por un lado, el agua puede circular por la superficie (sistema de drenaje mayor) y/o por la red de colectores (sistema de drenaje menor). Por otro lado, el flujo puede ser en lámina libre o a presión. En el sistema de drenaje menor, en los colectores y en las uniones de una red de alcantarillado ocurren fenómenos hidráulicos complejos. La interacción entre el sistema de drenaje mayor (tejados, calles, áreas verdes, áreas grises, canales naturales y canales artificiales, etc.) y el sistema de drenaje menor (colectores, uniones, compuertas, etc.), conocida como drenaje urbano dual, se realiza en los puntos físicos correspondientes a la localización espacial de los elementos de captación (fuentes y sumideros).

Desde el punto de vista numérico, para la resolución de las ecuaciones de flujo en lámina libre en una dimensión, las ecuaciones de flujo en presión (tanto en una dimensión como en dos) y las ecuaciones de aguas poco profundas en dos dimensiones, se emplean esquemas numéricos en volúmenes finitos, capaces de calcular discontinuidades y frentes de onda en la solución, sin necesidad de técnicas adicionales [16].

3.1.1 Ecuaciones de la hidrodinámica

Para el flujo en zona urbana se utilizan las ecuaciones de Saint Venant, tanto en una dimensión como en dos dimensiones. Las ecuaciones bidimensionales se utilizan sobre todo en superficie, pero también para elementos que forman parte de la red de colectores, como pueden ser depósitos de retención o uniones de colectores. Las ecuaciones unidimensionales se utilizan principalmente en los colectores.

3.1.2 Ecuaciones de Saint Venant en 2D

En el módulo de drenaje urbano, los cálculos en dos dimensiones se realizan en base a las ecuaciones de aguas someras, o ecuaciones de Saint Venant 2D, que se describen en detalle en el Manual de Referencia Básico. Se repiten aquí por completitud, incluyendo los términos que representan las componentes de precipitación e infiltración y con la fórmula de Manning para representar la fricción, siendo el sistema de ecuaciones el siguiente:

${\begin{matrix}{\frac {\partial h}{\partial t}}+{\frac {\partial {q}_{x}}{\partial x}}+{\frac {\partial {q}_{y}}{\partial y}}=R-f\\{\frac {\partial {q}_{x}}{\partial t}}+{\frac {\partial }{\partial x}}\left({\frac {{q}_{x}^{2}}{h}}+g{\frac {{h}^{2}}{2}}\right)+{\frac {\partial }{\partial y}}\left({\frac {{q}_{x}{q}_{y}}{h}}\right)=-gh{\frac {\partial {z}_{b}}{\partial x}}-g{\frac {{n}^{2}}{{h}^{\frac {7}{3}}}}\left|{\mathit {\boldsymbol {q}}}\right|{q}_{x}\\{\frac {\partial {q}_{y}}{\partial t}}+{\frac {\partial }{\partial x}}\left({\frac {{q}_{x}{q}_{y}}{h}}\right)+{\frac {\partial }{\partial y}}\left({\frac {{q}_{y}^{2}}{h}}+g{\frac {{h}^{2}}{2}}\right)=-gh{\frac {\partial {z}_{b}}{\partial y}}-g{\frac {{n}^{2}}{{h}^{\frac {7}{3}}}}\left|{\mathit {\boldsymbol {q}}}\right|{q}_{y}\end{matrix}}$

(1)

En ellas, ${\textstyle h}$ es el calado de agua, ${\textstyle {q}_{x}}$ y ${\textstyle {q}_{y}}$ son las dos componentes del caudal unitario, ${\textstyle \vert {\mathit {\boldsymbol {q}}}\vert }$ es su módulo, ${\textstyle R}$ es la intensidad de precipitación, ${\textstyle f}$ es la tasa de infiltración, ${\textstyle {z}_{b}}$ es la cota del fondo, ${\textstyle g}$ es la aceleración de la gravedad y ${\textstyle n}$ es el coeficiente de Manning.

Estas ecuaciones son las que se resuelven mediante el módulo hidrodinámico de Iber desde sus inicios, y se encuentran ampliamente documentadas, tanto las propias ecuaciones como los métodos de resolución (p.ej. en Toro [17]).

3.1.3 Ecuaciones de Saint Venant en 1D

La deducción de las ecuaciones de Saint Venant en 1D se realiza directamente aplicando las leyes de conservación de la masa y cantidad de movimiento a un volumen de control considerando una sección arbitraria y canal no prismático [18]. El resultado es un sistema de ecuaciones que, en forma conservativa, se puede escribir como sigue:

${\begin{matrix}{\frac {\partial A}{\partial t}}+{\frac {\partial Q}{\partial x}}=0\\{\frac {\partial Q}{\partial t}}+{\frac {\partial }{\partial x}}\left({\frac {{Q}^{2}}{A}}+g{I}_{1}\right)=g{I}_{2}+gA\left({S}_{0}-{S}_{f}\right)\end{matrix}}$

(2)

En las expresiones anteriores, ${\textstyle A}$ es el área de la sección mojada, ${\textstyle Q}$ el caudal circulante, ${\textstyle {S}_{0}}$ la pendiente del terreno, ${\textstyle {S}_{f}}$ la pendiente de fricción o motriz, ${\textstyle {I}_{1}}$ la fuerza debida a la presión del agua en una sección e ${\textstyle {I}_{2}}$ la contribución de las fuerzas de presión del contorno. Las tres últimas se expresan respectivamente como:

${S}_{f}={\frac {{v}^{2}{n}^{2}}{{R}_{h}^{4/3}}}$	(3)
${I}_{1}=\int _{0}^{h}\left(h-\eta \right)b\left(x,\eta \right)d\eta$	(4)
${I}_{2}=\int _{0}^{h}\left(h-\eta \right){\frac {\partial b\left(x,\eta \right)}{\partial x}}d\eta$	(5)

donde ${\textstyle n}$ es el coeficiente de fricción de Manning, ${\textstyle {R}_{h}^{}}$ el radio hidráulico. Para canales prismáticos, aunque tengan una sección cualquiera, el término ${\textstyle {I}_{2}}$ es idénticamente igual a cero.

3.1.4 Esquema numérico para el flujo en colectores

La aplicación de un método numérico conlleva la discretización espacial del dominio en celdas o elementos de cálculo. La unión de estos elementos forma una malla. Iber, en todos sus módulos, utiliza el Método de los Volúmenes Finitos para resolver las ecuaciones, tal como se detalla en el Manual de Referencia Hidráulico, por lo que cada elemento será a su vez un Volumen Finito. Un volumen finito unidimensional es un segmento, mientras que un volumen finito bidimensional es una superficie (en Iber tienen forma de triángulo o cuadrilátero).

En el módulo de drenaje urbano conviven dominios unidimensionales y bidimensionales y, por lo tanto, también habrá mallas en 1D y mallas en 2D (Fig. 10). La malla puede ser estructurada (Fig. 10b) o no estructurada (Fig. 10c). En las mallas estructuradas, cada elemento de la misma está inequívocamente identificado por unos índices en coordenadas cartesianas. Por su parte, en mallas no estructuradas los elementos no tienen un orden particular. Los elementos de una malla estructurada en 2D por lo general son cuadriláteros, aunque también pueden ser triángulos equiláteros. En cambio, en una malla no estructurada pueden ser triángulos, cuadriláteros o una combinación de éstos. Las mallas no estructuradas tienen la ventaja de adaptarse mejor a geometrías complejas como las que presentan la batimetría y fronteras de estudios hidráulicos. Sin embargo, éstas requieren de una estructura de datos compleja y unos requerimientos de memoria mayores.

Fig. 10. Tipos de malla: (a) malla en una dimensión; (b) malla estructurada en dos dimensiones; (c) malla no estructurada en dos dimensiones.

En los dominios unidimensionales, la discretización empleada para la resolución en el tiempo de las ecuaciones de Saint Venant con el método de los volúmenes finitos se puede escribir de forma genérica como:

${\mathit {\boldsymbol {U}}}_{i}^{n+1}={\mathit {\boldsymbol {U}}}_{i}^{n+1}-{\frac {\Delta t}{\Delta x}}\left({\mathit {\boldsymbol {F}}}_{i+{\frac {1}{2}}}^{\ast }-\right.$$\left.{\mathit {\boldsymbol {F}}}_{i-{\frac {1}{2}}}^{\ast }\right)+{\frac {\Delta t}{\Delta x}}{\mathit {\boldsymbol {H}}}_{i}^{\ast }$

(6)

En esta nomenclatura, los subíndices enteros hacen referencia a un volumen finito o elemento de la malla (segmento asociado a una sección transversa), mientras que los subíndices ${\textstyle i+}$${\frac {1}{2}}$ e ${\textstyle i-{\frac {1}{2}}}$ hacen referencia a los puntos de contacto entre segmentos. En Iber el cálculo de los flujos numéricos se realiza con el esquema de Roe [19] y su expresión resulta:

${\mathit {\boldsymbol {F}}}_{i+{\frac {1}{2}}}^{\ast }={\frac {1}{2}}\left({\mathit {\boldsymbol {F}}}_{i}+\right.$$\left.{\mathit {\boldsymbol {F}}}_{i+1}\right)-{\frac {1}{2}}\left[\sum _{j=i}^{2}{\tilde {\alpha }}_{j}{\varphi }_{j}{\tilde {e}}_{j}+\right.$$\left.\sum _{j=i}^{2}{\tilde {\gamma }}_{j}signo\left({\tilde {\lambda }}_{j}\right){\tilde {\mathit {\boldsymbol {e}}}}_{j}\right]$

(7)

En la ecuación anterior, y las siguientes, los términos con el superíndice ${\textstyle {\tilde {}}}$ son valores promediados en los contornos de los elementos. Los coeficientes constantes o fuerzas de cada onda ${\textstyle {\tilde {\alpha }}_{j}}$, los valores y vectores propios ${\textstyle {\tilde {\lambda }}_{j}}$ y ${\textstyle {\tilde {\mathit {\boldsymbol {e}}}}_{j}}$ de la matriz jacobiana definida como ${\textstyle {\mathit {\boldsymbol {A}}}=}$$\partial {\mathit {\boldsymbol {F}}}/\partial {\mathit {\boldsymbol {U}}}$, se pueden expresar respectivamente como:

${\tilde {\alpha }}_{1,2}={\frac {\left({\tilde {c}}\mp {\tilde {u}}\right)\Delta A\pm \Delta Q}{2{\tilde {c}}}}$	(8)
${\tilde {\lambda }}_{1,2}={\tilde {c}}\mp {\tilde {u}}$	(9)
${\tilde {\mathit {\boldsymbol {e}}}}_{1,2}=\left[{\begin{matrix}1\\{\tilde {c}}\pm {\tilde {u}}\end{matrix}}\right]$	(10)

siendo,

$\Delta A={A}_{i+1}-A$

(11)

$\Delta Q={Q}_{i+1}-Q$	(12)
${\tilde {u}}={\frac {{\frac {{Q}_{i}}{\sqrt {{A}_{i}}}}+{\frac {{Q}_{i+1}}{\sqrt {{A}_{i+1}}}}}{{\sqrt {{A}_{i}}}+{\sqrt {{A}_{i+1}}}}}$	(13)

La corrección de entropía de Harten y Hyman ${\textstyle {\varphi }_{j}}$, según Toro [17], es:

${\varphi }_{j}=\left\{{\begin{matrix}\left|{\tilde {\lambda }}_{j}\right|\,si\,\left|{\tilde {\lambda }}_{j}\right|\geq {\epsilon }_{j}\\{\epsilon }_{j}\,si\,\left|{\tilde {\lambda }}_{j}\right|<{\epsilon }_{j}\end{matrix}}\right.$

(14)

donde,

${\epsilon }_{j,i+1/2}=max\left[0,\left({\tilde {\lambda }}_{j,i+1/2}-{\tilde {\lambda }}_{j,i}\right),\left({\tilde {\lambda }}_{j,i+1}-\right.\right.$$\left.\left.{\tilde {\lambda }}_{j,i+1/2}\right)\right]$

(15)

Para reproducir correctamente el salto de las fuerzas de presión, se utiliza una expresión que mantiene el significado físico de la celeridad, entendido como la variación de las fuerzas de presión en cada sección transversal con respecto a la variación del área de flujo en cada sección [20]:

${\tilde {c}}^{2}=\left\{{\begin{matrix}g{\frac {\left({I}_{i+1}-{I}_{i}\right)-{\left({\Delta {I}_{i}\vert }_{\tilde {A}}\right)}_{i+1/2}}{{A}_{i+1}-{A}_{i}}}\,si\,{A}_{i}\not ={A}_{i+1}\\{\frac {{c}_{i}^{2}+{c}_{i+1}^{2}}{2}};\,o\,bien\,g{\frac {\tilde {A}}{\tilde {B}}}\,;{\acute {o}}\,bien\,{\frac {\Delta I}{\Delta A}}\,si\,{A}_{i}={A}_{i+1}\end{matrix}}\right.$

(16)

Cuando ${\textstyle {A}_{i}={A}_{i+1}}$, Iber emplea la segunda opción. En las anteriores expresiones:

${\tilde {A}}={\frac {2{A}_{i}{A}_{i+1}}{{A}_{i}+{A}_{i+1}}}$	(17)
${\tilde {B}}={\frac {{B}_{i}+{B}_{i+1}}{2}}$	(18)

La variación de las fuerzas de presión ${\textstyle {\left({\Delta {I}_{i}\vert }_{\tilde {A}}\right)}_{i+1/2}}$ en volúmenes finitos para un área constante ${\textstyle {A}_{i+1/2}}$ es:

${\left({\Delta {I}_{1}\vert }_{\tilde {A}}\right)}_{i+1/2}={\left({{I}_{i}\vert }_{\tilde {A}}\right)}_{i+1}-$${\left({\Delta {I}_{1}\vert }_{\tilde {A}}\right)}_{i}$

(19)

donde ${\textstyle {\left({{I}_{i}\vert }_{\tilde {A}}\right)}_{i}}$ y ${\textstyle {\left({{I}_{i}\vert }_{\tilde {A}}\right)}_{i+1}}$ son las fuerzas de presión en el volumen finito ${\textstyle i}$ e ${\textstyle i+}$$1$ para un valor del área ${\textstyle {\tilde {A}}}$.

Finalmente, el término ${\textstyle {\tilde {\gamma }}_{j}}$ se define como:

${\begin{matrix}\,{\tilde {\gamma }}_{1}={\frac {1}{2{\tilde {c}}}}g\Delta \left({{I}_{i}\vert }_{\tilde {A}}\right)\\{\tilde {\gamma }}_{2}=-{\tilde {\gamma }}_{1}\end{matrix}}$

(20)

3.1.5 Flujo mixto en conductos cerrados

La ocurrencia simultánea de flujo en lámina libre y flujo en presión en un mismo conducto se conoce comúnmente como flujo mixto. Este fenómeno se puede encontrar, frecuentemente, en los colectores de aguas pluviales, en túneles y conductos de obras de toma de centrales hidroeléctricas, llenado/vaciado de conductos, colectores de almacenamiento, etc.

El proceso de transición de flujo en lámina libre a flujo en presión puede ser inducido, entre otras causas, por la variación de caudales de entrada, bloqueo o reducción de la capacidad del conducto, sumersión de la salida del conducto, operación de elementos de control (compuertas, vertederos, etc.), fallo de estaciones de bombeo, presencia de pozos de caída y falta de aire, y puede producir una variación importante de los calados, cargas de presión, velocidades, entrada o expulsión de aire.

Fig. 11. Flujo mixto en un conducto: (a) transición gradual; (b) transición brusca.

Los frentes de onda generados en la transición de flujo en lámina libre a flujo en presión pueden ser de dos tipos: gradual, que consiste de un frente de onda en lámina libre con condiciones de flujo en lámina libre en ambos lados del frente (Fig. 11a); o brusco, que es un frente de onda a tubo lleno con condiciones de flujo en lámina libre delante y flujo en presión atrás del mismo (Fig. 11b). Además, estos frentes de onda de presión se pueden clasificar como: interfase positiva (el frente de onda de presión se mueve hacia el flujo en lámina libre), interfase negativa (el frente de onda de presión se mueve hacia el flujo en presión), o interfase con cambio de dirección (el frente de onda se mueve de una interfase positiva hacia una interfase negativa o viceversa).

Fig. 12. Método de la ranura de Preissmann: (a) flujo en lámina libre; (b) flujo en presión; (c) flujo en presión negativa.

En el módulo de drenaje urbano en Iber, para el cálculo de los flujos a presión, tanto en una dimensión como en dos, se ha implementado el método de la ranura de Preissmann [21]. El método de la ranura de Preissmann consiste en implementar una ranura hipotética de ancho ${\textstyle {T}_{s}}$ en la clave de los conductos cerrados (Fig. 12b), proporcionando la presión hidrostática extra mientras que se permite el uso de las ecuaciones de flujo en lámina libre para ambos tipos de flujo (flujo en lámina libre y flujo en presión). El ancho de la ranura se calcula bajo la consideración de que la celeridad de la onda de gravedad ${\textstyle c}$ es igual a la celeridad de la onda de presión ${\textstyle a}$. Esencialmente, esta aproximación explota la similitud de la ecuación de onda del flujo en lámina libre y flujo en presión. Así, a partir de la definición de la celeridad de la onda en lámina libre, se tiene:

$c={\sqrt {g{\frac {A}{B}}}}$

(21)

Despejando el ancho de la superficie libre del agua de la expresión anterior, que al emplear el método de la ranura de Preissmann se simboliza con ${\textstyle {T}_{s}}$ (Fig. 12b), utilizando el área del colector a tubo lleno denominada área máxima ${\textstyle {A}_{max}}$ y empleando la celeridad de la onda de presión ${\textstyle a}$, se obtiene:

${T}_{s}={\frac {g{A}_{max}}{{a}^{2}}}$

(22)

Cuando en un colector ocurre flujo en presión, el calado calculado es la carga de presión. Así, con el método de la ranura de Preissmann clásico, el área en presión es ${\textstyle {A}_{p}=A-}$${A}_{max}$ y la carga debida a la presión ${\textstyle {h}_{p}}$ se obtiene como:

${h}_{p}={\frac {{A}_{p}}{{T}_{s}}}$

(23)

Para mitigar los problemas relacionados con las inestabilidades numéricas debidas a los esquemas numéricos cuando el flujo cambia rápidamente de flujo en lámina libre a flujo en presión, se plantea el uso de una transición geométrica gradual entre el conducto y el ancho de la ranura de Preissmann. Barnett [22] sugiere una transición geométrica gradual de tipo exponencial, la cual ha sido modifican ligeramente, resultando:

$w={T}_{s}{e}^{k\left(1,5{h}_{max}-h\right)}$

(24)

donde ${\textstyle h}$ es el calado o carga de presión, ${\textstyle {h}_{max}}$ es el calado máximo, ${\textstyle w}$ el ancho de la ranura en la zona de la transición geométrica, ${\textstyle k}$ un parámetro que determina la forma de la transición; valores más pequeños corresponden a transiciones más elongadas. En Iber se emplea la transición geométrica gradual a partir de una carga de presión de ${\textstyle h=}$$0,98{h}_{max}$ para conductos con sección transversal convergente en su parte superior y de ${\textstyle h=}$${h}_{max}$ para otro tipo de secciones transversales, ambos hasta un calado de ${\textstyle \,h=}$$1,5{h}_{max}$.

La principal limitación de este método es que asume flujo a superficie libre en todo el sistema. Por consiguiente, si la carga piezométrica cae por debajo de la coronación del conducto, se desarrolla flujo en lámina libre y, por tanto, sería incapaz de simular flujo en presión negativa. Con el fin de simular flujo en presión subatmosférica, Kerger et al. [23] ampliaron este concepto y desarrollaron la ranura de Preissmann negativa que consiste en extender la ranura por debajo de la coronación del conducto (Fig. 12c). De esta forma, cuando el flujo está en presión, para cada valor de calado menor al calado máximo existen dos valores de área: uno para superficie libre y otro para flujo en presión negativa.

Cabe decir que el método de la ranura de Preissmann podría ser numéricamente inestable cuando el ancho de la ranura es pequeña, necesaria para representar grandes niveles de agua o una celeridad de onda en presión importante. En Iber se considera correctamente la calidad de la onda de presión en la reducción del incremento de tiempo por la condición de Courant.

3.1.6 Uniones de colectores y dominios 1D y 2D

En una red de colectores puede haber zonas adecuadas para simular en 1D, pero otras que conviene simular en 2D, como pueden ser depósitos de retención, cámaras de disipación de energía, uniones complejas en colectores, entre otras. Por ello se ha implementado una formulación para la conexión entre dominios calculados en 1D y dominios calculados en 2D.

Para uniones entre dominios 1D y 2D se utiliza el método completamente conservativo presentado por Bladé et al. [16]. En él, una unión ${\textstyle k}$ se puede discretizar en ${\textstyle i}$ elementos o volúmenes finitos en 2D y ${\textstyle j}$ elementos frontera o ficticios; cada uno de los ${\textstyle i}$ colectores ( ${\textstyle p}$ colectores de entrada y ${\textstyle q}$ colectores de salida) se puede descomponer en ${\textstyle n}$-2 celdas o volúmenes finitos en 1D y 2 celdas frontera o ficticios (volúmenes finitos 1 y ${\textstyle n}$). La conexión entre los dominios en 1D (conductos) y en 2D (unión) se realiza a través de los volúmenes finitos ubicados en la frontera de cada dominio, para lo cual existe un conjunto de elementos (volúmenes finitos en 2D que representan la unión) conectados a una celda (volumen finito en 1D que representa un tramo de colector, cuyo centro alberga los valores de las variables hidráulicas y las propiedades geométricas de la sección transversal del conducto).

3.1.7 Uniones de colectores

Las uniones son los elementos que conectan los tramos de ríos, las calles en una zona urbana, los canales en una red de riego o los colectores en una red de alcantarillado. En una red de alcantarillado, un pozo de registro al inicio de un colector se puede entender como una unión con un único colector de salida. Este tipo de unión, además, puede recibir un volumen de agua del escurrimiento superficial. Uniones comunes son aquellas formadas por dos, tres o cuatro conductos, pero también existen uniones con cinco o más conductos.

La unión de colectores presenta una problemática específica, ya que en una unión pueden concurrir colectores con alineaciones verticales y horizontales muy diferentes, producirse cambios de sección y pendiente, así como notables variaciones de caudal entre colectores de entrada y salida [24]. Estas singularidades geométricas del colector y de las uniones, pueden provocar que en los colectores se presenten múltiples combinaciones de los diferentes tipos de flujo.

Para el cálculo de las uniones de colectores simples (es decir, donde no se requiere la caracterización detallada del flujo en 2D), la estrategia implementada en Iber es considerar la unión como un único elemento 2D. Ello se hace internamente, sin información al usuario, y la única opción que existe es si se considera continuidad de flujo a través de la unión, o no. En el caso de considerar continuidad de flujo, que sería para el caso donde los colectores de entrada están bien alineados, los caudales de entrada y salida están relacionados por la conservación de la cantidad de movimiento. En caso contrario, que correspondería a la falta de alineación entre colectores de entrada y salida, la presencia de obstáculos, o directamente de disipadores de energía, internamente se considera que la velocidad en el colector es nula.

3.1.8 Drenaje dual

El drenaje urbano integra dos componentes distintos: uno superficial, o sistema de drenaje mayor, compuesto por calles, cunetas, y canales naturales o artificiales; y un sistema subsuperficial, o sistema de drenaje menor, compuesto por colectores, uniones y estructuras de control como compuertas, vertederos, etc. (Fig. 13). El sistema de drenaje menor habitualmente se diseña para funcionar en lámina libre para un caudal asociado a un periodo de retorno, relacionado con los eventos más frecuentes. Cuando se presenta un evento de lluvia con periodo de retorno superior al de diseño, los colectores (sistema de drenaje menor) pueden entrar en carga, incluso el agua puede llegar a salir a la superficie urbana y ser trasportada por el sistema de drenaje mayor (calles). La conexión entre los dos sistemas es bidireccional, es decir, el flujo puede ir del sistema de drenaje mayor al sistema de drenaje menor, y viceversa.

Fig. 13. Elementos y procesos hidráulicos en el drenaje urbano ante un evento de lluvia.

Los sumideros se diseñan para recolectar el agua de lluvia que circula sobre las calles, así como de otras áreas urbanas, y conducirla hacia los colectores, o lo que es lo mismo, del sistema de drenaje mayor al sistema de drenaje menor. En ocasiones puede existir conducción del flujo del sistema de drenaje menor al sistema de drenaje mayor, cuando la carga piezométrica en el colector es mayor a la superficie libre del agua en la calle en el punto donde se ubica dicho elemento. Los elementos que funcionan como fuentes son, por lo general, las tapas de los pozos de registro, pero los sumideros también pueden desempeñar dicha función.

Según lo descrito anteriormente, hay una interacción entre el sistema de drenaje mayor y el sistema de drenaje menor, que da lugar a un intercambio de flujo en ambas direcciones a través de las fuentes y los sumideros, el cual debe ser estimado adecuadamente para modelar el drenaje urbano dual convenientemente. De esta forma, en los puntos de localización de las fuentes y sumideros el proceso de interacción de flujo da lugar a tres casos:

a) La carga piezométrica en el colector es menor que el nivel del terreno de la superficie urbana, por lo que ésta no interviene sobre el flujo de intercambio y el flujo captado del sistema de drenaje mayor entra al sistema de drenaje menor libremente (Fig. 14a).

b) La carga piezométrica del colector se encuentra entre el nivel del terreno de la superficie urbana y la superficie libre del agua en la calle, por lo que el flujo aún entra al sistema de drenaje menor, pero bajo la influencia de la carga hidráulica del colector (Fig. 14b).

c) La carga hidráulica en el colector es mayor a la superficie libre del agua en la calle, por lo tanto, el flujo escapa del sistema de drenaje menor hacia el sistema de drenaje mayor (Fig. 14c).

Fig. 14. Esquema de flujo en los puntos de interacción (fuente/sumidero): (a) ‘sumidero libre’; (b) ‘sumidero anegado’; (c) ‘fuente’.

La capacidad hidráulica o de intercepción (cantidad de agua interceptada por un sumidero bajo ciertas condiciones) gobierna el caudal de intercambio (caudal interceptado por los sumideros) en los casos a) y b). El caudal interceptado es función de la posición, tamaño y tipo de reja, condiciones de flujo, rugosidad de la superficie, geometría de la calzada, y el fenómeno de colmatación [25,26]. Esta definición, válida para los casos a) y b), asume que la funcionalidad de los sumideros no es afectada por las condiciones existentes en los pozos de registro y en los colectores.

Para estimar los caudales de intercambio, existen formulaciones semiempíricas tipo orificio o vertedero y formulaciones experimentales basadas en la definición de la eficiencia hidráulica, las características del flujo y la geometría de los sumideros. La primera formulación se puede utilizar en los casos a), b) y c), mientras la segunda, solamente es aplicable al caso a).

Formulación orificio/vertedero

Este tipo de formulación considera que el flujo de intercambio es función del calado de aproximación. Así, el caudal de intercambio para el caso a) puede ser estimado mediante una formulación tipo vertedero para calados pequeños y una formulación tipo orificio para calados grandes. Como la transición de los regímenes de flujo de vertedero a orificio y viceversa no es muy clara, se recomienda utilizar la formulación que estime el caudal mínimo. En el caso b), el flujo puede ser caracterizado con una formulación tipo orificio o tipo compuerta. Finalmente, para el caso c) una ecuación tipo orificio es la más apropiada.

Las fórmulas tradicionales para estimar el caudal a través de un orificio ${\textstyle {Q}_{o}}$ y a través de un vertedero ${\textstyle {Q}_{v}}$ son:

${Q}_{o}={C}_{d,o}{A}_{int}{\sqrt {2gH}}$	(25)
${Q}_{v}={C}_{d,v}{L}_{int}{H}^{3/2}$	(26)

donde ${\textstyle {A}_{int}}$ es el área de captación, ${\textstyle {L}_{int}}$la longitud de captación, ${\textstyle H}$ la carga hidráulica, ${\textstyle {C}_{d,o}}$ el coeficiente de descarga del orificio y ${\textstyle {C}_{d,v}}$ el coeficiente de descarga del vertedero. Valores entre 0,60 y 0,80 son recomendados para el coeficiente de descarga del orificio y valores entre 1,70 y 3,20 para el coeficiente de descarga del vertedero.

Uno de los inconvenientes del cálculo del caudal a través de un vertedero consiste en determinar la longitud de intercepción correcta a emplea. En cambio, una formulación tipo orificio emplea el área ${\textstyle {A}_{int}}$ (corresponde al área de huecos). Cuando el intercambio de flujo es del sistema de drenaje mayor al sistema de drenaje menor y el calado cubre totalmente el elemento de captación, esta longitud se caracteriza con el perímetro total del elemento. En caso contrario, esta longitud será menor al perímetro total del elemento de captación.

Metodología UPC

Esta metodología permite estimar el caudal interceptado por un sumidero a partir su eficiencia hidráulica ${\textstyle E}$, definida como la proporción o porcentaje entre el caudal interceptado ${\textstyle {Q}_{int}}$ y el caudal total de aproximación al sumidero ${\textstyle Q}$.

$E={\frac {{Q}_{int}}{Q}}$

(27)

Los resultados de los ensayos experimentales, realizados en el laboratorio de hidráulica del Departament d’Enginyeria Hidràulica, Marítima i Ambiental (DEHMA) de la Escola Tècnica Superior d’Enginyers de Camins, Canals i Ports (ESTCCPB) de la Universitat Politècnica de Catalunya (UPC), para analizar el comportamiento hidráulico y capacidad de captación de diferentes tipos de rejas y macro-rejas [26], determinaron una ecuación de tipo potencial que relaciona la eficiencia de captación de rejas y macro-rejas con el calado y caudal de aproximación a la reja y de dos parámetros característicos de la geometría las mismas, definida como:

$E=A{\left({\frac {Q}{h}}\right)}^{-B}$

(28)

donde,

$A={\frac {1,988{A}_{g}^{0,403}}{{p}^{0,190}{\left({n}_{t}+1\right)}^{0,088}{\left({n}_{l}+1\right)}^{0,012}{\left({n}_{d}+1\right)}^{0,082}}}$

(29)

$B=1,346{\frac {{L}^{0,179}}{{W}^{0,394}}}$

(30)

donde ${\textstyle h}$ es el calado de aproximación, ${\textstyle A}$ y ${\textstyle B}$ los parámetros geométricos característicos de cada reja. Estos parámetros se obtienen a partir del área mínima que engloba los hueco ${\textstyle \,{A}_{g}}$, el porcentaje entre el área de huecos ${\textstyle \,{A}_{H}}$ y el área mínima que engloba los huecos ${\textstyle \,p=}$$100\left({A}_{H}/{A}_{g}\right)$ , el número de barras transversales ${\textstyle {n}_{t}}$, el número de barras longitudinales ${\textstyle {n}_{l}}$, el número de barras diagonales ${\textstyle {n}_{d}}$, la longitud de la reja ${\textstyle L}$ y ancho de la reja ${\textstyle W}$.

Es importante enfatizar que esta ecuación se obtuvo experimentalmente para el caudal que circula por un ancho de calzada de 3 m, por lo que se tuvieron las precauciones necesarias al momento de implementarla en el modelo numérico. Para anchos de calzada distintos, la Ecuación (28) puede generalizarse considerando una distribución uniforme de la velocidad [26], y así se hace en Iber [27].

De forma similar, para el caso de rejas transversales continuas, Gómez y Russo [28] determinaron una ley que relaciona la eficiencia hidráulica de captación de estas estructuras con el número de Froude ${\textstyle \,{F}_{r}}$, calado de aproximación y dos coeficientes característicos de la geometría ${\textstyle \,\alpha }$ y ${\textstyle \beta }$ , resultando en:

$E=\alpha \left[\,{F}_{r}{\left({\frac {h}{L}}\right)}^{0,812}\right]+\,\beta$

(31)

En este caso, la eficiencia hidráulica de captación de rejas continuas puede ser expresada en términos de una eficiencia unitaria, definida como el cociente entre el caudal interceptado y el caudal de aproximación, ambos por unidad de ancho.

Una vez determinada la eficacia con cualquiera de las dos fórmulas anteriores, el caudal interceptado ${\textstyle \,{Q}_{int}}$ se puede determinar como:

${Q}_{int}=EQ$

(32)

Esta metodología tiene la ventaja de poder extrapolar la formulación a otras rejas con características similares sin necesidad de ser ensayadas [26].

3.2 Interfaz de usuario

Como se ha indicado, para el funcionamiento del módulo de drenaje urbano se sigue la misma estrategia que para los módulos existentes: el cálculo se realiza con una malla de cálculo que tiene una serie de propiedades y condiciones, y dónde cada volumen finito corresponde a un elemento de la malla. La malla se suele generar a partir de una geometría previa, trasfiriendo posteriormente las propiedades asignadas en ella a la malla de cálculo. Es posible, sin embargo, asignar propiedades directamente a la malla, o modificar las que ya tiene.

A continuación, se describen los diferentes menús y ventanas desde las que se pueden activar los procesos descritos en las secciones anteriores y definir los correspondientes parámetros. El módulo se activa desde el menú Herramientas Iber >> Plug-ins >> Drenaje urbano. Esto permite al/la usuario/a utilizar e implementar únicamente variables relacionadas con los módulos que interaccionan con el drenaje urbano dual.

3.2.1 Formulaciones, opciones y parámetros generales del modelo

Las opciones de cálculo se especifican en Datos >> Datos del problema >> Drenaje urbano. Una vez activado el módulo de drenaje urbano, automáticamente aparece la pestaña ‘Drenaje Urbano’ en la ventana de Datos del Problema (Fig. 15), desde donde se deben definir el tipo de flujo (superficial, subterráneo o dual) y el tipo de método de intercambio entre el flujo superficial y el flujo subterráneo, ya sea según el método ‘UPC’ [26] o como un orificio/vertedero.

Fig. 15. Datos generales del módulo de drenaje urbano.

3.2.2 Red de drenaje

La red de drenaje puede estar compuesta por uniones, conductos (que pueden unirse mediante uniones) y zonas donde el flujo subterráneo puede fluir o almacenarse (superficie 2D). Los elementos que conectan el flujo superficial con el flujo en la red de drenaje subterránea son las rejas (o inlets en inglés), y son los que permiten el intercambio de caudales entre la parte superficial (p.ej. calles) y la parte subterránea (p.ej. la red de drenaje).

Esta red puede construirse de forma manual, empleando las herramientas disponibles en Iber para la creación de geometrías (líneas, superficies), o bien importarse desde archivos externos, pudiendo ser o no necesaria la implementación de ciertas condiciones a posteriori.

Definición manual de la red

Para la definición completa de la red es necesario especificar qué conductos la forman (líneas), las uniones entre conductos (puntos) y, en su caso, qué zonas existen para ser calculadas en dos dimensiones (superficies). A la definición de la red de drenaje se accede mediante el menú Datos >> Drenaje urbano >> Red de drenaje. Desde aquí es posible definir las características de las uniones, los conductos y las superficies de la red de drenaje.

Las Uniones se definen mediante sus características geométricas (Fig. 16a): la cota de fondo, su diámetro (o diámetro equivalente si no son circulares), el calado inicial y su cota superior (para láminas de agua superiores a la cota superior la unión funciona a presión). Además, es posible indicar si se desea permitir la continuidad del flujo entre conducciones de entrada y salida (opción Continuidad del flujo en la unión) o no. En caso de no existir continuidad, se supone que la velocidad del agua en la unión es nula, o lo que es lo mismo, que no hay continuidad entre la línea de energía de los conductos entrantes y salientes. En las uniones es posible asignar, también, un caudal de entrada.

En cuanto a las Conducciones, pueden ser de tipo circular o rectangular (Fig. 16b), y deben asignarse a líneas o elementos lineales. Los conductos pueden trabajar en presión, por lo que es necesario definir la celeridad de la onda (en metros por segundo).

Fig. 16. Definición e implementación de la red de drenaje: (a) uniones; (b) conducciones; (c) zonas de la subsuperficie.

Finalmente, está la posibilidad de que existan zonas en la red de drenaje que tengan una geometría más compleja y requieran un cálculo en 2D. Para ello se dispone de las Áreas de la subsuperficie (Fig. 16c). Aquí tan solo se tiene que indicar qué superficie corresponde a este tipo de estructura, y el programa automáticamente calculará la hidrodinámica en 2D teniendo en cuenta las entradas/salidas de flujo que haya en sus contornos (más adelante se indica cómo conectar dichas superficies con el resto de la red de drenaje).

Importación de la red

La importación, no solo de la geometría de la red, sino también de las características de la red (uniones, conductos, áreas no conectadas con la superficie) puede realizarse desde el menú Datos >> Drenaje urbano >> Importar SHP. Para ello es necesario disponer de un archivo ESRI Shapefile (*.shp) que contenga la definición geométrica de la red, ya sea en 2D o 3D.

Fig. 17. Opción ‘Crear condiciones’ para que Iber lea la base de datos del *.shp (archivo *.dbf asociado) y asigne las condiciones a la geometría.

Si el archivo *.shp va acompañado de una base de datos (*.dbf), se podrá emplear para asignar las condiciones de la red (tipo de unión, dimensiones del conducto, rugosidad, etc.) también de forma automática (opción ‘Crear condiciones’ activada, Fig. 17).

Este archivo debe tener una estructura concreta para que la importación de todas las condiciones sea conforme. En caso contrario, tan solo se importará la geometría. La estructura de la base de datos para la asignación automática de las condiciones de la red en Iber para elementos tipo Uniones se detalla en la Tabla 1:

Tabla 1. Estructura de la base de datos (*.dbf) que acompaña al archivo de geometría (*.shp) de las Uniones.

Col. 1	Col. 2	Col. 3	Col. 4	Col. 5	Col. 6	Col. 7
Elevación superior	Elevación del fondo	Calado	Diámetro	Continuidad flujo	Entrada/salida flujo	Tabla

Las Columnas 5 y 6 tienen que ser un número entero, pudiendo tomar 0, en caso de que no se desee habilitar dicha opción, o 1, en caso que se desee habilitar dicha opción. Si la Columna 6 toma valor 1, se espera que en la Columna 7 se defina la tabla de intercambio de entrada/salida de flujo en la unión. Esta tabla se tiene que definir de la siguiente manera:

• Tipo: texto
• Campo como sigue: #N# XX A1 A2 B1 B2 C1 C2 … n1 n2

El valor del campo de la Columna 7, separado mediante espacios, queda definido en primer lugar con “#N#”, seguido del número de elementos que tiene la tabla (“XX”) y luego los valores de la tabla agrupados por parejas, es decir, tiempo-caudal. Así, por ejemplo, un hidrograma de entrada de forma triangular simétrico con un caudal punta de 100 L/s a la 1h, deberá introducirse en la Columna 7 como: #N# 6 0 0 3600 0.1 7200 0.

En este sentido, la estructura de la base de datos para la asignación automática de las condiciones de la red en Iber para elementos tipo Conducciones es la que se detalla en la Tabla 2:

Tabla 2. Estructura de la base de datos (*.dbf) que acompaña al archivo de geometría (*.shp) de las Conducciones.

Col. 1	Col. 2	Col. 3	Col. 4	Col. 5	Col. 6	Col. 7
Núm. conducción	Tipo	Diámetro	Ancho	Alto	Manning	Celeridad

La Columna 2 define el tipo de conducción, por lo que debe ser un campo tipo Texto cuyo valor debe ser “Rectangular” o “Circular”. En caso que la Columna 2 sea tipo “Circular”, la Columna 3 debe tomar el valor del diámetro del conducto, mientras que las Columnas 4 y 5 pueden tomar cualquier valor, pero deben existir. En cambio, en caso que la Columna 2 sea tipo “Rectangular”, la Columna 3 podrá tomar cualquier valor, pero deberá existir, mientras que las Columnas 4 y 5 deben tomar el valor del ancho y el alto del conducto, respectivamente.

Una vez importada la geometría del modelo, es posible modificar los parámetros mediante los menús comentados para la definición manual de la misma.

3.2.3 Sumideros

Definición manual de los sumideros

La definición/modificación de los sumideros se realiza desde el menú Datos >> Drenaje urbano >> Rejas. La Fig. 18a muestra la ventana para la implementación de rejas, que pueden ser tipo longitudinal, transversal o como una reja circular (similar a un pozo de registro). La Fig. 18b muestra, a modo de ejemplo, cómo se representaría la reja (resaltada en amarillo) y su conexión con la red de drenaje.

Para cada reja o sumidero lo primero que hay que hacer es crearla con el botón correspondiente y luego introducir sus características. Los parámetros que hay que introducir son:

• Posición en superficie: Coordenadas
• Tipo: puede ser “To network” (a la red) o “Sink” (sumidero sin que el agua se incorpore a ningún colector)
• Posición en la red: Coordenadas del punto de conexión a la red
• Tipo: Transversal o longitudinal
• A: Parámetro A del método UPC
• B: Parámetro B del método UPC
• Longitud
• Anchura
• Área efectiva (no obstruida)
• Coeficiente de desagüe de vertedero para la metodología orificio/vertedero
• Coeficiente de desagüe de orificio para la metodología orificio/vertedero

En el caso de que exista flujo entrando a la reja, además del que intercepta de la superficie, es posible incorporarlo a través de un hidrograma.

Fig. 18. Definición e implementación de rejas: (a) ventana para la introducción de datos; (b) representación de la reja y su conexión con la red de drenaje.

Importación de sumideros

También, desde el menú Datos >> Drenaje urbano >> Importar archivo de rejas, se posibilita la importación de la localización y propiedades de las rejas. Se trata de un archivo de texto que debe contener el identificador de la reja, su ubicación, el tipo (longitudinal, transversal o pozo de registro), sus características geométricas y propiedades hidráulicas. También es posible indicar un hidrograma de entrada/salida en ese elemento.

La Fig. 19 presenta la estructura del archivo requerido. Cada fila corresponde a una reja con los siguientes parámetros:

• Número de reja
• 1 si es visible (se representa en ventana), 0 en caso contrario
• Coordenadas (x,y,z) en superficie
• 1 si es un sumidero, 2 si hay conexión a la red
• Coordenadas (x,y,z) de la conexión al colector
• 1 si es longitudinal, 2 si es transversal 3 si es pozo de registro
• A
• B
• Longitud
• Anchura
• Porcentaje de área efectiva
• C_d Vertedero
• C_d Orificio
• 0 si no hay hidrograma asociado, 1 en caso que haya
• Nombre del sumidero (sin espacios)
• Número de valores del hidrograma (1 en caso de que no haya hidrograma)
• Tantas filas como valores del hidrograma, con parejas de datos tiempo-caudal

Fig. 19. Estructura del archivo de importación automática de rejas.

El formato del archivo de importación corresponde con el formato del archivo “Iber_Sink&Sources.dat”, que se crea dentro de la carpeta del proyecto al ejecutar un cálculo con sumideros introducidos manualmente. De esta manera, leyendo este archivo de un proyecto ya creado, se pueden importar los sumideros a un proyecto nuevo.

3.2.4 Condiciones de contorno y conexiones

La red de drenaje puede presentar diferentes tipos de condiciones de contorno, así como conexiones directas con la superficie (además de las Rejas). Para ello se han generado nuevas ventanas donde se pueden implementar estas condiciones.

Fig. 20. Ventana para la introducción de datos de las condiciones de contorno de entrada de la red: (a) régimen supercrítico; (b) crítico/subcrítico.

Por un lado, se encuentran las condiciones de contorno de la propia red. Su implementación se realiza desde el menú Datos >> Hidrodinámica >> Condiciones de contorno de la red. Se pueden indicar las condiciones de contorno de entrada, que pueden ser en régimen supercrítico (Fig. 20a) o en régimen crítico/subcrítico (Fig. 20b). En cada caso se muestran el tipo de condiciones que se pueden implementar, ya sea una tabla Tiempo-Caudal-Elevación/Calado o una tabla Tiempo-Caudal (hidrograma). Para las condiciones de salida, se sigue una filosofía análoga, cambiando tan solo el tipo de salida según el régimen hidráulico. Estas condiciones deben imponerse en un contorno de la red (puntos extremos de los colectores).

Por otro lado, se encuentran las condiciones de contorno para elementos de tipo superficie, que son análogos a los ya existentes en Iber para la modelización de flujo en superficie. La ventana es accesible desde el menú Datos >> Hidrodinámica >> Condiciones de contorno, y ha sido adaptada para contener opciones de conexión entre elementos 1D (conducciones de la red, Fig. 21a) y elementos 2D (Zonas de la subsuperfice o zonas de la superficie, Fig. 21b).

La existencia de condiciones de contorno, tanto en superficie como en la red de drenaje, requiere conectar ambos tipos de flujo. Por ejemplo, es posible que el flujo que sale por una condición 2D sea el flujo que entra en una condición 1D, y viceversa. Por ese motivo existen las opciones “Conexión 1D” y “Conexión 2D” en las condiciones de contorno 2D y 1D, respectivamente.

Fig. 21. Condiciones de contorno de la red: (a) opciones para la conexión entre condiciones salida/entrada 1D; (b) condiciones entrada/salida 2D.

En las condiciones de contorno de la red de drenaje (1D), se permite la conexión con entidades 2D. En ese caso, se solicita el punto de la red que va a ser considerado como una entrada/salida en una entidad 2D del modelo. Del mismo modo, en una condición de contorno 2D se permite la conexión con la red de drenaje 1D, es necesario indicar el punto de la red en la que va a entrar/salir flujo. Para facilitar su selección, se permite la selección del punto de la red a través de la interfaz.

4 Módulo de hidrología distribuido

4.1 Fundamentos teóricos

El módulo de procesos hidrológicos de Iber permite calcular la transformación lluvia-escorrentía en toda una cuenca, haciendo posible su utilización como modelo hidrológico distribuido en cuencas de pequeño y mediano tamaño. Es completamente bidimensional y distribuido en el espacio, y trabaja sobre la misma malla de cálculo que el módulo hidrodinámico.

Este módulo se basa en las ecuaciones hidrodinámicas de aguas poco profundas bidimensionales, las cuales se utilizan para resolver el movimiento de la escorrentía superficial sobre el terreno. Se le han añadido una serie de nuevos procesos como son, el aporte de agua de lluvia, la infiltración de agua en el terreno, el flujo subsuperficial, la evapotranspiración y la percolación.

4.1.1 Procesos considerados

La base del módulo hidrológico es el módulo hidrodinámico, al que se le han añadido los procesos de aporte de agua por precipitación y la infiltración del agua en el terreno. De esta forma, el exceso de precipitación se mueve por el terreno y por los cauces de acuerdo con las ecuaciones de aguas someras bidimensionales. El agua infiltrada pasa a formar parte del subsuelo, donde se realiza un balance de masa de agua distribuido en el espacio a nivel de elemento de la malla de cálculo. En dicho balance de agua en el terreno se consideran los procesos de infiltración, evapotranspiración, percolación y flujo en el subsuelo. Los tres primeros implican un movimiento de agua en la columna vertical, mientras que el flujo en el subsuelo se considera bidimensional y paralelo al terreno. Todos estos procesos, representados de forma esquemática en la Fig. 22, pueden ser activados o desactivados por el usuario.

Fig. 22. Esquema de los procesos considerados en el módulo hidrológico.

4.1.2 Escorrentía superficial

El cálculo de la escorrentía superficial, y su propagación a lo largo de la cuenca hidrográfica (incluyendo laderas y cauces de agua), se realiza mediante la resolución de las ecuaciones de aguas poco profundas bidimensionales, incluyendo términos que representan las componentes de precipitación e infiltración (Ecuación (1)).

Este sistema de ecuaciones es el mismo utilizado en Iber para calcular la hidrodinámica en un curso fluvial. La única diferencia es que se han incluido 2 nuevos esquemas numéricos para su resolución, que suelen dar resultados numéricamente más estables en aplicaciones hidrológicas. Estos modelos son los siguientes:

• DHD
• DHD_basin

La descripción detallada del modelo DHD se presenta en el artículo Cea y Bladé [29]. En cuanto al modelo DHD_basin, es una variante del modelo DHD que se diferencia únicamente en el tratamiento de los frentes que separan las zonas secas de las mojadas (frente seco-mojado). En el esquema DHD, la cota de la lámina de agua en los frentes seco-mojado se calcula utilizando la ecuación (15) de Cea y Bladé [29]. En el esquema DHD_basin, la altura de la lámina en los frentes seco-mojado se calcula mediante una interpolación lineal de la cota de agua en el elemento mojado y la cota del terreno en el elemento seco.

Se recomienda la utilización de uno de estos dos esquemas (preferiblemente el DHD_basin) en aplicaciones que impliquen transformación lluvia-escorrentía en la totalidad de una cuenca hidrológica. Sin embargo, estos esquemas producen resultados menos precisos que los esquemas de Roe de orden 1 y orden 2 para el cálculo de flujo en ríos, canales o estuarios, por lo que no deben utilizarse en este tipo de aplicaciones.

4.1.3 Precipitación

La precipitación puede definirse a partir de series temporales de intensidad de precipitación observadas en pluviómetros de coordenadas conocidas (hietogramas), o a través de ficheros raster georreferenciados. Estos últimos se definen como campos espaciales de intensidad o profundidad de precipitación en diferentes instantes de tiempo. En ambos casos, la precipitación puede ser variable en tiempo y espacio.

En el caso de imponer la lluvia a partir de hietogramas, el usuario puede definir tantos pluviómetros como desee. Dichos pluviómetros pueden asignarse manualmente a diferentes zonas del modelo, con el fin de definir manualmente la distribución espacial de la lluvia. En el caso de que no se asignen manualmente, el programa realizará una interpolación espacial de la lluvia utilizando el método del vecino más cercano (natural neighbour), utilizando para ello las coordenadas espaciales de cada pluviómetro.

En el caso de definir la lluvia a partir de ficheros raster, el usuario debe suministrar una serie de ficheros georreferenciados (en el mismo sistema de referencia que el resto del modelo) en formato ASCII Grid de ArcInfo. Cada raster define la intensidad de precipitación o la profundidad de precipitación durante un intervalo de tiempo específico. Todos los ficheros raster deben tener el mismo número de columnas y filas, el mismo tamaño de pixel y el mismo origen. Los ficheros raster pueden tener cualquier nombre, pero todos ellos deben copiarse en una subcarpeta llamada Rain_rasters, localizada dentro de la carpeta del proyecto (proyecto.gid/Rain_rasters). La carpeta Rain_rasters debe crearla el usuario.

4.1.4 Infiltración

Se incluyen los cuatro modelos de infiltración siguientes:

• Lineal
• Green-Ampt
• Número de Curva (SCS)
• Horton

Estas formulaciones son las más extendidas para el cálculo de la infiltración, brindando suficiente versatilidad al usuario a la hora de definir y parametrizar las pérdidas por infiltración. En todo caso, la infiltración de agua en un terreno real es un proceso complejo que puede, en ocasiones, no estar correctamente representado por ninguna de estas formulaciones. Asimismo, la determinación de los parámetros de estas formulaciones suele requerir su calibración mediante datos de campo.

El modelo lineal es el más sencillo, y permite introducir unas pérdidas iniciales seguidas de una infiltración potencial constante en el tiempo. Puede ser adecuado para eventos cortos en los cuales puede considerarse la infiltración más o menos constante en el tiempo. También puede utilizarse en casos en los que la capacidad de infiltración del suelo sea muy elevada inicialmente (superior a la precipitación), hasta que llega un momento en el que se satura el suelo y la infiltración disminuye drásticamente a un valor aproximadamente constante en el tiempo.

El modelo de Green-Ampt utiliza una formulación basada en parámetros físicos que pueden relacionarse con ciertas propiedades del terreno, como son la porosidad, la permeabilidad saturada, la succión o el contenido de humedad del suelo. Genera una infiltración potencial elevada inicialmente (dependiente de la succión y de la humedad inicial del suelo), que va decayendo de forma aproximadamente exponencial en el tiempo, hasta llegar a un valor mínimo de infiltración igual a la permeabilidad saturada del terreno.

El modelo del Número de Curva del SCS es una formulación ampliamente conocida y utilizada en múltiples modelos hidrológicos. A unas pérdidas iniciales le sigue una capacidad de infiltración variable, que disminuye de forma aproximadamente exponencial en el tiempo. Su principal ventaja es que depende de un único parámetro, cuyo valor puede ser relacionado con la pendiente del terreno, el tipo de suelo y los usos del suelo. Existe, además, numerosa bibliografía para estimar su valor.

El modelo de Horton es similar al de Green-Ampt en cuanto a la capacidad de infiltración generada, pero sus parámetros no están directamente relacionados con las características del terreno.

Modelo de infiltración lineal

En el modelo lineal, la tasa de infiltración potencial ( ${\textstyle {f}_{pot}}$) se asume constante en el tiempo. El usuario introduce su valor en mm/h como un parámetro del modelo, que puede ser variable en espacio. La infiltración real ( ${\textstyle f}$) en cada instante de tiempo es igual a la infiltración potencial siempre y cuando la intensidad de lluvia sea superior a ${\textstyle {f}_{pot}}$, o si el calado sobre la superficie es suficiente para satisfacer su valor.

Además, el modelo lineal considera unas pérdidas iniciales ( ${\textstyle {I}_{a}}$) definidas por el usuario en mm. Las pérdidas iniciales se satisfacen con la lluvia inicial o con el calado de agua sobre la superficie del terreno.

Modelo de infiltración de Green-Ampt

En el modelo de Green-Ampt [30] se asume que existe un frente de saturación que separa una capa de suelo saturado, situada inmediatamente bajo la superficie del terreno, de una capa inferior de suelo no saturado. A medida que el agua se infiltra en el terreno, el frente saturado desciende y el espesor de la capa saturada ( ${\textstyle b}$) aumenta.

En la implementación del modelo de Green-Ampt en Iber, la tasa de infiltración potencial ( ${\textstyle {f}_{pot}}$) y la infiltración real ( ${\textstyle f}$) se calculan como:

${f}_{pot}={k}_{s}\left(1+{\frac {\left(h+\psi \right)}{b}}\right)\quad \quad f=$$\mathrm {min} \,\left({f}_{pot},\,{\frac {h}{\Delta t}}\right)$

(33)

donde ${\textstyle {k}_{s}}$ es la permeabilidad saturada del suelo en la dirección vertical, ${\textstyle h}$ es el calado de agua sobre el terreno, ${\textstyle \psi }$ es la succión en la capa de suelo no saturada y ${\textstyle b}$ es el espesor de la capa de suelo saturada. La succión ( ${\textstyle \psi }$ ) y la permeabilidad saturada ( ${\textstyle {k}_{s}}$) son parámetros del modelo que debe introducir el usuario, mientras que el espesor de la capa de suelo saturada se actualiza a medida que se infiltra más agua en el terreno mediante las siguientes expresiones:

$b={b}_{0}+{\frac {F}{\phi }}\quad \quad \,F=\int _{0}^{t}f\,dt$

(34)

donde ${\textstyle {b}_{0}}$ es el espesor inicial de la capa de suelo saturada, que depende del contenido de humedad inicial del suelo al comienzo de la simulación, y ${\textstyle \phi }$ es la porosidad del suelo.

Dado que Iber es un modelo 2D en el que no se resuelve el flujo en dirección vertical, el espesor de la capa de suelo saturada ( ${\textstyle b}$) se relaciona con el contenido de humedad promediado verticalmente en la capa de suelo ( ${\textstyle \theta }$ ), mediante la siguiente relación:

${\frac {b}{{d}_{s}}}={\frac {\theta }{\phi }}\quad \quad {b}_{0}={\theta }_{0}{\frac {{d}_{s}}{\phi }}$

(35)

donde ${\textstyle {d}_{s}}$ es el espesor total de la capa de suelo, ${\textstyle \theta }$ es el contenido de humedad promediado verticalmente en la capa de suelo. Es decir, el valor de ${\textstyle b}$ se obtiene agrupando todo el contenido de agua de la capa de suelo en una capa horizontal situada justo bajo el terreno (Fig. 23). Los valores del espesor total de la capa de suelo ( ${\textstyle {d}_{s}}$), su porosidad ( ${\textstyle \phi }$ ) y la humedad inicial ( ${\textstyle {\theta }_{0}}$), son parámetros del modelo que debe introducir el usuario.

Fig. 23. Esquema de la distribución de humedad en la dirección vertical en una capa de suelo.

Cuando el espesor de la capa saturada alcanza el espesor total de la capa de suelo ( ${\textstyle b=}$${d}_{s}$), ya no es posible infiltrar más agua. Esta situación se corresponde a la de un suelo totalmente saturado ( ${\textstyle \theta =}$$\phi$ ).

El modelo de Green-Ampt permite, asimismo y de forma opcional, definir unas pérdidas iniciales ( ${\textstyle {I}_{a}}$) definidas por el usuario en mm. Las pérdidas iniciales se satisfacen con la lluvia inicial o con el calado de agua sobre la superficie del terreno.

En la implementación del modelo de Green-Ampt se consideran, por tanto, un máximo de seis parámetros: ${\textstyle {k}_{s}}$, ${\textstyle \psi }$ , ${\textstyle \phi }$ , ${\textstyle {d}_{s}}$, ${\textstyle {\theta }_{0}}$, ${\textstyle {I}_{a}}$.

La principal ventaja del modelo de Green-Ampt es que se puede recuperar capacidad de infiltración realizando un balance de masa de agua en el suelo, incorporando los procesos de percolación hacia capas más profundas, de evapotranspiración y de flujo subsuperficial. En efecto, teniendo en cuenta la relación entre la humedad del suelo y el espesor de la capa saturada, la infiltración potencial puede escribirse como:

${f}_{pot}={k}_{s}\left(1+{\frac {\left(h+\psi \right)}{\theta }}\,{\frac {\phi }{{d}_{s}}}\right)$

(36)

De esta forma, si la percolación y la evapotranspiración superan a la infiltración, la humedad del suelo disminuirá y, por tanto, la capacidad de infiltración aumentará. En el caso de que no se consideren los procesos de percolación o evapotranspiración, el contenido de humedad siempre aumentará y la capacidad de infiltración disminuirá de forma continua.

El modelo puede simplificarse reduciendo el número de parámetros. Si el usuario define las pérdidas iniciales y la humedad inicial igual a cero ( ${\textstyle {I}_{a}=}$$0$ y ${\textstyle {\theta }_{0}=0}$), y un espesor de suelo suficientemente elevado ( ${\textstyle {d}_{s}\rightarrow \infty }$ ), los parámetros del modelo se reducen a los tres utilizados en la implementación clásica del modelo de Green-Ampt: ${\textstyle {k}_{s}}$, ${\textstyle \psi }$ , ${\textstyle \phi }$ , y la infiltración se calculará como:

${f}_{pot}={k}_{s}\left(1+{\frac {\left(h+\psi \right)\,\phi }{F}}\right)$

(37)

Modelo de infiltración del SCS

En el método del número de curva del SCS [31-33], se define la capacidad máxima de retención del suelo ( ${\textstyle S}$) como:

$S={\frac {25400}{CN}}-254$

(38)

donde ${\textstyle CN}$ es el Número de Curva, a definir por el usuario en función del tipo de suelo, usos del suelo y pendientes del terreno. La abstracción inicial ( ${\textstyle {I}_{a}}$) se define como un porcentaje de la capacidad máxima de retención, como:

${I}_{a}=\alpha \,S$

(39)

En la implementación clásica del método del número de curva, la abstracción inicial es un 20% de la capacidad máxima de retención ( ${\textstyle {I}_{a}=}$$0,2\,S$). Sin embargo, diversos estudios han encontrado que en muchos casos las pérdidas iniciales son inferiores a dicho valor, por lo que en la implementación en Iber se introduce el valor de ${\textstyle \alpha }$ como un parámetro a definir por el usuario.

A partir de los parámetros anteriores, la precipitación neta acumulada en cada instante se calcula como:

${P}_{n}={\frac {{\left(P-\alpha \,S\right)}^{2}}{P+\left(1-\alpha \right)\,S}}\quad P>\alpha \,S$

(40)

donde ${\textstyle P}$ es la profundidad de precipitación bruta acumulada, y ${\textstyle {P}_{n}}$ es la profundidad de precipitación neta acumulada. La ecuación anterior es válida una vez que se han satisfecho las pérdidas iniciales ( ${\textstyle P>{I}_{a}}$). Para valores de precipitación bruta inferiores a ${\textstyle {I}_{a}}$, la precipitación neta es igual a cero.

La infiltración acumulada se calcula en cada paso de tiempo como la diferencia entre la precipitación neta y la bruta:

$F=P-{P}_{n}$

(41)

A partir de la infiltración acumulada en dos pasos de tiempo consecutivos, la tasa de infiltración se calcula como:

$f={\frac {{F}_{(t+\Delta t)}-{F}_{\left(t\right)}}{\Delta t}}$

(42)

Una diferencia con la implementación clásica del método del número de curva en la implementación de las pérdidas iniciales en Iber, es que pueden satisfacerse con la precipitación inicial o con la escorrentía superficial que haya sobre el terreno antes de que comience la precipitación. Es decir, puede haber infiltración antes de que comience a llover. Una vez que se han satisfecho las pérdidas iniciales, las pérdidas por infiltración se satisfacen únicamente con la precipitación (a diferencia de los otros métodos de infiltración, que pueden continuar infiltrando agua incluso aunque no llueva).

Modelo de infiltración del SCS continuo

Esta es una variante del modelo SCS que permite recuperar capacidad de infiltración durante períodos de tiempo en los que no llueve. La formulación implementada en Iber está basada en la formulación propuesta en Kannan et al. [34] para el cálculo de infiltración a nivel diario. La principal diferencia es que en Iber se implementa con un paso de tiempo igual al paso de tiempo de cálculo, mientras que en Kannan et al. se plantea únicamente para valores diarios de precipitación e infiltración.

El método se basa en recalcular, en cada paso de tiempo, la capacidad de retención del suelo ( ${\textstyle S}$) a partir de la siguiente ecuación:

${S}^{n}={S}^{n-1}+\left({e}^{n-1}+{p}^{n-1}-{f}^{n-1}\right)\cdot \Delta t$

(43)

donde ${\textstyle S}$ es la capacidad de retención del suelo en un instante determinado, ${\textstyle f}$ es la infiltración real, ${\textstyle e}$ es la evapotranspiración y ${\textstyle p}$ es la percolación. Los superíndices ${\textstyle n}$ y ${\textstyle n-}$$1$ hacen referencia a los instantes de tiempo en los que se calcula cada término. Por lo tanto, la capacidad de retención del suelo disminuye a medida que infiltra el agua, y aumenta a medida que evapotranspira o que percola. Si el usuario no activa la percolación ni la evapotranspiración, el modelo no recuperará capacidad de infiltración.

A partir del valor de la capacidad de retención del suelo, se calcula la infiltración en cada paso de tiempo como:

$f={\frac {h-{h}^{\ast }}{\Delta t}}\,\,$

(44)

donde ${\textstyle h}$ es el calado sobre el terreno antes de producirse la infiltración, y ${\textstyle {h}^{\ast }}$ es el calado restante una vez satisfecha la capacidad de infiltración en este paso de tiempo, el cual se evalúa como:

${h}^{\ast }=\left\{{\begin{matrix}{\frac {{\left(h-\alpha \,{S}_{\Delta t}\right)}^{2}}{h+\left(1-\alpha \right)\,{S}_{\Delta t}}}&&h>\alpha \,{S}_{\Delta t}\\&&\\0&&h<\alpha \,{S}_{\Delta t}\end{matrix}}\right.$

(45)

con:

${S}_{\Delta t}=S\,{\frac {\Delta t}{86400}}$

(46)

siendo ${\textstyle \Delta t}$ el paso de tiempo de cálculo expresado en segundos. El valor de la capacidad de retención ${\textstyle S}$ se actualiza a cada paso de tiempo y se limita entre 0 y ${\textstyle {S}_{max}}$, siendo ${\textstyle {S}_{max}}$ la máxima capacidad de retención del terreno, calculada a partir del número de curva para condiciones antecedentes de humedad secas ( ${\textstyle C{N}_{dry}}$), como:

${\begin{matrix}{S}_{max}={\frac {25400}{C{N}_{dry}}}-254\\C{N}_{dry}={\frac {4,2\,CN}{10-0,058\,CN}}\end{matrix}}$

(47)

Modelo de infiltración de Horton

En el modelo de Horton la tasa de infiltración potencial se calcula como:

${\begin{matrix}{f}_{pot}={f}_{r}+\left({f}_{0}-{f}_{r}\right)\,{e}^{-kt}\\\,\quad \quad f=\mathrm {min} \,\left({f}_{pot},\,{\frac {h}{\Delta t}}\right)\end{matrix}}$

(48)

donde ${\textstyle t}$ es el tiempo de comienzo de la precipitación, ${\textstyle {f}_{o}}$ es la tasa de infiltración inicial (máxima), ${\textstyle {f}_{r}}$ es la tasa de infiltración una vez que el suelo está completamente saturado (mínima) y ${\textstyle k}$ es un coeficiente de decaimiento que define la variación temporal de la tasa de infiltración.

La ecuación anterior no se puede aplicar directamente en un modelo como Iber en el que la precipitación puede ser intermitente, ya que la capacidad de infiltración disminuiría de forma constante con el tiempo incluso durante períodos en los cuales no hubiese precipitación ni infiltración. Por ello es necesario implementar una ecuación en la cual la capacidad de infiltración disminuya a medida que el agua se infiltra en el terreno.

En Iber se implementa de la siguiente forma. En cada paso de tiempo se calcula la infiltración potencial como la integral de ${\textstyle {f}_{pot}}$ entre ${\textstyle t}$ y ${\textstyle t+}$$\Delta t$, y se actualiza el valor de ${\textstyle {f}_{0}}$ según la siguiente ecuación:

$\Delta {F}_{pot}={f}_{r}\,\Delta t+{\frac {\left({f}_{0}-{f}_{r}\right)}{k}}\left(1-\right.$$\left.{e}^{-k\,\Delta t}\right)$

(49)

Si ${\textstyle \Delta {F}_{pot}<h}$:

${\begin{matrix}\Delta F=\,\Delta {F}_{pot}\\{f}_{o}={f}_{r}+\left({f}_{0}-{f}_{r}\right)\,{e}^{-k\,\Delta t}\end{matrix}}$

(50)

Si ${\textstyle \Delta {F}_{pot}>h}$:

${\begin{matrix}\Delta F=h\\{f}_{o}={f}_{r}+\left({f}_{0}-{f}_{r}\right)\,{e}^{-k\,{\frac {\Delta F}{f}}}\end{matrix}}$

(51)

Este proceso se repite en cada paso de tiempo de cálculo. En el caso de que siempre haya agua suficiente en la superficie para satisfacer la capacidad de infiltración, esta implementación se corresponde con una discretización temporal de la ecuación clásica de Horton, con un decaimiento exponencial continuo en el tiempo de la capacidad de infiltración. Sin embargo, en el caso de que no haya agua disponible para infiltrarse porque no llueve y el terreno está seco ( ${\textstyle h=}$$0$), entonces ${\textstyle \Delta F=0}$ y ${\textstyle {f}_{0}=}$$cte$. Es decir, en esta situación no se infiltra agua y la tasa de infiltración potencial deja de disminuir porque no está aumentando la humedad del suelo.

Resumen de los modelos de infiltración implementados

En la Tabla 3 se resumen las principales características de los modelos de infiltración implementados en Iber.

Tabla 3. Formulaciones implementadas en Iber expresadas como capacidad potencial de infiltración ( ${\textstyle {\mathit {\boldsymbol {f}}}_{}^{\mathit {\boldsymbol {pot}}}}$).

Modelo de infiltración	Parámetros	Notas
Lineal	2	Incluye abstracción inicial. Capacidad de infiltración constante en el tiempo.
Green-Ampt	Entre 3 y 6	Incluye abstracción inicial. Si se activan la percolación y/o evapotranspiración, se considera la recuperación de capacidad de infiltración en el tiempo. En caso contrario la infiltración potencial disminuye en el tiempo.
SCS	2	Incluye abstracción inicial. Se incluye un parámetro para calibrar las pérdidas iniciales. La infiltración potencial disminuye en el tiempo.
SCS continuo	3	Incluye abstracción inicial. Se incluye un parámetro para calibrar las pérdidas iniciales. Si se activan la percolación y/o evapotranspiración, se considera la recuperación de capacidad de infiltración en el tiempo. En caso contrario la infiltración potencial disminuye en el tiempo.
Horton	3	No incluye pérdidas iniciales. La infiltración potencial disminuye en el tiempo.

4.1.5 Evapotranspiración

La evapotranspiración potencial se calcula mediante la formulación propuesta por Blaney y Criddle, posteriormente modificada por Doorenbos y Pruitt [35]. Se ha optado por esta formulación por su sencillez y porque únicamente necesita un parámetro. La evapotranspiración potencial se calcula mediante la siguiente ecuación:

$ETP=-2+b\,\left[\xi \,\left(0,46\,{T}_{a}+8,13\right)\right]$

(52)

donde ${\textstyle ETP}$ es la evapotranspiración potencial en mm/día, ${\textstyle {T}_{a}}$ es la temperatura del aire en grados centígrados, ${\textstyle \xi }$ es el porcentaje de horas diurnas totales en el día considerado sobre el total de horas diurnas en un año (365 x 12) y ${\textstyle b}$ es un parámetro del modelo que debe calibrarse. Como término de referencia para el valor de ${\textstyle \xi }$ , en el equinocio de primavera u otoño (12 horas de día), se obtiene un valor de ${\textstyle \xi =}$$0,27$. En la versión actual de Iber valor de ${\textstyle \xi }$ está prefijado a 0,27 (12 horas de día de sol), y no puede ser modificado por el usuario.

En cuanto a los valores habituales de ${\textstyle b}$, suelen estar comprendidos entre 0,7 y 2. El usuario puede modificar este valor a su conveniencia, con el fin de calibrar el modelo.

Evapotranspiración real con el modelo de infiltración de Green-Ampt

En el caso de que se seleccione el modelo de Green-Ampt la evapotranspiración real, en mm/día, se calcula a partir de la potencial como:

$e=ETP{\frac {\theta }{\phi }}$

(53)

siendo ${\textstyle \phi }$ la porosidad del suelo y ${\textstyle \theta }$ la humedad del suelo, que debe estar comprendida entre 0 y ${\textstyle \phi }$ . Es decir, si el suelo está completamente saturado, la evapotranspiración real es igual a la potencial, y si el suelo está totalmente seco, la evapotranspiración real es igual a cero. Entre estos dos estados extremos, se establece una relación lineal.

Además, se incluye una formulación para considerar el estrés hídrico en el cálculo de la evapotranspiración real. Cuando la humedad del subsuelo es menor a un determinado valor umbral, las plantas tienen una gran dificultad para extraer el agua disponible y se dice que están bajo estrés hídrico. En estas condiciones se produce una disminución de la tasa de evapotranspiración. Si se activa la opción de ‘Estrés Hídrico’, se utiliza el método conocido como Coeficiente Único del Cultivo (single crop coefficient) para calcular la evapotranspiración como:

$e={K}_{s}\,{K}_{c}\,ETP$

(54)

dónde ${\textstyle {K}_{s}}$ es el coeficiente de estrés hídrico, un factor adimensional comprendido entre 0 y 1 que es función del contenido en agua del suelo, y ${\textstyle {K}_{c}}$ es un coeficiente que depende del tipo de cultivo y del estado del cultivo (época del año), comprendido entre 0,3 y 1,2.

El valor de ${\textstyle {K}_{c}}$ lo debe introducir el usuario en función del tipo de cultivo, y puede ser variable en el tiempo. El coeficiente ${\textstyle {K}_{s}}$ se evalúa en función del punto de marchitez (wilting point) y de la capacidad de campo (field capacity), como:

${K}_{s}=\left\{{\begin{matrix}0&si&\theta \leq {\theta }_{wp}\\{\frac {\theta -{\theta }_{wp}}{\left({\theta }_{fc}-{\theta }_{wp}\right)\left(1-p\right)}}&si&\theta <{\theta }_{t}\\1&si&\theta \geq {\theta }_{t}\end{matrix}}\right.\,\,$

(55)

siendo ${\textstyle {\theta }_{t}={\theta }_{fc}\left(1-p\right)+{\theta }_{wp}\,p}$, la humedad del suelo por debajo de la cual comienza a haber estrés hídrico, ${\textstyle {\theta }_{wp}}$ el punto de marchitez, ${\textstyle {\theta }_{fc}}$ la capacidad de campo, y ${\textstyle p}$ el factor de agotamiento (depletion factor), que debe estar comprendido entre 0 y 1.

Evapotranspiración real con el modelo de infiltración de SCS continuo

En el caso de que se active el modelo de infiltración del SCS continuo, la evapotranspiración real se calcula como:

$e=ETP\,\cdot \,exp\left(-B{\frac {S}{{S}_{max}}}\right)\,$

(56)

donde ${\textstyle B}$ es un coeficiente de agotamiento a definir por el usuario, ${\textstyle S}$ es la capacidad de retención del terreno en el instante considerado, y ${\textstyle {S}_{max}}$ es la máxima capacidad de retención del terreno, calculada a partir del número de curva ( ${\textstyle CN}$) para condiciones antecedentes de humedad secas ( ${\textstyle C{N}_{dry}}$), como:

${\begin{matrix}{S}_{max}={\frac {25400}{C{N}_{dry}}}-254\\C{N}_{dry}={\frac {4,2\,CN}{10-0,058\,CN}}\end{matrix}}$

(57)

donde ${\textstyle CN}$ es el número de curva introducido por el usuario como parámetro del modelo SCS continuo. El coeficiente de agotamiento ${\textstyle B}$ suele estar comprendido en el rango ${\textstyle \left[0,2\right]}$ .

Evapotranspiración real con los modelos de infiltración Lineal, SCS y Horton

Cuando se activa un modelo de infiltración diferente a Green-Ampt o a SCS continuo, únicamente se puede calcular la evapotranspiración si al mismo tiempo se activa el flujo subsuperficial y se selecciona el modelo agregado. En caso de que no se active el flujo subsuperficial no tiene sentido calcular la evapotranspiración con los modelos de infiltración Lineal, SCS u Horton, ya que ninguno de ellos hace un seguimiento de la humedad del suelo ni puede recuperar capacidad de infiltración.

En caso de que se active el modelo agregado de flujo subsuperficial en combinación con los modelos de infiltración Lineal, SCS o Horton, la evapotranspiración real se calcula a partir de la potencial como:

$e=ETP{\frac {{V}_{s}}{{V}_{s,max}\,}}$

(58)

siendo ${\textstyle {V}_{s}}$ el volumen de agua en el subsuelo, y ${\textstyle {V}_{s,max}}$ el volumen de agua máximo en el subsuelo. Para realizar el cálculo de la evapotranspiración, el cálculo del volumen de agua en el subsuelo se realiza a nivel de depósito subterráneo, tal como se detalla en la descripción del modelo agregado de flujo subsuperficial.

4.1.6 Percolación

Al igual que la evapotranspiración, la percolación se calcula de forma diferente en función del modelo de infiltración seleccionado.

Percolación con el modelo de infiltración de Green-Ampt

Si se selecciona el modelo de infiltración de Green-Ampt, la percolación se calcula mediante la ecuación propuesta por Famiglietti y Wood [36]:

$p={k}_{p}{\left({\frac {\theta }{\phi }}\right)}^{\left(3+{\frac {2}{\lambda }}\right)}$

(59)

donde ${\textstyle p}$ es la tasa de percolación, ${\textstyle {k}_{p}}$ es la permeabilidad de las capas inferiores del suelo en las cuales se produce la percolación, ${\textstyle \phi }$ la porosidad del suleo, ${\textstyle \theta }$ la humedad del suelo (que debe estar comprendida entre 0 y ${\textstyle \phi }$ ) y ${\textstyle \lambda }$ es un parámetro del modelo (soil size distribution index). Los valores habituales de ${\textstyle \lambda }$ están comprendidos entre 0,05 y 0,5.

Percolación con el modelo de infiltración de SCS continuo

Si se selecciona el modelo de infiltración SCS continuo, la percolación se calcula mediante la siguiente ecuación:

$p={k}_{p}{\left(1-{\frac {S}{{S}_{max}}}\right)}^{\left(3+{\frac {2}{\lambda }}\right)}$

(60)

donde ${\textstyle S}$ es la capacidad de retención del terreno en el instante considerado, y ${\textstyle {S}_{max}}$ es la máxima capacidad de retención del terreno, calculada a partir del número de curva para condiciones antecedentes de humedad secas.

Percolación con los modelos de infiltración Lineal, SCS y Horton

Cuando se activa un modelo de infiltración diferente a Green-Ampt o a SCS continuo, únicamente se puede calcular la percolación si al mismo tiempo se activa el flujo subsuperficial y se selecciona el modelo agregado. En caso de que no se active el flujo subsuperficial no tiene sentido calcular la percolación con los modelos de infiltración Lineal, SCS o Horton, ya que ninguno de ellos hace un seguimiento de la humedad del suelo ni puede recuperar capacidad de infiltración.

En el caso de que se active el modelo agregado de flujo subsuperficial en combinación con los modelos de infiltración Lineal, SCS o Horton, la percolación se calcula como:

$p={k}_{p}{\left({\frac {{V}_{s}}{{V}_{s,max}}}\right)}^{\left(3+{\frac {2}{\lambda }}\right)}$

(61)

siendo ${\textstyle {V}_{s}}$ el volumen de agua en el subsuelo, y ${\textstyle {V}_{s,max}}$ el volumen de agua máximo en el subsuelo. Para realizar el cálculo de la percolación, el cálculo del volumen de agua en el subsuelo se realiza a nivel de depósito subterráneo, tal como se detalla en la descripción del modelo agregado de flujo subsuperficial.

4.1.7 Flujo subsuperficial

En el módulo hidrológico de Iber se incluyen dos modelos de flujo subsuperficial: un modelo agregado y un modelo distribuido. El modelo distribuido debe utilizarse conjuntamente con el modelo de infiltración de Green-Ampt, ya que hereda de dicho modelo varios de sus parámetros. El modelo agregado puede utilizarse con cualquiera de las formulaciones de infiltración definidas en la sección 4.1.4.

Modelo agregado

Este modelo es más sencillo que el modelo distribuido. En muchos casos, el modelo agregado puede ser más atractivo de cara al usuario que quiere considerar el flujo subterráneo en el cálculo de forma sencilla, sin necesidad de resolver un modelo complejo que depende de nuevos parámetros variables espacialmente y, en general, difíciles de caracterizar.

El modelo agregado realiza un balance global de agua en el subsuelo de la cuenca considerando los aportes por infiltración desde el terreno y las pérdidas por percolación, evapotranspiración y exfiltración. La componente de exfiltración representa la salida del agua subterránea hacia la red de cauces fluviales, mientras que las componentes de percolación y evapotranspiración representan agua que se pierde del sistema. De esta forma, el agua infiltrada en la cuenca pasa a formar parte de un depósito subterráneo desde el cual puede percolar a acuíferos inferiores (en cuyo caso sale del sistema y se pierde su seguimiento), perderse por evapotranspiración, o salir a los cauces fluviales de la cuenca (exfiltración), en cuyo caso pasa a aumentar el caudal de los cauces superficiales.

Se puede considerar toda la cuenca como un único depósito, o puede dividirse la superficie de la cuenca en varios depósitos con propiedades diferenciadas. Independientemente del número de depósitos en los que se divida la cuenca, el balance de masa de agua en cada depósito se calcula como:

${\frac {d{V}_{s}}{dt}}={Q}_{inf}-{Q}_{per}-{Q}_{evap}-{Q}_{out}$

(62)

donde ${\textstyle {V}_{s}}$ es el volumen de agua en el depósito (subsuelo), ${\textstyle {Q}_{out}}$ es el caudal que sale desde el subsuelo hacia los cauces fluviales (exfiltración), ${\textstyle {Q}_{inf}}$ es el caudal infiltrado desde el terreno al subsuelo, ${\textstyle {Q}_{per}}$ es el caudal que percola desde el subsuelo hacia capas más profundas del terreno, ${\textstyle {Q}_{evap}}$ es el caudal de agua que se pierde por evapotranspiración.

El caudal de infiltración ( ${\textstyle {Q}_{inf}}$) se calcula a partir de cualquiera de los modelos de infiltración presentados en la sección 4.1.4. Los caudales de percolación ( ${\textstyle {Q}_{per}}$) y de evapotranspiración ( ${\textstyle {Q}_{evap}}$) se calculan a partir de sus respectivos valores por unidad de superficie (formulaciones descritas en la secciones 4.1.5 y 4.1.6), simplemente multiplicándolos por la superficie de la cuenca asociada al depósito correspondiente ( ${\textstyle {S}_{tot}}$).

El caudal exfiltrado ( ${\textstyle {Q}_{out}}$) se calcula según la siguiente expresión:

${Q}_{out}=A\,\left({V}_{s}-B\,{S}_{tot}\right)$

(63)

donde ${\textstyle {S}_{tot}}$ es la superficie de la cuenca asociada al depósito correspondiente y ${\textstyle A}$ y ${\textstyle B}$ son parámetros asociados a dicho depósito. El caudal exfiltrado se aporta a la escorrentía superficial, en tramos de ríos definidos previamente por el usuario.

Para aplicar el modelo agregado el usuario debe definir la extensión espacial de cada depósito (a partir de la cual el modelo calculará ${\textstyle {S}_{tot}}$), los parámetros ${\textstyle A}$ y ${\textstyle B}$ para cada depósito, y la zona de exfiltración de cada depósito (en la cual se repartirá el caudal exfiltrado ${\textstyle {Q}_{out}}$).

Además, el usuario debe definir, para cada depósito, el volumen máximo de agua en el subsuelo a efectos de flujo subsuperficial y el contenido inicial de agua. Estos valores se definen como una profundidad máxima de subsuelo ( ${\textstyle {d}_{max}}$) y una saturación inicial ( ${\textstyle {S}_{r,0}}$). A partir de estos valores se calcula el volumen inicial y máximo de agua en cada depósito como:

${\begin{matrix}{V}_{s,max}={S}_{tot}\,{d}_{max}\\{V}_{s,0}={S}_{r,0}{V}_{s,max}\end{matrix}}$

(64)

donde ${\textstyle {V}_{s,max}}$ es el volumen máximo de agua en el depósito y ${\textstyle {V}_{s,0}}$ el volumen inicial de agua en el depósito.

Modelo distribuido

El modelo de flujo subterráneo distribuido se basa en la ecuación de Boussinesq para flujo 2D en medio poroso saturado no confinado. Para evaluar el flujo de agua subsuperficial y la distribución espacio-temporal del agua en el suelo se resuelve la siguiente ecuación diferencial de conservación de la masa:

$\phi {\frac {\partial b}{\partial t}}={\frac {\partial }{\partial x}}\left(T{\frac {\partial H}{\partial x}}\right)+$${\frac {\partial }{\partial y}}\left(T{\frac {\partial H}{\partial y}}\right)+f-p-e-$$g$

(65)

siendo ${\textstyle \phi }$ la porosidad del suelo, ${\textstyle b}$ el espesor saturado de suelo, ${\textstyle T}$ la transmisividad horizontal, ${\textstyle H}$ la carga hidráulica, ${\textstyle f}$ la infiltración de agua desde el terreno al subsuelo, ${\textstyle p}$ la percolación, ${\textstyle e}$ la evapotranspiración y ${\textstyle g}$ la exfiltración de agua del subsuelo a la superficie.

Dado que en Iber no se resuelve el flujo en dirección vertical, el espesor de la capa de suelo saturada ( ${\textstyle b}$) se relaciona con el contenido de humedad promediado verticalmente en la capa de suelo ( ${\textstyle \theta }$ ), mediante la siguiente relación:

$b=\theta {\frac {{d}_{s}}{\phi }}$

(66)

donde ${\textstyle {d}_{s}}$ es el espesor total de la capa de suelo. Es decir, el valor de ${\textstyle b}$ se obtiene agrupando todo el contenido de agua de la capa de suelo en una capa horizontal situada justo bajo el terreno (Fig. 23). Esta simplificación es únicamente válida para espesores de suelo pequeños.

La transmisividad horizontal se considera variable en función del contenido de humedad del suelo, calculándose mediante la siguiente ecuación:

$T={k}_{h}\,{d}_{s}\,{e}^{-\left({\frac {\phi -\theta }{m}}\right)}$

(67)

donde ${\textstyle {k}_{h}}$ es la permeabilidad saturada horizontal del suelo y ${\textstyle m}$ es un parámetro que en la versión actual se toma constante e igual a 0,1. Los valores de porosidad del suelo ( ${\textstyle \phi }$ ) y de espesor total de la capa de suelo ( ${\textstyle {d}_{s}}$) se heredan del modelo de Green-Ampt, mientras que la humedad del suelo ( ${\textstyle \theta }$ ) es una variable que se va actualizando a medida que se resuelve la ecuación de conservación de la masa de agua en el subsuelo, y que también se comparte con el modelo de Green-Ampt.

En un suelo completamente saturado ( ${\textstyle \theta =\phi }$ ), la transmisividad es máxima e igual a la permeabilidad saturada por el espesor de suelo ${\textstyle T=}$${k}_{h}\,{d}_{s}$.

La carga hidráulica se calcula como la cota inferior de la capa de suelo más el espesor saturado:

$H={z}_{b}-{d}_{s}+b$

(68)

donde ${\textstyle {z}_{b}}$ es la cota del terreno.

Cuando se utiliza el modelo de flujo subsuperficial distribuido, la tasa de infiltración ( ${\textstyle f}$) debe calcularse siempre con la formulación de Green-Ampt descrita anteriormente, de la cual hereda los siguientes parámetros: porosidad del suelo ( ${\textstyle \phi }$ ), espesor total de la capa de suelo ( ${\textstyle {d}_{s}}$) y humedad inicial del suelo ( ${\textstyle {\theta }_{0}}$). La permeabilidad saturada del suelo utilizada para el cálculo de la transmisividad horizontal ( ${\textstyle {k}_{h}}$) no se hereda del modelo de Green-Ampt ya que la infiltración es un fenómeno fundamentalmente en la dirección vertical, mientras que el flujo subsuperficial representado por la ecuación de Boussinesq se produce fundamentalmente en dirección horizontal (paralela al terreno). Debido a la anisotropía de muchos suelos, se permite utilizar un valor diferente de permeabilidad saturada en el modelo de Green-Ampt y en la ecuación de Boussinesq. Evidentemente, el usuario puede utilizar el mismo valor si lo considera oportuno.

Las tasas de evapotranspiración ( ${\textstyle e}$) y percolación ( ${\textstyle p}$) se calculan según las formulaciones descritas en las secciones 4.1.5 y 4.1.6 para el modelo de Green-Ampt. La exfiltración ( ${\textstyle e}$) se calcula imponiendo que el espesor de la capa de suelo saturada debe ser inferior o igual al espesor total de la capa de suelo ( ${\textstyle b\leq {d}_{s}}$), lo cual equivale a decir que la humedad del suelo no puede ser superior a su porosidad ( ${\textstyle \theta \leq \phi }$ ). De esta forma, si en algún momento del cálculo se obtiene un valor de ${\textstyle b}$ superior a ${\textstyle {d}_{s}}$, la diferencia entre ambos se impone como un caudal de exfiltración que sale del subsuelo hacia los cauces fluviales, al mismo tiempo que el valor de ${\textstyle b}$ se hace igual a ${\textstyle {d}_{s}}$.

La ecuación de conservación de la masa de agua en el subsuelo se resuelve mediante el método de volúmenes finitos en la misma malla no estructurada en la que se resuelven las ecuaciones de aguas someras.

Fig. 24. Esquema de los procesos considerados en el modelo distribuido de flujo subsuperficial.

En cuanto a la implementación y acoplamiento del flujo subsuperficial y del flujo de escorrentía superficial, el proceso de cálculo (Fig. 24) en cada paso de tiempo es el siguiente:

1. Resolver las ecuaciones de aguas someras 2D (escorrentía superficial) sin considerar aportes de lluvia ni de infiltración. Actualización de los calados y velocidades del agua en superficie.

2. Resolver la ecuación de flujo subsuperficial 2D. Actualización de la humedad media en la capa de suelo ( ${\textstyle \theta }$ ) y del espesor saturado equivalente ( ${\textstyle b}$).

3. Si el espesor saturado es superior al espesor de la capa de suelo ( ${\textstyle b\geq \,{d}_{s}}$) se está produciendo exfiltración (agua subsuperficial saliendo a la superficie), en cuyo caso se actualizan las variables de flujo de la siguiente manera:

• ${\textstyle b\,=\,{d}_{s}}$
• ${\textstyle h\,=\,h\,+\,\left(b-{d}_{s}\right)\,\phi }$
• El caudal de exfiltración en este caso sería igual a ${\textstyle g=}$${\frac {\left(b-{d}_{s}\right)\,\phi }{\Delta t}}$

4. Añadir el aporte por lluvia

• ${\textstyle h\,=\,h\,+\,R\,\Delta t}$

5. Añadir las pérdidas por infiltración según el modelo de Green-Ampt

• ${\textstyle {f}_{pot}={k}_{s}\left(1+{\frac {\left(h+\psi \right)}{b}}\right)}$
• ${\textstyle {F}_{max}=\mathrm {min} \,\left[h,\,\left({d}_{s}-b\right)\phi \right]}$
• ${\textstyle F=\mathrm {min} \,\left[{f}_{pot}\,\Delta t,\,{F}_{max}\right]\,}$
• ${\textstyle h=h-F}$
• ${\textstyle b=b+F/\phi }$

6. Añadir las pérdidas por evapotranspiración y por percolación

4.1.8 Resultados de las variables hidrológicas

Además de las variables hidrodinámicas de calado y velocidad del agua superficial (escorrentía superficial), al activar los procesos hidrológicos se escriben los siguientes resultados:

• Hidrogramas: el usuario puede definir diferentes contornos de salida de agua, así como secciones de aforo internas. Los hidrogramas calculados a través de dichos contornos de salida y de las secciones internas de aforo se escribirán como series temporales en el fichero Q_hydrograph.rep.
• Series temporales de variables hidrológicas: las series temporales de las siguientes variables se escriben en el fichero Hydrology_timeseries.rep:

• Profundidad de escorrentía superficial (mm) instantánea y promediada en toda la cuenca
• Profundidad de precipitación (mm) acumulada y promediada en toda la cuenca
• Intensidad de precipitación (mm/h) instantánea y promediada en toda la cuenca
• Profundidad de infiltración (mm) acumulada y promediada en toda la cuenca
• Tasa de infiltración (mm/h) instantánea y promediada en toda la cuenca
• Volumen de agua en el subsuelo (m³) instantáneo
• Espesor saturado del subsuelo (mm) instantáneo y promediado en toda la cuenca
• Saturación del subsuelo (-) instantánea y promediada en toda la cuenca
• Profundidad de evapotranspiración (mm) acumulada y promediada en toda la cuenca
• Tasa de evapotranspiración (mm/h) instantánea y promediada en toda la cuenca
• Profundidad de percolación (mm) acumulada y promediada en toda la cuenca
• Tasa de percolación (mm/h) instantánea y promediada en toda la cuenca
• Caudal exfiltrado (m³/s) instantáneo

• Campos espaciales de variables hidrológicas: en la interfaz de post-proceso, se muestra un nuevo grupo de variables (Hydrology) en el que se incluyen los campos espaciales de las siguientes variables:

• Intensidad de precipitación (mm/h)
• Profundidad de precipitación (mm) acumulada desde el inicio del cálculo
• Tasa de infiltración (mm/h)
• Profundidad de infiltración (mm) acumulada desde el inicio del cálculo
• Espesor saturado del suelo (m)
• Altura piezométrica en el subsuelo (m)
• Caudal unitário subsuperficial (m²/s)

4.2 Interfaz de usuario

En esta sección se describen los diferentes menús y ventanas desde las que se pueden activar los procesos hidrológicos descritos en las secciones anteriores y definir los correspondientes parámetros.

4.2.1 Precipitación

La precipitación puede definirse de forma variable en espacio y tiempo, ya sea mediante asignación manual de hietogramas a diferentes zonas del modelo o de forma automática mediante la definición de estaciones pluviométricas o de rasters de precipitación.

La activación y el modo de la lectura pluviométrica se realiza a través del menú Datos >> Datos del Problema >> Hidrología >> Precipitación (Fig. 25). Dentro de dicho menú, se dispone de 4 opciones:

• Sin precipitación
• Hietograma
• Rasters
• Interpolación raster

Aunque se definan y se asignen los hietogramas o los rasters de lluvia, éstos nos serán utilizados en el cálculo si se mantiene la opción “Sin precipitación”.

Fig. 25. Ventana de selección del modo de lectura de datos pluviométricos.

Lluvia definida por hietogramas

La definición de hietogramas se realiza a través del menú Datos >> Procesos hidrológicos >> Lluvia >> Definición de hietogramas (Fig. 26a). En la ventana correspondiente, el usuario puede introducir las series temporales de datos pluviométricos a través del instante temporal (s) y la precipitación (mm/h). Además, es necesario establecer la ubicación de la estación pluviométrica en la que se registran dichos datos (Fig. 26a). La ubicación de la estación se selecciona a través de la interfaz de GiD; para indicar la ubicación fuera de la superficie del modelo, es necesario escribir las coordenadas que se desean manualmente.

Fig. 26. Definición y asignación de hietogramas: (a) ventana de definición de hietogramas; (b) ventana de selección del modo de lectura de datos pluviométricos.

Una vez definidos los hietogramas, el usuario dispone de dos opciones:

• Se puede proceder a la asignación manual de cada hietograma en diferentes regiones del modelo, a través del menú Datos >> Procesos hidrológicos >> Lluvia >> Asignación de hietogramas (Fig. 26b). De esta manera, se asigna la precipitación definida en el hietograma de forma homogénea en las regiones del modelo seleccionadas. Si ya se ha generado la malla, también es posible asignar los hietogramas directamente a los elementos de la malla.
• En las regiones en las cuales no se haya asignado manualmente un hietograma, el programa atribuirá automáticamente la pluviometría a cada zona del modelo a través de una interpolación del tipo Nearest Neighbour, considerando para ello la ubicación de cada pluviómetro.

Lluvia definida por rasters

Si seleccionamos esta opción, la lluvia se definirá a partir de archivos raster en formato ASCII Grid de ArcInfo (Fig. 27). En cada archivo se definirá un campo de precipitación distribuido espacialmente. Todos los rasters deben figurar en una carpeta llamada Rain_rasters, ubicada dentro de la carpeta del proyecto, para garantizar su correcta lectura. Se aclara, por tanto, que esta carpeta debe ser creada por el usuario y que no se genera por defecto con el modelo. Además, todos los ficheros rasters deben tener el mismo número de columnas y filas, el mismo tamaño de pixel y el mismo origen.

Dentro de la definición de la lluvia por rasters se tienen dos alternativas:

• Rasters. Se toma el mismo valor de lluvia entre cada dos rasters, constante en el tiempo. Esta es la forma análoga a la definición temporal de la lluvia en hietogramas, i.e. la variación temporal de la precipitación consiste en bloques de intensidad constante. Si se selecciona esta opción es posible elegir el tipo de datos de lluvia incluido en los rasters, que pueden ser valores de intensidad de precipitación (mm/h) o de profundidad de precipitación (mm) entre cada dos rasters.
• Rasters interpolados. Se realiza una interpolación lineal en el tiempo de la precipitación a partir de los datos de dos archivos raster consecutivos. En este caso, los rasters presentan siempre valores de intensidad de lluvia (mm/h).

En cualquiera de los casos anteriores, el usuario debe definir una tabla en la que se asigna cada archivo raster a un instante de tiempo. Los archivos raster pueden tener cualquier nombre, y en la tabla se debe incluir el nombre completo del archivo (incluyendo la extensión).

Fig. 27. Ventana de definición de series temporales de archivos raster de precipitación.

4.2.2 Modelo y parámetros de infiltración

Iber permite la asignación de parámetros de infiltración de manera manual o automática a partir de un archivo raster. En ambos casos se permite elegir entre diferentes modelos de infiltración (Lineal, SCS, Horton, Green-Ampt, SCS continuo) y los parámetros de cada modelo pueden ser variables espacialmente.

Como siempre en Iber, la asignación automática de propiedades se realiza directamente sobre la malla de cálculo, que por lo tanto deberá estar creada con anterioridad. Asimismo, cada vez que se cree la malla de cálculo, deben realizarse todas las asignaciones automáticas de nuevo.

La elección del modelo de pérdidas por infiltración se realiza desde el menú Datos >> Datos del Problema >> Hidrología >> Pérdidas por infiltración (Fig. 28). Los parámetros relativos a cada modelo de infiltración se pueden asignar de forma manual o automática, según se detalla a continuación.

Fig. 28. Ventana de selección del modelo de infiltración.

Asignación manual de parámetros de infiltración

Para realizar una asignación manual de los parámetros de infiltración se debe seleccionar la opción Activar en Datos >> Datos del Problema >> Hidrología >> Pérdidas por infiltración. Una vez hecho esto aparecerán los modelos de infiltración disponibles.

Una vez seleccionado un modelo, la asignación manual de parámetros se realiza desde el menú Datos >> Procesos Hidrológicos >> Pérdidas por infiltración >> Asignación manual. En esta ventana el usuario puede introducir los parámetros correspondientes a cada modelo de infiltración (Fig. 29a) y asignarlos manualmente a diferentes superficies de la geometría o elementos de la malla. Desde la ventana que se muestra en la Fig. 29a pueden asignarse los parámetros de cualquier modelo de infiltración, pero Iber calculará las pérdidas únicamente con el modelo seleccionado en Datos >> Datos del Problema >> Hidrología >> Pérdidas por infiltración.

Fig. 29. (a) Ventana de selección manual del modelo de infiltración. (b) Ventana de definición de automática de zonas de infiltración.

Asignación automática a partir de archivos raster con los valores de cada parámetro

Iber permite la asignación automática de parámetros de infiltración a partir de archivos raster georreferenciados en formato ASCII Grid de ArcInfo. Se debe preparar un archivo raster por cada parámetro, que defina la distribución espacial del valor de dicho parámetro en la zona de estudio.

Para realizar la asignación automática a partir de dichos archivos raster se debe seleccionar la opción Activar en Datos >> Datos del Problema >> Hidrología >> Pérdidas por infiltración. Una vez hecho esto, aparecerán los modelos de infiltración disponibles. Una vez seleccionado un modelo de infiltración, la asignación automática de parámetros se realiza desde el menú Datos >> Procesos Hidrológicos >> Pérdidas por infiltración >> Asignación automática por parámetros.

Se presenta, en primer lugar, la opción de escoger el modelo de infiltración para el cual queremos asignar parámetros. La asignación automática puede realizarse para los modelos Lineal, SCS y Green-Ampt. Una vez que se elige el modelo de infiltración, se abre un directorio para escoger el archivo raster correspondiente al parámetro que se introducirá. A continuación, se describen los casos posibles:

a) Si se elige el modelo SCS, únicamente se pedirá cargar un archivo raster, correspondiente al número de curva. En dicho archivo se especificarán directamente en cada celda los valores de número de curva correspondientes.

b) Si se elige el modelo Lineal, se pedirá cargar dos archivos raster en el siguiente orden:

• Archivo correspondiente a las pérdidas iniciales (mm)
• Archivo de tasas de infiltración potencial constantes en tiempo (mm/h)

c) Si se elige el modelo Green-Ampt, se pedirá cargar seis archivos raster, en el siguiente orden:

• Succión en región no saturada (mm)
• Porosidad total
• Saturación inicial
• Permeabilidad saturada del suelo en la dirección vertical (mm/h)
• Pérdidas iniciales (mm)
• Espesor total de la capa de suelo (m)

Al igual que en la asignación manual de parámetros, pueden asignarse de forma automática los parámetros de cualquier modelo de infiltración, pero Iber calculará las pérdidas únicamente con el modelo seleccionado en Datos >> Datos del Problema >> Hidrología >> Pérdidas por infiltración.

Asignación automática de parámetros de infiltración por zonas

Adicionalmente a la asignación automática de parámetros de infiltración descrita en la sección anterior, Iber también permite asignar los parámetros de infiltración de forma distribuida espacialmente por zonas homogéneas. De esta forma, en lugar de definir cada parámetro por separado, se definen manualmente zonas con un conjunto de parámetros homogéneos (en cada zona) y posteriormente se asignan estas zonas al modelo mediante un archivo raster en el cual se incluye la distribución espacial de dichas zonas (y no de cada parámetro independiente). Esta forma de asignar parámetros se llama en Iber Asignación automática por zonas.

Para asignar los parámetros de esta manera, en primer lugar, es necesario seleccionar la opción Encendido [por zonas] en Datos >> Datos del Problema >> Hidrología >> Pérdidas por infiltración. Una vez hecho esto aparecerán los modelos de infiltración disponibles, debiendo seleccionar uno de ellos.

Una vez seleccionado el modelo de infiltración, se definen las zonas de infiltración desde el menú Datos >> Procesos Hidrológicos >> Pérdidas por infiltración >> Asignación automática por zonas >> Definición de zonas de infiltración. Desde aquí se accede a una ventana en la que se pueden crear tantas zonas como se quiera. Para cada zona hay que seleccionar el modelo de infiltración y los parámetros correspondientes (Fig. 29b). Una vez definidas las zonas de infiltración se deben asignar al modelo mediante el menú Datos >> Procesos Hidrológicos >> Pérdidas por infiltración >> Asignación automática por zonas >> Asignación de zonas de infiltración. Al seleccionar esta opción, se abrirá una ventana desde la cual será necesario seleccionar un archivo raster en el que se incluya la distribución espacial de las zonas de infiltración. De esta forma, cada zona del modelo tendrá asignada una zona de infiltración con sus parámetros correspondientes.

4.2.3 Evapotranspiración

Solo puede activarse la evapotranspiración si se activa alguna de las siguientes opciones en Datos >> Datos del Problema >> Hidrología:

• Modelo de infiltración de Green-Ampt
• Modelo de infiltración de SCS continuo
• Modelo de flujo subterráneo agregado

Desde esta ventana también se puede determinar el coeficiente de corrección para el cálculo de la evapotranspiración potencial (comprendido entre 0,7 y 2) y, en el caso de haber elegido el modelo de infiltración de Green-Ampt, activarse también la formulación de estrés hídrico (Soil Water Stress).

Para el cálculo de la evapotranspiración potencial es necesario definir la temperatura del aire. La definición de series temporales de temperatura y su asignación a diferentes zonas del modelo puede realizarse desde los menús Datos >> Atmosféricas >> Variables Atmosféricas y Datos >> Atmosféricas >> Asignación Var Atmosféricas.

Estrés hídrico

La formulación de estrés hídrico (Soil Water Stress) solo está disponible en el caso de que se seleccione Green-Ampt como modelo de infiltración y se active la evapotranspiración. Éstas se seleccionan desde Datos >> Datos del Problema >> Hidrología.

Para utilizar esta formulación, una vez activada en Datos >> Datos del Problema >> Hidrología, es necesario definir uno o varios tipos de cultivo desde el menú Datos >> Procesos Hidrológicos >> Estrés hídrico del suelo >> Definición de cultivo. Desde aquí se accede a la ventana de definición de cultivos, en la que se pueden crear tantos cultivos como se quiera. Para cada cultivo creado hay que definir los siguientes parámetros: capacidad de campo, punto de marchitez, factor de agotamiento y ${\textstyle {K}_{c}}$. Dado que ${\textstyle {K}_{c}}$ es un coeficiente que depende del tipo de cultivo y de su estado, y por lo tanto puede ser variable en función de la época del año, el programa permite definir una serie temporal de dicho parámetro, con el fin de poder incluir su variación estacional.

Una vez definidos los cultivos que sean necesarios, se deben asignar al modelo mediante el menú Datos >> Procesos Hidrológicos >> Estrés hídrico del suelo >> Asignación de cultivo. Esta opción da acceso a la ventana de asignación de cultivos, desde la cual podemos seleccionar cualquiera de los cultivos creados y asignarlo manualmente a diferentes zonas del modelo.

4.2.4 Percolación

Solo puede activarse la percolación si se activa alguna de las siguientes opciones en Datos >> Datos del Problema >> Hidrología:

• Modelo de infiltración de Green-Ampt
• Modelo de infiltración de SCS continuo
• Modelo de flujo subterráneo agregado

El cálculo de la percolación se activa o desactiva a través del menú Datos >> Datos del Problema >> Hidrología, donde también puede determinarse el índice de distribución de poros (comprendido entre 0,05 y 0,5) y la permeabilidad para el cálculo de la percolación.

4.2.5 Flujo subsuperficial

Iber dispone de dos modelos para evaluar el flujo de agua a través del subsuelo: un modelo agregado y un modelo distribuido. El modelo distribuido solo puede utilizarse si se activa el modelo de Green-Ampt para calcular la infiltración, ya que hereda sus parámetros de dicho modelo. El modelo agregado puede utilizarse con cualquiera de las formulaciones de infiltración definidas en la sección 4.1.4.

En todo caso, es necesario tener activado algún modelo de infiltración para poder calcular el flujo subsuperficial. La elección del modelo de flujo subsuperficial se realiza en el menú Datos >> Datos del Problema >> Hidrología.

Si se selecciona el modelo de infiltración de Green-Ampt, en la ventana Datos >> Datos del Problema >> Hidrología aparecerá una nueva opción llamada Flujo Subsuperficial que permite elegir entre el modelo agregado o el distribuido (o desactivar el flujo subsuperficial). Si se selecciona otro modelo de infiltración (SCS, Lineal o Horton) aparecerá una opción llamada Groundwater flow que únicamente permitirá seleccionar el modelo agregado (o desactivar el flujo subsuperficial).

Una vez seleccionado un modelo de flujo subsuperficial, la asignación de los parámetros correspondientes a cada modelo se realiza a través de los menús Datos >> Procesos hidrológicos >> Flujo subsuperficial >> Modelo agregado y Datos >> Procesos hidrológicos >> Flujo subsuperficial >> Modelo distribuido.

Modelo agregado de flujo subsuperficial

Para aplicar el modelo agregado el usuario debe definir la extensión espacial de cada depósito subterráneo, los parámetros de cada depósito (constantes ${\textstyle A}$ y ${\textstyle B}$, profundidad máxima y saturación inicial), y la zona de exfiltración de cada depósito.

La creación de depósitos y la definición de sus parámetros se realiza a través del menú Datos >> Procesos hidrológicos >> Flujo subsuperficial >> Modelo agregado >> Definición de depósitos, desde donde se pueden crear tantos depósitos como queramos con parámetros diferentes. Las superficies de infiltración y exfiltración de dichos depósitos pueden asignarse a través de Datos >> Procesos hidrológicos >> Flujo subsuperficial >> Modelo agregado >> Asignación de depósitos, haciendo referencia a los depósitos creados anteriormente. Es necesario aclarar que se debe asignar tanto las áreas de infiltración, como de exfiltración correspondientes a cada depósito.

Modelo distribuido de flujo subsuperficial

El modelo distribuido de flujo subsuperficial hereda la mayoría de sus parámetros del modelo de Green-Ampt. Únicamente es necesario definir adicionalmente la permeabilidad saturada horizontal del suelo, lo cual puede hacerse desde Datos >> Procesos hidrológicos >> Flujo subsuperficial >> Modelo distribuido >> Parámetros (Fig. 30a). De esta forma se puede asignar de forma manual una permeabilidad saturada a diferentes zonas del modelo (i.e. distribuida espacialmente).

Asimismo, es necesario definir los contornos por donde puede entrar o salir el agua de la zona de estudio a través del subsuelo. Si no se define ninguna condición de contorno para el modelo distribuido el agua no podrá salir del dominio de estudio a través del subsuelo (únicamente a través de la superficie). En dichos contornos se puede imponer como condición la descarga unitaria o la carga hidráulica. Dichas condiciones de contorno se definen a través del menú Datos >> Procesos hidrológicos >> Flujo subsuperficial >> Modelo distribuido >> Condiciones de contorno, siendo necesario asignarlas a los contornos del modelo correspondientes (Fig. 30b).

Fig. 30. (a) Ventana de asignación de los parámetros del modelo distribuido. (b) Ventana de asignación de condiciones de contorno del modelo distribuido.

Las condiciones iniciales de humedad del subsuelo se heredan de los parámetros definidos en el modelo de Green-Ampt. Concretamente se utiliza la saturación inicial y la porosidad definidas en Green-Ampt para determinar la humedad inicial del subsuelo.

4.2.6 Secciones de aforo y series temporales

Iber permite definir secciones de aforo y sondas temporales, a través de las cuales obtener salidas de resultados de forma localizada. En las secciones de aforo, Iber calcula el caudal que pasa a través de ellas, generando una salida de hidrogramas en cada sección de aforo definida antes del cálculo.

Las secciones de aforo se definen a través del menú Datos >> Hidrodinámica >> Condiciones Internas. En esta ventana, en la opción Tipo de condición deberemos elegir Estación de aforo (Fig. 31a). Es necesario asignarle un número de estación diferente a cada estación de aforo.

Los resultados se escriben en forma de series temporales en el fichero Q_hydrograph.rep (Fig. 31b). Este archivo que también incluye los hidrogramas correspondientes a otras condiciones internas que se hayan definido en el modelo, como puede ser el caso de puentes, vertederos y/o compuertas.

La definición de sondas temporales se realiza desde el menú Datos >> Datos del Problema >> Resultados >> Sondas. Desde esta ventana el usuario puede introducir las coordenadas de las sondas, así como el intervalo temporal de salida de resultados (que no tiene por qué coincidir con el intervalo de salida temporal de resultados en toda la malla de cálculo). En las sondas definidas se almacenan las variables hidráulicas interpoladas linealmente a partir de su valor en las celdas de la malla de cálculo. Las series temporales resultantes se escriben en los ficheros res_gauge_NN.rep, donde NN es el número de la sonda.

Fig. 31. (a) Ventana de la función de asignación de secciones de aforo. (b) Ejemplo de formato de archivo Q_hydrograph.rep.

5 Módulo distribuido de erosión de suelos en laderas por escorrentía

5.1 Fundamentos teóricos

El módulo de erosión de suelos es un modelo distribuido y bidimensional, que resuelve una serie de ecuaciones con base física dependientes de variables utilizadas habitualmente en hidrología. La ecuación principal es la de conservación de la masa de sedimento en el suelo, la cual se resuelve en cada elemento de la malla de cálculo y para cada clase de sedimento que se defina (en el caso de granulometrías distribuidas). Además, se resuelve la ecuación de conservación de sedimento en suspensión, también para cada clase de sedimento. Ambas ecuaciones están relacionadas por un término que modela la sedimentación y resuspensión de partículas de sedimento.

5.1.1 Procesos considerados

En la Fig. 32 se representa de forma esquemática la conceptualización de la estructura del suelo que se realiza en el módulo de erosión de suelos. La capa superior representa el agua de escorrentía (o corriente de un río) que transporta sedimento en suspensión, la cual se extiende entre la cota del terreno ( ${\textstyle {z}_{b}}$) y la superficie de la lámina de agua ( ${\textstyle {z}_{s}}$), con un calado ${\textstyle h=}$${z}_{s}-{z}_{b}$. La concentración de sedimentos en suspensión en dicha capa (promediada en profundidad) es igual a ${\textstyle C}$, y en ella existe un transporte de sedimentos en suspensión igual ${\textstyle {q}_{s}=}$$C\,h\,U$.

Por debajo de la capa de escorrentía existe una primera capa de suelo formada por sedimento no-cohesivo que ya ha sido erosionada previamente de la matriz original de suelo. Las partículas de sedimento de dicha capa pueden ponerse en suspensión y pasar a la capa superior, ya sea por el impacto de las gotas de lluvia o por la fricción de fondo ejercida por la escorrentía superficial. Estos dos procesos se representan en la Fig. 32 por los términos ${\textstyle {D}_{rdrd}}$ y ${\textstyle {D}_{fdre}}$, respectivamente. Al mismo tiempo, una parte de las partículas de suelo que se encuentran en suspensión pueden depositarse en la capa de suelo no-cohesivo debido a su peso. Este proceso se representa en la Fig. 32 por el término ${\textstyle {D}_{dep}}$.

Debajo de la capa de suelo no-cohesivo puede existir una segunda capa de suelo correspondiente a la matriz original de suelo que aún no ha sido erosionado. Suele corresponder a un suelo con fuerzas de cohesión importantes, por ello nos referiremos a dicha capa como suelo cohesivo, por contraposición a la capa de suelo no-cohesivo. El impacto de las gotas de lluvia y la fuerza de fricción de fondo generada por la escorrentía superficial pueden poner en suspensión el sedimento de la capa cohesiva. Estos procesos se representan en la Fig. 32 por los términos ${\textstyle {D}_{rdd}}$ y ${\textstyle {D}_{fde}}$ respectivamente, y son análogos a los términos ${\textstyle {D}_{rdrd}}$ y ${\textstyle {D}_{fdre}}$.

Fig. 32. Esquema de la conceptualización de la estructura del suelo considerada en el módulo de erosión de suelos.

Nótese en la Fig. 32 que el proceso de deposición ( ${\textstyle {D}_{dep}}$) siempre aporta sedimento a la capa no-cohesiva, ya que no existen fuerzas de cohesión entre las partículas en suspensión que se depositan.

En cada una de las dos capas de suelo pueden existir distintas clases de sedimento, en general diferenciadas entre ellas por su diámetro. Esto permite considerar granulometrías distribuidas o diferentes tipos de suelo. La proporción de cada clase de sedimento en la capa de suelo erosionada se denotará con el factor ${\textstyle {f}_{k}}$, que representa el tanto por uno de la clase ${\textstyle k}$ en la capa no-cohesiva. Al representar la proporción de cada clase con respecto al sedimento total, la suma de estos factores debe ser igual a 1:

$\sum _{k=1}^{N}{f}_{k}=1\,$

(69)

siendo ${\textstyle N}$ el número de clases de sedimento consideradas. La distribución de clases de sedimento en la capa erosionada (factores ${\textstyle {f}_{k}}$) se actualiza durante el cálculo en función del sedimento erosionado y depositado en dicha capa.

5.1.2 Definición de variables utilizadas en varias formulaciones

Las siguientes variables se utilizan en varias de las formulaciones implementadas en el modelo:

• Densidad relativa sumergida de las partículas del sedimento:

$R={\frac {{\rho }_{s}-\rho }{\rho }}$

(70)

• Diámetro adimensional de las partículas de la clase ${\textstyle k}$:

${D}_{\ast ,k}={D}_{k}{\left({\frac {Rg}{{\nu }^{2}}}\right)}^{0,33}$

(71)

• Velocidad de fricción de fondo:

${u}_{\ast }={\sqrt {\frac {{\tau }_{b}}{\rho }}}$

(72)

• Tensión de fondo adimensional:

${\tau }_{k}^{\ast }={\frac {{u}_{\ast }^{2}}{Rg{D}_{k}}}={\frac {\tau }{\rho Rg{D}_{k}}}$

(73)

• Velocidad de sedimentación de las partículas según van Rijn [37]:

${\begin{matrix}{D}_{k}<{10}^{-4}m&\rightarrow &{w}_{s,k}={\frac {Rg{D}_{k}^{2}}{18\,\nu }}\\{10}^{-4}m<{D}_{k}<{10}^{-3}m&\rightarrow &{w}_{s,k}={\frac {10\,\nu }{{D}_{k}}}\left({\sqrt {1+0,01{D}_{\ast ,k}^{3}}}-1\right)\\{D}_{k}>{10}^{-3}m&\rightarrow &{w}_{s,k}=1,1\,{\sqrt {Rg{D}_{k}}}\end{matrix}}$

(74)

siendo ${\textstyle \nu }$ la viscosidad cinemática del agua ( ${\textstyle {10}^{-6}}$ m²/s).

5.1.3 Ecuaciones de la hidrodinámica

Para definir los campos de profundidad y velocidad se utilizan las ecuaciones de aguas poco profundas bidimensionales, incluyendo los términos que representan las componentes de precipitación e infiltración (Ecuación (1)). La simulación de la fricción de fondo se realiza mediante la inclusión de la ecuación de Manning.

5.1.4 Transporte por carga en suspensión

Conservación de la masa de sedimento en suspensión

El cálculo de la concentración de sedimento en suspensión se realiza resolviendo la siguiente ecuación de transporte en suspensión promediada en profundidad, para cada clase ( ${\textstyle k}$) de sedimento considerada en el modelo:

${\frac {\partial h{C}_{k}}{\partial t}}+\left({\frac {\partial {q}_{x}\,{C}_{k}}{\partial x}}+\right.$$\left.{\frac {\partial {q}_{y}\,{C}_{k}}{\partial y}}\right)={D}_{rdd,k}+{D}_{rdrd,k}+$${D}_{fde,k}+{D}_{fdre,k}-{D}_{dep,k}$

(75)

donde ${\textstyle {C}_{k}}$ [kg/m³] es la concentración en suspensión de la clase de sedimento ${\textstyle k}$. Es decir, se resuelve una ecuación de transporte en suspensión para cada clase de sedimento. Los términos a la derecha de la igualdad representan cada uno de los procesos de puesta en suspensión y deposición considerados en el modelo (Fig. 32). Sus unidades son [kg/m²/s]:

• ${\textstyle {D}_{rdd}}$ representa la puesta en suspensión de partículas de suelo procedentes de la capa cohesiva debido al impacto de las gotas de lluvia (rainfall-driven detachment).
• ${\textstyle {D}_{rdrd}}$ representa la puesta en suspensión de partículas de suelo procedentes de la capa no-cohesiva debido al impacto de las gotas de lluvia (rainfall-driven re-detachment).
• ${\textstyle {D}_{fde}}$ representa la puesta en suspensión de partículas de suelo procedentes de la capa cohesiva debido a la fricción de fondo ejercida por la escorrentía superficial (flow-driven entrainment).
• ${\textstyle {D}_{fdre}}$ representa la puesta en suspensión de partículas de suelo procedentes de la capa no-cohesiva debido a la fricción de fondo ejercida por la escorrentía superficial (flow-driven re-entrainment).
• ${\textstyle {D}_{dep}}$ representa la deposición de partículas en suspensión en la capa no-cohesiva, debido a su velocidad de sedimentación.

A continuación, se detallan las formulaciones empleadas para calcular los términos anteriores.

Producción de sedimento en suspensión por impacto de lluvia

Los términos de puesta en suspensión de sedimento debido al impacto de las gotas de lluvia ( ${\textstyle {D}_{rdd}}$ y ${\textstyle {D}_{rdrd}}$) expresados en [kg/m²/s] se calculan como:

${\begin{matrix}{D}_{rdd,k}={\alpha }_{d,k}\,R\,{\sigma }_{h}\left(1-\,\epsilon \right)&\rightarrow &{D}_{rdrd,k}={\alpha }_{rd,k}\,R\,{\sigma }_{h}\,\epsilon \,{f}_{k}\end{matrix}}$

(76)

donde ${\textstyle R}$ es la intensidad de lluvia [m/s], ${\textstyle {\alpha }_{d,k}}$ y ${\textstyle {\alpha }_{rd,k}}$ [kg/m³] son los coeficientes de erosionabilidad para la capa cohesiva y no-cohesiva respectivamente (detachment y redetachment), ${\textstyle {\sigma }_{h}}$ es un factor que considera el efecto protector de la lámina de escorrentía sobre el suelo, y ${\textstyle \epsilon }$ es un factor adimensional de escudo que representa el efecto protector que la capa ya erosionada (no-cohesiva) ejerce sobre la matriz original de suelo (capa cohesiva). El término ${\textstyle {D}_{rdrd,k}}$ incluye la proporción de la clase de sedimento ${\textstyle k}$ en la capa erosionada. De esta forma, si a lo largo del cálculo en algún momento se erosiona toda una clase de sedimento, su proporción será igual a cero ( ${\textstyle {f}_{k}=}$$0$) y, por lo tanto, dicha clase no contribuirá a generar más sedimento en suspensión ( ${\textstyle {D}_{rdrd,k}=}$$0$).

En las ecuaciones anteriores, también se puede expresar la lluvia en [mm/h] y los coeficientes de erosionabilidad en [(kg/s/m²)/(mm/h)]. De esta forma se puede interpretar el coeficiente de erosionabilidad como los kg de sedimento erosionado por unidad de tiempo (segundo), por unidad de superficie (m²) y por unidad de intensidad de precipitación (mm/h).

El factor ${\textstyle {\sigma }_{h}}$ controla el efecto protector de la lámina de escorrentía sobre la erosión del suelo. Si el calado de la lámina de escorrentía es elevado, las gotas de lluvia no serán capaces de poner el sedimento en suspensión ya que su energía será disipada por la lámina de agua. El término ${\textstyle {\sigma }_{h}}$ se calcula como:

${\sigma }_{h}={\left({\frac {{z}_{m}}{h}}\right)}^{0,8}\quad \,\,{\sigma }_{h}\in \left(0,1\right)$

(77)

siendo ${\textstyle {z}_{m}}$ un calado umbral por debajo del cual no existe efecto protector, que debe ser definido como parámetro de entrada del modelo. En efecto, si el calado de la escorrentía es inferior a dicho calado umbral ( ${\textstyle h<{z}_{m}}$) el valor de ${\textstyle {\sigma }_{h}=}$$1$. A medida que el calado aumenta por encima de ${\textstyle {z}_{m}}$, el valor de ${\textstyle {\sigma }_{h}}$ comienza a disminuir, tendiendo a cero para calados elevados.

Por su parte, el coeficiente protector ${\textstyle \epsilon }$ , cuyo valor está comprendido entre 0 y 1, se puede interpretar como un factor de reparto de la energía de la lluvia entre la capa cohesiva y la no cohesiva. Se calcula considerando una relación lineal en base a la masa de sedimento por unidad de superficie en la capa erosionada ( ${\textstyle {M}_{s}}$) y a un valor crítico de dicha masa de sedimento ( ${\textstyle {M}_{s,cr}}$), expresados ambos en [kg/m²]:

$\,\epsilon =min\left[{\frac {{M}_{s}}{{M}_{s,cr}}},1\right]$

(78)

Cuando ${\textstyle {M}_{s}}$ alcanza o supera el valor crítico ( ${\textstyle {M}_{s}\geq {M}_{s,cr}}$), el efecto protector es máximo ( ${\textstyle \epsilon =}$$1$), siendo nulo el material erosionado en la capa cohesiva ( ${\textstyle {D}_{rdd}=}$$0$). Cuando no existe material en la capa no-cohesiva ( ${\textstyle {M}_{s}=}$$0)$ el efecto protector es nulo ( ${\textstyle \epsilon =0}$).

Producción de sedimento en suspensión por escorrentía superficial (fricción de fondo)

El término de puesta en suspensión desde la capa cohesiva por efecto de la fricción de fondo ( ${\textstyle {D}_{fde}}$) se calcula como:

${D}_{fde,k}=\left\{{\begin{matrix}{K}_{k}\,\left(\tau -{\tau }_{crc}\right)\left(1-\epsilon \right)&si&\,{\tau >\tau }_{cr}\\0&si&{\tau \leq \tau }_{cr}\end{matrix}}\right.$

(79)

donde ${\textstyle \tau }$ en [N/m²] es la tensión de fondo, ${\textstyle {\tau }_{crc}}$ en [N/m²] es la tensión crítica de la matriz original de suelo (cohesiva), por debajo del cual no se produce puesta en suspensión desde la capa cohesiva por fricción de fondo y ${\textstyle {K}_{k}}$ en [kg/s/N] es un coeficiente a definir por el usuario (flow-driven detachability coefficient). Nótese que el valor de ${\textstyle {\tau }_{crc}}$ es el mismo para todas las clases de sedimento, ya que mientras que la tensión de fondo no supere dicho valor no se producirá la erosión de las partículas de dicha capa. Una vez que se supera la tensión crítica de la capa cohesiva ( ${\textstyle \tau >{\tau }_{crc}}$), la contribución de cada clase de sedimento puede ser diferente, por lo que el coeficiente ${\textstyle {K}_{k}}$ puede ser diferente para cada clase.

Respecto al flujo de partículas desde la capa no-cohesiva a la columna de agua, su cálculo está ligado a los términos ${\textstyle {D}_{fdre}}$ y ${\textstyle {D}_{dep}}$, siendo el flujo neto de sedimento entre la capa de suelo erosionada y la columna de agua igual a la diferencia entre dichos términos. En Iber se contemplan dos formulaciones diferentes, siendo el usuario el que debe elegir cuál utilizar: Hairsine y Rose o van Rijn.

• Hairsine y Rose

${D}_{fdre,k}=\left\{{\begin{matrix}{\frac {{\rho }_{s}{r}_{f}{F}_{k}}{\left({\rho }_{s}-\rho \right)g}}\left({\frac {{\Omega -\Omega }_{cr}}{h}}\right)\epsilon \,{f}_{k}&si&\,{\Omega >\Omega }_{cr}\\0&si&{\Omega \leq \Omega }_{cr}\end{matrix}}\right.$

(80)

donde ${\textstyle \Omega }$ en [W/m²] es la potencia disipada por fricción de fondo y por unidad de superficie, ${\textstyle {\Omega }_{cr}}$ en [W/m²] es un valor crítico de ${\textstyle \Omega }$ , por debajo del cual no se produce puesta en suspensión desde la capa no-cohesiva por fricción de fondo, ${\textstyle {F}_{k}}$ es un coeficiente adimensional que expresa el tanto por uno de potencia disipada (en exceso sobre ${\textstyle {\Omega }_{cr}}$) que contribuye a la puesta en suspensión de sedimento desde la capa erosionada, ${\textstyle {\rho }_{s}}$ es la densidad de las partículas de sedimento (en torno a 2400 kg/m³, generalmente), ${\textstyle \rho }$ el la densidad del agua (1000 kg/m³), ${\textstyle h}$ [m] es el calado, ${\textstyle g}$ es la aceleración de la gravedad (9,81 m/s²), y ${\textstyle {r}_{f}}$ es el factor de Rouse que determina si el perfil de concentración en profundidad es más o menos homogéneo. Para flujos de escorrentía con calados muy bajos, se puede asumir un factor de Rouse igual a 1 ( ${\textstyle {r}_{f}=}$$1$).

En caso de utilizar el modelo de Hairsine y Rose, el término de deposición se calcula como:

${D}_{dep,k}={r}_{f}{w}_{e,k}{C}_{k}\quad \quad {w}_{e,k}=\mathrm {max} \,\left({w}_{s,k}-{u}_{\ast }\,,\,0\right)\,$

(81)

siendo ${\textstyle {w}_{s,k}}$ en [m/s] la velocidad de sedimentación de las partículas en suspensión y ${\textstyle {u}_{\ast }}$ la velocidad de fricción de fondo. La velocidad efectiva de sedimentación ${\textstyle {w}_{e,k}}$ tiene en cuenta el efecto de la turbulencia generada por fricción de fondo, que tiende a mantener en flotación a las partículas de sedimento. La velocidad de sedimentación de las partículas de sedimento se obtiene de su densidad y diámetro mediante la formulación de van Rijn [37].

El factor de Rouse se toma igual a 1 ( ${\textstyle {r}_{f}=1}$).

• van Rijn

${D}_{fdre,k}={\rho }_{s}\,{C}_{k}^{\ast }\,{w}_{s,k}\,\epsilon \,{f}_{k}$

(82)

donde ${\textstyle {\rho }_{s}}$ es la densidad de las partículas de sedimento, ${\textstyle {w}_{s,k}}$ en [m/s] es la velocidad de sedimentación de las partículas en suspensión calculada a partir de su diámetro y densidad mediante la formulación de van Rijn, y ${\textstyle {C}_{k}^{\ast }}$ es una concentración de equilibrio, la cual se puede calcular mediante alguna de las formulaciones siguientes:

• van Rijn [38]

${C}_{k}^{\ast }={\frac {0,005\,{T}_{k}^{1,5}}{{D}_{\ast ,k}^{0,3}}}\quad \quad {T}_{k}=$${\frac {{\tau }_{k}^{\ast }-{\tau }_{c,k}^{\ast }\,}{{\tau }_{c,k}^{\ast }}}\quad \,\quad {\tau }_{c,k}^{\ast }=$${\frac {{w}_{s,k}^{2}}{Rg{D}_{k}}}$

(83)

• Smith [39]

${C}_{k}^{\ast }={\frac {0,00156\,{T}_{k}}{1+0,0024{T}_{k}}}\quad \quad {T}_{k}={\frac {{\tau }_{k}^{\ast }-{\tau }_{c,k}^{\ast }\,}{{\tau }_{c,k}^{\ast }}}\,\,\quad \,\,{\tau }_{c,k}^{\ast }=$${\frac {{w}_{s,k}^{2}}{Rg{D}_{k}}}$

(84)

• García y Parker [40]

${C}_{k}^{\ast }={\frac {1,3\cdot {10}^{-7}\,{Z}_{k}^{5}}{1+4,33\cdot {10}^{-7}\,{Z}_{k}^{5}}}\quad \quad {Z}_{k}=$${\frac {{u}_{\ast }}{{w}_{s,k}}}{D}_{\ast ,k}^{0,9}$

(85)

• Zyserman y Fredsoe [41]

${C}_{k}^{\ast }={\frac {0,331\,{\left({\tau }_{k}^{\ast }-0,045\right)}^{1,75}}{1+0,72\,{\left({\tau }_{k}^{\ast }-0,045\right)}^{1,75}}}$

(86)

En caso de utilizar el modelo de van Rijn, el término de deposición se calcula como:

${D}_{dep,k}={r}_{f}{w}_{s,k}{C}_{k}$

(87)

siendo ${\textstyle {w}_{s,k}}$ en [m/s] la velocidad de sedimentación de las partículas en suspensión calculada mediante la formulación de van Rijn [37]. El factor de Rouse se toma igual a 1 ( ${\textstyle {r}_{f}=}$$1$).

5.1.5 Transporte por carga de fondo

El transporte por carga de fondo representa el movimiento de partículas de sedimento que se arrastran por la superficie del terreno como consecuencia de la fuerza de fricción ejercida por el agua sobre ellas. No son partículas que se transporten en suspensión en el seno de la columna de agua. No debe confundirse el transporte por carga de fondo con la puesta en suspensión de partículas por fricción de fondo, representado por los términos ${\textstyle {D}_{fde}}$ y ${\textstyle {D}_{fdre}}$ definidos anteriormente. En efecto, la fricción de fondo puede poner en suspensión ciertas partículas (las más finas y menos pesadas) y posteriormente estas partículas transportarse flotando en el seno de la columna de agua. En este caso estaríamos hablando de transporte en suspensión, y se calcularía con las formulaciones detalladas en el apartado 5.1.4. Pero cuando la fricción de fondo no tiene suficiente energía para poner en suspensión las partículas más pesadas, éstas se arrastrarán por la superficie del terreno. En este caso el transporte sería por carga de fondo.

El transporte por carga de fondo no suele considerarse en los modelos de erosión por lluvia-escorrentía, ya que generalmente se asume que este tipo de fenómenos producen básicamente el transporte de sedimento muy fino y ligero, que por lo tanto se transporta en suspensión. La mayor parte del suelo erosionado por lluvia-escorrentía son limos o arcillas. En el caso de que exista cierto transporte de fondo de las partículas más gruesas (arenas), la masa de suelo erosionada de estas fracciones de suelo suele ser mucho menor en magnitud que la erosionada por transporte en suspensión.

De todas formas, en el modelo de erosión implementado en Iber, sí que se ha considerado la componente de transporte de fondo, la cual puede ser activada o desactivada por el usuario. De esta manera se consigue una mayor versatilidad en el rango de aplicaciones del modelo, así como su posible aplicación a estudios de investigación en los que sea necesario cuantificar el transporte de fondo.

Condiciones de transporte de fondo de equilibrio y no-equilibrio

Para calcular el transporte de partículas de sedimento por carga de fondo se contemplan dos opciones: una formulación de equilibrio y una formulación de no-equilibrio. El usuario debe elegir qué formulación utilizar, siendo sensiblemente más sencilla la formulación de equilibrio.

• Formulación de equilibrio

En el caso de utilizar la formulación de equilibrio, se asume que las formulaciones de capacidad de transporte de fondo expuestas en el apartado siguiente son válidas para condiciones no-estacionarias (de no-equilibrio). En este caso, se asume que en todo momento el caudal sólido de fondo real es igual a la capacidad de transporte de fondo, simplificándose la resolución del modelo.

${q}_{b,k}={q}_{pot,k}$

(88)

donde ${\textstyle {q}_{b,k}}$ es el transporte de fondo de las partículas de la clase ${\textstyle k}$, y ${\textstyle {q}_{pot,k}}$ es la capacidad de transporte de fondo de la clase de sedimento ${\textstyle k}$.

• Formulación de no-equilibrio

En el caso de utilizar la formulación de no-equilibrio, el caudal sólido de fondo se calcula para cada clase de sedimento ( ${\textstyle {q}_{b,k}}$) mediante la siguiente ecuación:

$-{\frac {\partial }{\partial x}}\left({q}_{bx,k}\right)-{\frac {\partial }{\partial y}}\left({q}_{by,k}\right)+$${\frac {{q}_{b,k}-{q}_{pot,k}}{{L}_{s}}}=0$

(89)

siendo,

${L}_{s}=3{D}_{50}{D}_{\ast }^{0,6}{T}^{0,9}\quad \quad {D}_{\ast }={D}_{50}{\left({\frac {gR}{{\nu }^{2}}}\right)}^{\frac {1}{3}}\quad \,\,T=$${\frac {{\tau }^{\ast }-{\tau }_{c}^{\ast }}{{\tau }_{c}^{\ast }}}$

(90)

donde ${\textstyle {q}_{bx},{q}_{by}}$ son las proyecciones en los ejes X e Y del caudal sólido por transporte de fondo, y ${\textstyle {L}_{s}}$ es la longitud de adaptación a las condiciones de equilibrio. Un valor pequeño de ${\textstyle {L}_{s}}$ implica que el transporte de fondo se adapta muy rápidamente a los cambios en las condiciones de flujo, en cuyo caso ${\textstyle {q}_{b,k}=}$${q}_{pot,k}$ y la capacidad de transporte de fondo depende únicamente de las condiciones de flujo locales. Al contrario, un valor grande de ${\textstyle {L}_{s}}$ implica que la adaptación de la carga de fondo a las condiciones de flujo es relativamente lenta, evitándose transiciones bruscas en la capacidad de transporte de fondo de un punto a otro del dominio. Como se ha comentado, la utilización de la formulación de no-equilibrio complica de forma considerable el cálculo del transporte por carga de fondo.

Capacidad de transporte de fondo. Formulaciones

Tanto si se utiliza la formulación de equilibrio como la de no-equilibrio, es necesario calcular la capacidad de transporte de fondo en condiciones de equilibrio ( ${\textstyle {q}_{pot,k}}$) para cada clase de sedimento. Para ello, se implementan las siguientes formulaciones, siendo el usuario el que debe elegir cuál utilizar:

• Genérica

${q}_{pot,k}^{\ast }=A\,{\left({\tau }_{k}^{\ast }-{\tau }_{c,k}^{\ast }\right)}^{B}$

(91)

donde ${\textstyle A}$ y ${\textstyle B}$ son parámetros impuestos por el usuario y ${\textstyle {\tau }_{c,k}^{\ast }}$ es la tensión crítica adimensional, que puede ser calculada mediante el ábaco de Shields o la formulación de Parker et al. [42].

• Meyer-Peter y Müller [43]

${q}_{pot,k}^{\ast }=8\,{\left({\tau }_{k}^{\ast }-{\tau }_{c,k}^{\ast }\right)}^{1,5}\quad \quad {\tau }_{c,k}^{\ast }=$$0,047$

(92)

donde ${\textstyle {\tau }_{c,k}^{\ast }=0,047}$ es la tensión crítica adimensional.

• Wong y Parker [44]

${q}_{pot,k}^{\ast }=3,97\,{\left({\tau }_{k}^{\ast }-{\tau }_{c,k}^{\ast }\right)}^{1,5}\quad \quad {\tau }_{c,k}^{\ast }=$$0,0495$

(93)

donde ${\textstyle {\tau }_{c,k}^{\ast }=0,0495}$ es la tensión crítica adimensional.

• Einstein y Brown [45]

${\begin{matrix}{\tau }_{k}^{\ast }<0,18&{q}_{pot,k}^{\ast }=2,15\mathrm {exp} \,\left({\frac {-0,391}{{\tau }_{k}^{\ast }}}\right)\\0,18<{\tau }_{k}^{\ast }<0,52&{q}_{pot,k}^{\ast }=40\,{\left({\tau }_{k}^{\ast }\right)}^{3}\\0,52<{\tau }_{k}^{\ast }&{q}_{pot,k}^{\ast }=15\,{\left({\tau }_{k}^{\ast }\right)}^{1,5}\end{matrix}}$

(94)

En la formulación de Einstein-Brown la tensión crítica adimensional es igual a cero ( ${\textstyle {\tau }_{c,k}^{\ast }=}$$0)$.

• van Rijn [37]

${\begin{matrix}{T}_{k}<0,3&{q}_{pot,k}^{\ast }=0,053{\frac {{\left({T}_{k}\right)}^{2,1}}{{D}_{\ast ,k}^{0,3}}}&\\&&siendo\,{T}_{k}={\frac {{\tau }_{k}^{\ast }-{\tau }_{c,k}^{\ast }\,}{{\tau }_{c,k}^{\ast }}}\\{T}_{k}>0,3&{q}_{pot,k}^{\ast }=0,100{\frac {{\left({T}_{k}\right)}^{1,5}}{{D}_{\ast ,k}^{0,3}}}&\end{matrix}}$

(95)

donde ${\textstyle {\tau }_{c,k}^{\ast }}$ es la tensión crítica adimensional, a definir por el usuario según alguna de las formulaciones expuestas al final de esta sección.

• Engelund y Hansen [45]

${q}_{pot,k}^{\ast }={\frac {0,1\,{\left({\tau }_{k}^{\ast }\right)}^{2,5}}{\frac {2\,{u}_{\ast }^{2}\,{D}_{k}}{h\,{V}^{2}}}}$

(96)

En la formulación de Engelund-Hansen la tensión crítica adimensional es igual a cero ( ${\textstyle {\tau }_{c,k}^{\ast }=}$$0)$.

• Ashida y Michiue [45]

${q}_{pot,k}^{\ast }=17\left({\tau }_{k}^{\ast }-{\tau }_{c,k}^{\ast }\right)\left({\sqrt {{\tau }_{k}^{\ast }}}-\right.$$\left.{\sqrt {{\tau }_{c,k}^{\ast }}}\right)$

(97)

donde ${\textstyle {\tau }_{c,k}^{\ast }}$ es la tensión crítica adimensional, a definir por el usuario según alguna de las formulaciones expuestas al final de esta sección.

Tensión crítica para transporte de fondo

Algunas de las formulaciones anteriores (Meyer-Peter y Müller, Wong y Parker, Einstein y Brown, Engelund y Hansen) tienen un valor de la tensión crítica de fondo bien definido en la propia formulación, mientras que en otras (Genérica, van Rijn, Ashida y Michiue) este valor debe ser definido por el usuario (Tabla 4). En este último caso, el usuario puede elegir entre imponer un valor constante de la tensión crítica adimensional, calcularla mediante el ábaco de Shields o mediante la formulación de Parker et al. [42].

A estos efectos, el ábaco de Shields se implementa según la parametrización propuesta por Hoffmans y Verheij [46]:

${\begin{matrix}0<{D}_{\ast }<4&{\tau }_{c}^{\ast }={D}_{\ast }^{-1}\\4<{D}_{\ast }<10&{\tau }_{c}^{\ast }={D}_{\ast }^{-1}\\10<{D}_{\ast }<20&{\tau }_{c}^{\ast }={0,04\,D}_{\ast }^{-0,10}\end{matrix}}$ ${\begin{matrix}20<{D}_{\ast }<150&{\tau }_{c}^{\ast }={0,013\,D}_{\ast }^{0,29}\\{D}_{\ast }>150&{\tau }_{c}^{\ast }=0,055\end{matrix}}$

(98)

Alternativamente a la anterior parametrización del ábaco de Shields, puede utilizarse la siguiente formulación propuesta por Parker et al. para el cálculo de la tensión crítica, la cual proporciona valores de umbral del movimiento sensiblemente inferiores al ábaco de Shields original:

${\tau }_{c}^{\ast }=0,5\,\left(0,22\,{R}_{p}^{-0,6}+0,06\cdot {10}^{\left(-7,7\,{R}_{p}^{-0,6}\right)}\right)\quad \,\,{R}_{p}=$${\frac {\sqrt {Rg{D}^{3}}}{\nu }}$

(99)

Tabla 4. Tensión crítica adimensional para diferentes formulaciones de capacidad de transporte de fondo.

Formulación	Tensión crítica adimensional
Genérica	Constante, Shields o Parker et al. [42]
Meter-Peter y Müller	0,047
Wong y Parker	0,0495
Einstein y Brown	0,0
van Rijn	Constante, Shields o Parker et al. [42]
Engelund y Hansen	0,0
Ashida y Michiue	Constante, Shields o Parker et al. [42]

Factor de escondimiento (hiding factor)

En sedimentos con granulometría distribuida (multiclase), la tensión crítica adimensional ( ${\textstyle {\tau }_{c}^{\ast }}$) se calcula para el diámetro medio del sedimento. A partir de dicho valor, se calcula la tensión crítica para cada clase de partículas como:

${\tau }_{c,k}^{\ast }={\tau }_{c}^{\ast }{\left({\frac {{D}_{50}}{{D}_{k}}}\right)}^{\gamma }$

(100)

siendo ${\textstyle \gamma }$ un factor de ocultamiento (hiding factor) con un valor comprendido entre 0 y 1. Dicho factor representa la interacción entre partículas de diferente diámetro, de forma que las partículas gruesas protegen a las finas (aumentando su umbral de movimiento) y las finas exponen a las gruesas (facilitando su movilidad).

Teniendo en cuenta la definición de la tensión adimensional, la ecuación anterior se puede escribir como:

${\tau }_{c,k}={\tau }_{c}{\left({\frac {{D}_{k}}{{D}_{50}}}\right)}^{1-\gamma }$

(101)

De forma que si ${\textstyle \gamma =0}$ la interacción entre partículas es mínima (nula) y si ${\textstyle \gamma =}$$1$ la interacción entre partículas es máxima:

${\begin{matrix}Si\,\gamma =0&{\tau }_{c,k}^{\ast }={\tau }_{c}^{\ast }&{\tau }_{c,k}={\tau }_{c}{\frac {{D}_{k}}{{D}_{50}}}\\Si\,\gamma =1&{\tau }_{c,k}^{\ast }={\tau }_{c}^{\ast }{\frac {{D}_{50}}{{D}_{k}}}&{\tau }_{c,k}={\tau }_{c}\end{matrix}}$

(102)

Corrección del transporte de fondo por formas de fondo

La tensión de fondo total en el lecho de un río está generada por la rugosidad de grano del sedimento (la cual es proporcional al diámetro del sedimento) y por las formas de fondo (rizos, dunas o antidunas). Únicamente la tensión por grano contribuye al movimiento de sedimentos por carga de fondo. Si bien este fenómeno es principalmente importante en ríos, se ha incluido como opción en el módulo de erosión, de forma que el usuario puede activarlo o desactivarlo a su voluntad.

Si se activa esta opción, previamente al cálculo del caudal sólido de fondo se estima la tensión de fondo debida a grano utilizando para ello la partición de tensiones de Einstein, en la cual se calcula la tensión de grano a partir de la tensión total como:

${\tau }_{s}=\tau {\left({\frac {{n}_{s}}{n}}\right)}^{1,5}\quad \,\,{n}_{s}={\frac {{\left({K}_{s\,\left(m\right)}\right)}^{\frac {1}{6}}}{25}}\quad \quad \,{K}_{s}=$$\alpha \,{D}_{50}\quad \alpha \approx 2-3$

(103)

siendo ${\textstyle n}$ el coeficiente de Manning total (utilizado en el cálculo de la hidrodinámica e introducido por el usuario), ${\textstyle {n}_{s}}$ el coeficiente de Manning equivalente debido a grano, ${\textstyle {D}_{50}}$ el diámetro medio del sedimento expresado en metros, ${\textstyle {K}_{s}}$ la altura de rugosidad de grano expresada en metros (calculada a partir del diámetro del sedimento), ${\textstyle \tau }$ la tensión total de fondo, ${\textstyle {\tau }_{s}}$ la tensión de fondo debida a grano y ${\textstyle \alpha }$ un factor que relaciona la altura de rugosidad con el diámetro medio del sedimento (grain roughness factor). El factor ${\textstyle \alpha }$ lo debe introducir el usuario, si bien generalmente suelen utilizarse un valor en torno a 3.

Si se activa la corrección por formas de fondo, en el cálculo del transporte de fondo potencial ( ${\textstyle {q}_{pot,k}^{\ast }}$) mediante las ecuaciones presentadas anteriormente, se utiliza la tensión de grano ( ${\textstyle {\tau }_{s}}$) en lugar de la tensión total ( ${\textstyle \tau }$ ).

Corrección del transporte de fondo por pendiente del terreno

Las ecuaciones de transporte de fondo expuestas anteriormente son válidas únicamente para fondo horizontal. En el caso de que el fondo no sea horizontal es necesario introducir una corrección de pendiente que tenga en cuenta el efecto de la gravedad sobre la tensión crítica de movimiento y sobre el caudal sólido de fondo. El efecto de la gravedad puede ser positivo o negativo dependiendo de si esta actúa a favor o en contra de la velocidad del agua.

La corrección del caudal sólido de fondo por pendiente del terreno consiste, por un lado, en aumentar o minorar la tensión sobre el fondo del lecho (se mayora o minora dependiendo de si la fuerza de gravedad debida a la pendiente del fondo actúa a favor o en contra de la fuerza ejercida por el agua) y, por otro lado, en minorar la tensión crítica de movimiento (tensión máxima que puede soportar el sedimento sin que haya movimiento). De esta forma, se obtiene una tensión de fondo efectiva y una tensión crítica efectiva. En función de si la gravedad actúa a favor o en contra del movimiento de las partículas, la tensión efectiva sobre el fondo se mayora o se minora, mientras que la tensión crítica efectiva con fondo inclinado siempre es inferior a la tensión crítica con fondo horizontal.

La tensión crítica efectiva que es capaz de soportar el sedimento se evalúa como:

${\tau }_{c,e}={\tau }_{c}\mathrm {cos} \,\alpha$

(104)

siendo ${\textstyle \alpha }$ el ángulo que forma el terreno con la horizontal, y ${\textstyle {\tau }_{c}}$ la tensión crítica de arrastre con fondo horizontal. Por otro lado, la tensión efectiva sobre el fondo se calcula como:

${\tau }_{x,e}={\tau }_{x}-{\tau }_{c}{\frac {\mathrm {sin} \,{\alpha }_{x}}{\mathrm {tan} \,\varphi }}\quad \quad {\tau }_{y,e}=$${\tau }_{y}-{\tau }_{c}{\frac {\mathrm {sin} \,{\alpha }_{y}}{\mathrm {tan} \,\varphi }}\quad \,\,{\tau }_{e}=$${\sqrt {{\tau }_{x,e}^{2}+{\tau }_{y,e}^{2}}}$

(105)

siendo ${\textstyle {\tau }_{x},{\tau }_{y}}$ las tensiones de arrastre de fondo en las direcciones X e Y (tensión total o de grano en función de si se activa o no la corrección por formas de fondo), ${\textstyle {\tau }_{e}}$ la tensión de arrastre efectiva sobre el fondo, ${\textstyle {\alpha }_{x}}$, ${\textstyle {\alpha }_{y}}$ las proyecciones sobre los ejes X e Y del ángulo que forma el fondo con la horizontal y ${\textstyle \varphi }$ el ángulo de rozamiento interno del sedimento. Una vez calculadas la tensión crítica efectiva ( ${\textstyle {\tau }_{c,e}}$) y la tensión de arrastre efectiva ( ${\textstyle {\tau }_{e}}$), se utilizan estos valores en lugar de ${\textstyle {\tau }_{c}}$ y ${\textstyle \tau }$ en las ecuaciones de transporte correspondientes.

5.1.6 Evolución temporal de la masa erosionada

Para calcular la evolución temporal de la masa de sedimento en la capa no-cohesiva (erosionada) se resuelve la siguiente ecuación de balance de masa en la capa erosionada, para cada clase de sedimento considerada:

${\frac {\partial {M}_{s,k}}{\partial t}}=-{D}_{rdrd,k}-{D}_{fdre,k}+{D}_{dep,k}-{\rho }_{s}\left({\frac {\partial {q}_{bx,k}\,}{\partial x}}-\right.$$\left.{\frac {\partial {q}_{by,k}\,}{\partial y}}\right)$

(106)

siendo ${\textstyle {M}_{s,k}}$ la masa de sedimento de la clase ${\textstyle k}$ por unidad de superficie, expresada en [kg/m²]. Tanto ${\textstyle {M}_{s}}$, como los términos ${\textstyle {D}_{rdrd},\,{D}_{fdre}}$ y ${\textstyle {D}_{dep}}$, son variables en espacio y en tiempo.

En cada paso de tiempo se actualiza el valor de ${\textstyle {M}_{s,k}}$ para cada clase de sedimento, y a partir de estos valores la masa total de sedimento en la capa erosionada ( ${\textstyle {M}_{s}}$) y la nueva granulometría en dicha capa (factores ${\textstyle {f}_{k}}$), como:

${M}_{s}=\sum _{k}^{}{M}_{s,k}\quad \quad {f}_{k}={\frac {{M}_{s,k}}{{M}_{s}}}$

(107)

5.1.7 Evolución temporal de la topografía

La evolución temporal de la topografía debido a los procesos de erosión descritos anteriormente vendrá dada por la aplicación del principio de conservación de la masa para todas las clases de sedimento consideradas y en las dos capas de sedimento:

$\left(1-\phi \right){\frac {\partial {z}_{b}}{\partial t}}=\sum _{k}^{}\left(-\left({\frac {{D}_{rdd,k}+{D}_{rdrd,k}+{D}_{fde,k}+{D}_{fdre,k}-{D}_{dep,k}}{{\rho }_{s}\,}}\right)-\left({\frac {\partial {q}_{bx,k}\,}{\partial x}}-{\frac {\partial {q}_{by,k}\,}{\partial y}}\right)\right)$

(108)

siendo ${\textstyle \phi }$ la porosidad del sedimento.

En cada paso de tiempo, la nueva topografía obtenida a partir de esta ecuación se utiliza en las ecuaciones hidrodinámicas para asegurar un correcto acoplamiento entre los procesos de erosión y el flujo de agua.

5.1.8 Simplificaciones del modelo

Las ecuaciones y formulaciones implementadas en el módulo de erosión de suelos, descritas en las secciones anteriores, son relativamente complejas e incluyen un elevado número de parámetros. Esto puede complicar la aplicación del modelo en casos reales en los cuales sea difícil definir todos los parámetros y opciones del modelo. Aun disponiendo datos observados para la calibración del modelo, puede ser complicado ajustar todos sus parámetros debido a la alta no-linealidad de las formulaciones y a la posible equifinalidad entre parámetros [47]. Por último, puede ser interesante en algunas aplicaciones limitar la complejidad del modelo, por ejemplo, en su aplicación al lavado de sedimentos en entornos urbanos [48].

A continuación, se proponen algunas simplificaciones del modelo que pueden ser interesantes para su aplicación a casos reales.

Una única capa de sedimento y lavado de finos

La pérdida de suelo en cuencas rurales suele producirse principalmente por lavado de partículas finas del suelo. En este caso puede asumirse que el principal proceso de transporte de sólidos es por carga en suspensión, y despreciarse el transporte por carga de fondo ( ${\textstyle {q}_{b}=}$$0$), lo cual simplifica en gran medida el número de ecuaciones y parámetros del modelo.

Asimismo, puede simplificarse la estructura del suelo a una única capa (la capa erosionada), con el fin de evitar parametrizar dos capas de suelo con diferentes propiedades erosivas. En este caso los términos de producción de sedimento en suspensión ${\textstyle {D}_{rdd}}$ y ${\textstyle {D}_{fde}}$ son iguales a cero y no es necesario definir los parámetros correspondientes.

Si además:

• asumimos disponibilidad ilimitada de sedimento en el suelo, definiendo un valor inicial de masa de sedimento ( ${\textstyle {M}_{s}}$) elevado y un valor de masa crítica ( ${\textstyle {M}_{s,cr}}$) reducido, lo cual implica ${\textstyle \epsilon =}$$1$ siempre;
• despreciamos el factor protector de la lámina de escorrentía imponiendo un valor de ${\textstyle {z}_{m}}$ muy elevado, lo cual implica ${\textstyle {\sigma }_{h}=}$$1$;
• utilizamos la formulación de Harisine y Rose [49,50] para el lavado por escorrentía;

se obtiene un modelo con solamente 3 parámetros: ${\textstyle {\alpha }_{rd}}$, ${\textstyle {\Omega }_{cr}}$ y ${\textstyle F}$.

Podemos dar un paso más y despreciar la potencia crítica para la puesta en suspensión por escorrentía ( ${\textstyle {\Omega }_{cr}}$), argumentando que su valor suele ser varios órdenes de magnitud inferiores a ${\textstyle \Omega }$ , con lo cual nos quedaría un modelo de únicamente 2 parámetros: ${\textstyle {\alpha }_{rd}}$ y ${\textstyle F}$.

Este tipo de modelo fue utilizado en Cea et al. [47], proporcionando resultados satisfactorios para el cálculo del lavado de suelo en una parcela agrícola. Debe tenerse en cuenta que, en ciertas aplicaciones, es más interesante considerar la posible variabilidad espacial de estos coeficientes, que introducir parámetros adicionales que compliquen en exceso el modelo.

Una única capa de sedimento y lavado de finos por impacto de lluvia

Además de las simplificaciones realizadas anteriormente, en cuencas en las que no existe erosión en cárcavas (rill erosion o gully erosion) puede asumirse que el sedimento se pone en suspensión principalmente por acción de la lluvia, con lo cual podemos despreciar en primera aproximación el término ${\textstyle {D}_{fdre}}$.

De esta forma nos quedaría un modelo con ${\textstyle {\alpha }_{rd}}$ como único parámetro (que podría ser variable espacialmente).

Lavado de sólidos en zonas urbanas

Este tipo de modelos de erosión de suelos se puede aplicar, también, para la modelización del lavado de sólidos en zonas urbanas. Para este tipo de aplicación se asume una única capa de sedimento no-cohesivo que representaría los sólidos sobre el pavimento. El pavimento sería la capa inferior, que en este caso es no erosionable, por lo que no es necesario incluirla en el cálculo.

Bajo estas condiciones, los términos de producción desde la capa inferior son igual a cero ( ${\textstyle {D}_{rdd}=}$${D}_{fde}=0$). Si se utiliza la formulación de Harisine-Rose para el lavado por escorrentía se obtiene un modelo con 5 parámetros: ${\textstyle {\alpha }_{rd},\,{z}_{m},\,{M}_{s,cr},\,{\Omega }_{cr},\,F}$, similar al que se utilizó en Naves et al. [51]. En dicha publicación se pueden, además, encontrar rangos de variación realistas de dichos parámetros en un entorno urbano controlado en laboratorio.

5.2 Interfaz de usuario

A continuación, se describen los diferentes menús y ventanas desde las que se pueden activar los procesos descritos en las secciones anteriores y definir los correspondientes parámetros.

5.2.1 Formulaciones, opciones y parámetros generales del modelo

Para definir las opciones y parámetros generales del caso a modelar, debe accederse a la ventana de ‘Datos del problema’, desde el menú de Datos >> Datos del problema >> Erosión de suelos. Desde dicha ventana se puede activar el módulo de erosión de suelos y seleccionar los diferentes modelos y procesos a considerar en el cálculo, así como los resultados que se desean ver en la interfaz de post-proceso. Las principales opciones a las que se accede desde esta ventana son las siguientes (Fig. 33):

• Activar/desactivar el módulo de erosión de suelos y seleccionar la asignación manual o automática de las zonas de erosión.
• Estructura vertical del suelo o modelo de mezcla.
• Propiedades físicas del suelo y parámetros de transporte (para cada clase de sedimento si en el cálculo se considera granulometría distribuida).
• Activar/desactivar el transporte en suspensión, así como definir las formulaciones y parámetros correspondientes.
• Activar/desactivar el transporte de fondo, así como definir las formulaciones y parámetros correspondientes.
• Activar/desactivar diferentes correcciones a aplicar en las formulaciones utilizadas en el modelo (opcionales).
• Seleccionar las variables que se guardarán en los ficheros de resultados para poder verlas posteriormente en el post-proceso.

Fig. 33. Datos generales del módulo de erosión de suelos.

5.2.2 Asignación de zonas erosionables

Una vez definidas las formulaciones y datos generales que se utilizarán en el cálculo, deben definirse y asignarse las zonas del dominio susceptibles de ser erosionadas (zonas erosionables). Esto se hace desde el menú Datos >> Erosión de suelos >> Zonas erosionables.

A través de la ventana que se muestra en la Fig. 34a pueden crearse tantas zonas erosionables como se desee. Es recomendable darle un nombre a cada zona erosionable diferente al nombre que se genera por defecto. Para cada zona erosionable debe definirse la cantidad de material erosionable (kg/m²) disponible inicialmente para cada tipo de sedimento en la capa superior de suelo. Esta cantidad inicial de sedimento variará a lo largo del cálculo en función de los procesos de erosión y sedimentación que tengan lugar. Asimismo, se debe establecer si este material erosionable está situado sobre la cota definida por el MDT o por debajo de dicha cota. Por último, es necesario definir la concentración inicial de sólidos en suspensión en dicha zona (kg/m³). Lo más común en cálculos hidrológicos es que las zonas estén secas al comienzo del cálculo, en cuyo caso la concentración inicial de sólidos en suspensión será igual a cero.

Si un tipo de sedimento no está presente en dicha zona no es necesario incluirlo en la tabla.

Asignación manual

La asignación manual de zonas erosionables (ya deben estar creadas) se realiza mediante el menú Datos >> Erosión de suelos >> Asignación manual de zonas de erosión. Este menú da acceso a la ventana que se muestra en la Fig. 34b, en la cual se puede elegir una zona erosionable creada anteriormente y asignarla a diferentes regiones del modelo. Al igual que otras propiedades de Iber, la asignación puede realizarse en la geometría o en la malla. En general es preferible realizar la asignación sobre la geometría del modelo y posteriormente generar la malla de cálculo.

Fig. 34. (a) Definición de zonas erosionables. (b) Ventana de asignación manual de zonas erosionables.

Asignación automática

Para asignar las zonas erosionables de forma automática es necesario disponer de un archivo raster en formato ASCII Grid de ArcInfo (*.asc) que contenga la distribución espacial de las zonas en base a su identificador (ID). El identificador de cada zona puede conocerse cambiando el nombre de dicha zona (como se ha comentado anteriormente es recomendable asignarle un nombre a cada zona erosionable diferente al que se genera por defecto), apareciendo automáticamente a modo de prefijo del nombre asignado como puede verse en la Fig. 35a. En el ejemplo mostrado en la Fig. 35b se le ha dado el nombre Gravera, y al darle a Aceptar el programa la guarda como ID2 – Gravera (ya que existía previamente otra zona con el identificador ID1). Si se crea una nueva zona tendrá el prefijo ID3.

Fig. 35. (a) Ventana de cambio de nombre de una zona erosionable. (b) Ventana de asignación según el prefijo ID de la zona erosionable predefinida.

Sabiendo cual es el ID de cada zona y disponiendo del archivo ráster con su distribución espacial, accedemos a la asignación automática mediante el menú Datos >> Erosión de suelos >> Asignación automática de zonas de erosión, desde dónde se deberá seleccionar el archivo ráster. Como siempre en Iber, la asignación automática de propiedades se realiza directamente sobre la malla de cálculo, por lo que esta deberá estar creada con anterioridad. Asimismo, cada vez que se cree la malla de cálculo, deben realizarse todas las asignaciones automáticas de nuevo.

5.2.3 Condiciones de contorno

Una vez definidos los datos generales y asignado las zonas erosionables, deben asignarse las condiciones de contorno, tanto para el transporte en suspensión como para el transporte de fondo. Dichas condiciones deben asignarse únicamente en los contornos de entrada de caudal. En caso de no asignarse se considerará que no entra sedimento por dichos contornos.

Las condiciones de contorno se asignan desde el menú Datos >> Erosión de suelos >> Condiciones de contorno. Este menú da acceso a la ventana donde se puede elegir entre asignar las condiciones de contorno para el transporte en suspensión o de fondo.

Condiciones de contorno para el transporte en suspensión

Para asignar las condiciones de contorno de transporte en suspensión debe definirse la concentración de sólidos en suspensión (para cada clase de sedimento considerada) en los contornos de entrada de caudal. Para ello se define la concentración de sólidos totales en suspensión (kg/m³) como una serie temporal (Fig. 36a). En el caso de que se trabaje con varias clases de sedimentos, dicha concentración corresponde a la suma de todas las clases, por lo que también será necesario definir la fracción de cada clase de sedimento en dicho contorno, en la columna Fracción de material. La suma de los valores introducidos para las clases de sedimento en la columna Fracción de material debe ser igual a 1.

Fig. 36. Ventana de asignación de condiciones de contorno para transporte de sedimentos por: (a) carga en suspensión; (b) carga de fondo.

Condiciones de contorno para el transporte de fondo

Para el transporte de fondo se debe de asignar, en cada contorno de entrada de caudal, la evolución temporal del caudal sólido total que entra por dicho contorno (kg/s), y la fracción de dicho caudal que corresponde a cada tipo de sedimento (Fig. 36b). Al igual que en el caso de transporte en suspensión, la suma de las fracciones de sedimento (columna Fracción de material) debe ser igual a 1.

6 Módulo de transporte de sedimento de fondo para mezclas

6.1 Fundamentos teóricos

En este apartado se presenta el módulo de transporte de fondo para material no uniforme, es decir, cuando el sedimento transportado por un río no puede caracterizarse por un único diámetro. Es un hecho que un río que es alimentado con sedimento con distinto tamaño de grano, tiene la capacidad de clasificarlo. Esto es especialmente significativo en los ríos con lecho formado por una mezcla de gravas y arenas, pero también ocurre para otros tamaños, como pueden ser bolos, o gravas con distintos tamaños característicos. Esta clasificación acaba resultando en que, para muchos ríos, el fondo se encuentra formado por estratos con distintas granulometrías, con su interfaz sensiblemente horizontal, y una capa de sedimento grueso (coraza) en su parte superior.

Los aspectos relacionados con la variedad de diámetros de sedimento de un río tienen repercusiones importantes en distintas facetas de la ingeniería de ríos, como pueden ser las restauraciones fluviales, la preservación de hábitos, la propia evolución geomorfológica del río o en la interacción entre el material transportado por el fondo y las estructuras construidas en el mismo.

6.1.1 Concepto de la capa activa

En el transporte de sedimentos por arrastre de fondo para material distribuido, es decir, para mezclas de sedimento, es un fenómeno que se produce en un espesor determinado del fondo del río [52]. Es común explicar este fenómeno a través del concepto de la capa activa. Aunque existen distintas interpretaciones de la capa activa, la más utilizada es que se puede considerar como capa activa el estrato superior del fondo de un río que se mueve y que es donde la clasificación del material realmente tiene lugar (Fig. 37a). De esta manera, se suele llamar capa activa a la capa superior del fondo del río. Desde el punto de vista de la modelización numérica, su espesor se especifica mediante un parámetro, N, multiplicado por D₅₀ o D₉₀ del tamaño más grande utilizado en la simulación. La capa se suele considerar constante durante una simulación, aunque la composición de la capa activa en cada sección transversal es una función del tiempo e incorpora los efectos acumulativos de la erosión selectiva y la deposición [52,53].

Fig. 37. Teoría de la capa activa: (a) esquema de la conceptualización de la estructura del fondo de un río bajo la idea de transporte por capa activa (fuente: Spasojevic y Holly [54]); (b) división de la estructura del fondo de un río en capa activa y estratos.

6.1.2 Discretización del sustrato

Para poder calcular el transporte de fondo utilizando el concepto de capa activa del terreno erosionable, el sustrato se divide en una serie de estratos (Fig. 37b). Cada capa, en cada instante y en cada punto, tendrá su granulometría. La definición de la granulometría se define en base a una serie de clases o diámetros característicos del sedimento. Es decir, en cada punto de una capa o estrato, habrá un porcentaje (en peso) de cada diámetro característico.

6.1.3 Ecuaciones

El transporte sólido total (caudal solido) en cada punto ${\textstyle {q}_{b}\,}$ es la suma del caudal sólido para cada una de las clases o diámetros ${\textstyle {q}_{ks}}$ . La proporción o fracción de caudal de cada clase se puede expresar como:

$\,{\mathit {\boldsymbol {q}}}_{b}=\sum _{k=1}^{n}{q}_{{k}_{s}}\,\,$

(109)

donde ${\textstyle {q}_{{k}_{s}}}$ es la capacidad de transporte calculada con alguna fórmula de transporte sólido.

La ecuación de Exner, que representa la continuidad del caudal sólido en el fondo de un río, para el caso de mezclas, se escribe como:

$\left(1-p\right){\frac {\partial {z}_{b}}{\partial t}}+\sum _{{k}_{s}=1}^{n{^{\circ }}classes}{\frac {\partial {q}_{{k}_{s},x}}{\partial x}}+{\frac {\partial {q}_{{k}_{s},y}}{\partial y}}=0\,$

(110)

La ecuación que representa el balance de masa total en un punto de un río, si se realiza por cada una de las clases o diámetros, es:

$\left(1-p\right){\frac {\partial \left({\beta }_{{k}_{s}}{h}_{a}\right)}{\partial t}}+$${\frac {\partial {q}_{{k}_{s},x}}{\partial x}}+{\frac {\partial {q}_{{k}_{s},y}}{\partial y}}=$${\boldsymbol {F}}_{{k}_{s}}$

(111)

En las expresiones anteriores ${\textstyle {q}_{{k}_{s}}}$, es el caudal sólido unitario de la clase ${\textstyle {k}_{s}}$, ${\textstyle {\beta }_{{k}_{s}}}$ la fracción de dicha clase ${\textstyle {k}_{s}}$ en la capa active y ${\textstyle {F}_{{k}_{s}}}$ la fuente de sedimentos de la clase ${\textstyle {k}_{s}}$ hacia la capa activa procedente del substrato inferior y ${\textstyle p}$ representa la porosidad.

El intercambio de material entre la capa activa (móvil) y la capa inferior (estacionaria) se puede determinar según:

${F}_{{k}_{s}}=-\left(1-p\right){\frac {\partial \left[{\beta }^{'}\left({z}_{b}-{h}_{a}\right)\right]}{\partial t}}$

(112)

siendo ${\textstyle {\beta }^{'}}$:

• ${\textstyle {\beta }_{s{k}_{s}}}$, si el flujo es ascendente, del substrato hacia la capa activa
• ${\textstyle {\beta }_{{k}_{s}}}$, si el flujo es descendente, desde la capa activa al sustrato

${\textstyle {\beta }_{s{k}_{s}}}$, es la fracción de la clase ${\textstyle {k}_{s}}$ en el estrato activo (estrato en el que se encuentra la capa activa en cada instante y punto). Debe cumplirse que:

$\sum _{{k}_{s}=1}^{n{^{\circ }}\,classes}{\beta }_{{k}_{s}}=1\,$

(113)

El caudal sólido se puede calcular con cualquiera de las fórmulas de transporte que incluye Iber (Meyer-Peter y Müller, Engelund y Fredsoe, Engelund y Hansen, van Rijn, Recking, y una formulación ad-hoc donde el usuario puede modificar los coeficientes). En la versión 3 se ha añadido la de Ashida y Michiue [55], por ser una de las comúnmente utilizadas para mezclas. Con ella, el caudal sólido adimensional para una clase ${\textstyle i}$ se define como:

${q}_{i}^{\ast }=17\left({\tau }_{{s}_{i}}^{\ast }-{\tau }_{{sc}_{i}}^{\ast }\right)\left({\sqrt {{\tau }_{{s}_{i}}^{\ast }}}-\right.$$\left.{\sqrt {{\tau }_{{sc}_{i}}^{\ast }}}\right)\,$

(114)

A partir de allí, el caudal sólido de la clase es:

${q}_{i}={q}_{i}^{\ast }\left({F}_{i}{\sqrt {\gamma g{D}_{i}}}{D}_{i}\right)\,$

(115)

Siendo, por tanto, el caudal sólido total:

$q=\sum _{}^{}{q}_{i}\,$

(116)

6.1.4 Escondimiento (hiding)

En el cálculo del caudal sólido para mezclas de sedimentos es importante considerar el hecho de que una partícula de un diámetro característico no se mueve de la misma manera cuando forma parte de una mezcla, como lo haría si el resto de partículas tuvieran el mismo tamaño. Este efecto se conoce como escondimiento, ampliamente estudiado por Ashida y Michiue [55], y más adelante por otros muchos autores.

Por consistencia, en Iber se ha implementado la ecuación de escondimiento de los primeros, pero también la interpretación de la misma según otros autores, incluso la posibilidad de modificar los parámetros a criterio del usuario. Según Ashida y Michiue, la relación entre la tensión tangencial crítica adimensional para una clase ${\textstyle i}$, ${\textstyle {\tau }_{{sc}_{i}}^{\ast }}$, se relaciona con la tensión tangencial crítica adimensional que se tendría con un diámetro característico ${\textstyle {D}_{g}^{}}$, ${\textstyle {\tau }_{{sc}_{g}}^{\ast }.\,{D}_{g}^{}\,}$ corresponde a la media geométrica de la distribución granulométrica del estrato según Egiazaroff [56]:

${\frac {{\tau }_{{sc}_{i}}^{\ast }}{{\tau }_{{sc}_{g}}^{\ast }}}={\frac {{D}_{i}^{}}{{D}_{g}^{}}}\,{\left[{\frac {\mathrm {ln} \,(19)}{\mathrm {ln} \,\left(19{\frac {{D}_{i}^{}}{{D}_{g}^{}}}\right)}}\right]}^{2}$

(117)

Según Ashida y Michiue [55], la expresión sería:

${\begin{matrix}Si\,{\frac {{D}_{i}^{}}{{D}_{g}^{}}}<0,4&{\frac {{\tau }_{{sc}_{i}}^{\ast }}{{\tau }_{{sc}_{g}}^{\ast }}}={0,843\left({\frac {{D}_{i}^{}}{{D}_{g}^{}}}\,\right)}^{-1}\\Si\,{\frac {{D}_{i}^{}}{{D}_{g}^{}}}\geq 0,4&{\frac {{\tau }_{{sc}_{i}}^{\ast }}{{\tau }_{{sc}_{g}}^{\ast }}}={\frac {{D}_{i}^{}}{{D}_{g}^{}}}\,{\left[{\frac {\mathrm {ln} \,(19)}{\mathrm {ln} \,\left(19{\frac {{D}_{i}^{}}{{D}_{g}^{}}}\right)}}\right]}^{2}\end{matrix}}$

(118)

Según Parker, Klingeman y McLean [57]:

${\frac {{\tau }_{{sc}_{i}}^{\ast }}{{\tau }_{{sc}_{g}}^{\ast }}}={\left({\frac {{D}_{i}^{}}{{D}_{g}^{}}}\,\right)}^{-0,982}$

(119)

Según Wilcock y Crowe [58]:

${\tau }_{{c}_{i}}^{\ast }={\tau }_{{c}_{g}}^{\ast }{\left({\frac {{D}_{i}}{{D}_{g}}}~\right)}^{-~\left({\frac {0,67}{1-{e}^{1,5-{\frac {{D}_{i}}{{D}_{g}}}}}}\right)}$

(120)

Asimismo, se ha implementado la posibilidad de considerar una fórmula ad-hoc, con una expresión de forma similar a la anterior, pero con coeficientes ajustables por el usuario ( ${\textstyle \alpha }$ y ${\textstyle \beta }$ ):

${\frac {{\tau }_{{sc}_{i}}^{\ast }}{{\tau }_{{sc}_{g}}^{\ast }}}=\alpha {\left({\frac {{D}_{i}^{}}{{D}_{g}^{}}}\,\right)}^{\beta }$

(121)

6.1.5 Implementación de las ecuaciones en el esquema numérico

La ecuación de Exner presentada anteriormente no es nada más que la ecuación de conservación de la masa para el sedimento arrastrado por el fondo del río. Mediante el método de los volúmenes finitos, presentado en el manual de referencia hidráulico de Iber, resulta que la ecuación para actualizar la cota del fondo del río ${\textstyle {z}_{i}^{n}}$ en cada instante se puede escribir como:

${z}_{i}^{n+1}={z}_{i}^{n}-{\frac {\Delta t}{\left(1-p\right){V}_{i}}}\sum _{{k}_{s}=1}^{classes}\sum _{w=1}^{~lados}\left({\overline {q}}_{{k}_{s},i,w}{\overline {n}}_{i,w}\right){l}_{i,w}$

(122)

donde ${\textstyle {\overline {q}}_{{k}_{s},i,w}}$, es el caudal sólido unitario para la clase ${\textstyle {k}_{s}}$, en el lado (o pared) w del elemento (o volumen finito) ${\textstyle i}$. Para la capa activa, la actualización del porcentaje de cada clase, que se puede considerar como una ecuación de Exner local. Se expresará entonces como:

$\left({h}_{a}{\beta }_{{k}_{s}i}^{n+1}\right)=\left({h}_{a}{\beta }_{{k}_{s}i}^{n}\right)-$${\frac {\Delta t}{\left(1-p\right){V}_{i}}}\sum _{w=1}^{lados}\left({\overline {q}}_{{k}_{s},i,w}{\overline {n}}_{i,w}\right)-{F}_{{k}_{s}i}^{n}\,$

(123)

La fuente de fondo, o intercambio entre el estrato activo y la capa activa, responde a:

${F}_{{k}_{s}i}^{n+1}=-\left(1-p\right)\left[{\frac {\left({\beta }_{i}^{n+1}\left({z}_{b}^{n+1}-{h}_{a}^{n+1}\right)\right)-\left({\beta }_{i}^{n}\left({z}_{b}^{n}-{h}_{a}^{n}\right)\right)}{\Delta t}}\right]\,$

(124)

${\textstyle \beta }$ , desde el punto de vista numérico, toma la forma:

${\begin{matrix}{\beta }_{{sk}_{s\,}\,}\,\,si\quad \,\,{\frac {\left({z}_{b}^{n+1}-{h}_{a}^{n+1}\right)-\left({z}_{b}^{n}-{h}_{a}^{n}\right)}{\Delta t}}>0\,\,\\{\beta }_{{k}_{s\,}\,}\,\,si\quad \,\,{\frac {\left({z}_{b}^{n+1}-{h}_{a}^{n+1}\right)-\left({z}_{b}^{n}-{h}_{a}^{n}\right)}{\Delta t}}<0\,\end{matrix}}$

(125)

Para la actualización de la fracción de cada clase en el estrato activo, las expresiones resultantes son:

${\begin{matrix}Si\,{\,\left({z}_{b}-{h}_{a}\right)}^{n+1}<{\left({z}_{b}-{h}_{a}\right)}^{n}&\rightarrow &{\beta }_{{sk}_{s}}^{n+1}={\beta }_{{sk}_{s}}^{n}\\Si\,\,{\left({z}_{b}-{h}_{a}\right)}^{n+1}>{\left({z}_{b}-{h}_{a}\right)}^{n}&\rightarrow &{\beta }_{{sk}_{s}}^{n+1}={\beta }_{{sk}_{s}}^{n}{\left({z}_{b}-{h}_{a}-{z}_{{s}_{1}}\right)}^{n}\end{matrix}}$

(126)

siendo,

${\beta }_{{sk}_{s}}^{n+1}={\frac {{F}_{{k}_{s}}^{n+1}\Delta t+{\alpha }_{{sk}_{s}}^{n}}{{\left({z}_{b}-{h}_{a}-{z}_{{s}_{1}}\right)}^{n+1}}}\,$

(127)

6.2 Interfaz de usuario

A continuación, se describen los diferentes menús y ventanas desde las que se pueden activar los procesos descritos en las secciones anteriores, así como definir los correspondientes parámetros. Este nuevo módulo se ha creado en base a tablas dinámicas y, por tanto, requiere de un proceso secuencial para su correcto funcionamiento. Para evitar problemas para con la implementación de las condiciones, se han implementado mensajes de ayuda para que el usuario sepa cómo debe proceder.

Fig. 38. Datos generales del módulo de transporte de sedimentos (mezclas).

6.2.1 Formulaciones, opciones y parámetros generales del modelo

Para definir las opciones y parámetros generales del caso a modelar, debe accederse a la ventana de ‘Datos del problema’, desde el menú de Datos >> Datos del problema >> Sedimentos. Desde dicha ventana se puede activar el módulo de transporte de sedimentos, ya sea por suspensión como por fondo. En el caso de mezclas de sedimentos, esta opción requiere la activación del transporte de fondo y seleccionar en Tipo de sedimento la opción Mezclas. A continuación, se requiere la introducción de diferentes parámetros, que también son comunes para granulometría uniforme. Las principales opciones a las que se accede desde esta ventana son las siguientes (Fig. 38):

• Activar/desactivar el módulo de transporte de sedimentos.
• Selección del tipo de sedimento.
• Propiedades del sedimento, como son la porosidad, el ángulo de fricción y la densidad relativa del sedimento.
• Asimismo, se ofrece la posibilidad de acondicionar el terreno (pendiente) según las propiedades del sedimento a través de la activación del Modelo de Avalancha

6.2.2 Definición del número de estratos, clases y espesor de la capa activa

La definición del número de estratos y clases del sedimento, así como el espesor de la capa activa, es primordial para la implementación de las propiedades del sedimento, ya que de estos valores depende la generación del resto de ventanas.

A esta ventana se accede desde el menú Datos >> Transporte de sedimentos >> Transporte de fondo (mezclas) >> Estratos&Clases. En la ventana que se muestra en la Fig. 39a deben definirse el número de estratos y clases del sedimento, así como el espesor de la capa activa (en metros). Estos parámetros son fijos para todo el modelo.

Fig. 39. (a) Definición del número de estratos, clases y espesor de la capa activa. (b) Definición del diámetro de cada clase.

6.2.3 Definición de las clases

La definición del número y diámetro de cada clase se realiza mediante el menú Datos >> Transporte de sedimentos >> Transporte de fondo (mezclas) >> Definición de clases. Este menú da acceso a la ventana que se muestra en la Fig. 39b, en la cual se debe introducir la relación de cada clase y el diámetro de partícula. El valor del diámetro de cada clase es constante, pero no así el porcentaje de cada uno de ellos en cada capa, ya que variará a lo largo del proceso de cálculo.

6.2.4 Propiedades de las mezclas

A continuación, se deben definir las propiedades del sedimento, es decir, introducir el espesor de cada capa y la relación de porcentajes que definen la granulometría en base al número de clases que se ha predefinido. Esto se realiza desde el menú Datos >> Transporte de sedimentos >> Transporte de fondo (mezclas) >> Propiedades de las mezclas (Fig. 40a). En esta ventana se deben definir las propiedades de cada capa, indicando el identificador de la capa, el espesor de la capa (en metros) y el porcentaje (%) de cada clase, que está asociado al diámetro de la partícula previamente definido.

Fig. 40. (a) Definición de las propiedades de cada estructura del suelo. (b) Asignación de la estructura del suelo predefinida.

En esta ventana también se debe especificar la distancia entre la capa superior del estrato 1 y el fondo, y que determina la posición de los estratos. Tiene que ser un número mayor que cero. Si es cero, es que el fondo coincide con este estrato. Si es mayor que cero, quiere decir que en el momento inicial el fondo está erosionado, o que puede haber sedimentación, por lo que pueden haber estratos por encima de la cota de fondo inicial.

Se pueden crear tantas “estructuras de suelo” de sedimentos como se desee, siendo el número de estratos y de clases fijo para todo el modelo.

6.2.5 Asignación de las mezclas

Una vez definidas las “mezclas”, es posible asignarlas mediante el menú Datos >> Transporte de sedimentos >> Transporte de fondo (mezclas) >> Asignación de mezclas. Esto se puede hacer tanto en la geometría como en la malla de cálculo (Fig. 40b). Todas aquellas “mezclas” que se hayan creado en la ventana anterior, se mostrarán en una lista para poder asignarlas al modelo.

6.2.6 Condiciones de contorno

Una vez definidos los datos generales y definidas/asignadas las propiedades del sedimento, deben asignarse las condiciones de contorno. Dichas condiciones deben asignarse únicamente en los contornos de entrada de caudal. En caso de no asignarse, se considerará que no entra sedimento por dichos contornos.

Las condiciones de contorno se asignan desde el menú Datos >> Transporte de sedimentos >> Transporte de fondo (mezclas) >> Condiciones de contorno (mezclas). Este menú da acceso a la ventana que se muestra en la Fig. 41, desde donde se puede elegir entre asignar las condiciones de contorno en función del tipo de entrada: no-sedimentos (o agua clara), por capacidad de arrastre y como caudal mediante un sedimentograma.

No –Agua clara–

Esta opción introducirá caudal de agua, pero no de sedimentos. Es decir, no habrá aporte de sedimentos por los contornos de entrada de caudal.

Capacidad de arrastre

Con esta opción, el modelo considera la entrada de sedimentos en los contornos de entrada que se indique, siendo el tipo de sedimento (mezcla) idéntico al que hay en el elemento interior de dicho contorno. La cantidad de sedimento entrante dependerá, por tanto, del tipo de sedimento impuesto en la zona inmediatamente adyacente a la condición de contorno, así como a las propiedades del mismo.

Sedimentograma

Finalmente, se proporciona la posibilidad de imponer una condición de contorno ad hoc. En este caso se debe definir el caudal de entrada de material sólido, es decir, un sedimentograma. El caudal sólido de entrada total será el especificado. La granulometría será la asignada en el elemento junto a la entrada.

Fig. 41. Ventana de definición e implementación de las condiciones de contorno para mezclas.

6.2.7 Resultados

Para el caso de sedimento con granulometría no-uniforme, los resultados específicos aparecen dentro del análisis con el nombre de Sedimentos. Con las herramientas habituales de consulta de resultados se pueden analizar los resultados globales o por clases (Fig. 42). La particularidad de estos resultados, accesibles desde el análisis Sedimentos, es que se han agrupado por estratos (incluyendo un análisis propio para la capa activa). En cada estrato se puede consultar la proporción de cada clase de sedimento, hallándose cada uno de ellos a su correspondiente profundidad.

Fig. 42. Organización de los resultados para mezclas.

Asimismo, se ha creado una opción nueva, accesible a través del menú Herramientas Iber >> Transporte de Sedimentos >> Granulometría, que permite seleccionar unos elementos y ver la granulometría en ellos, tanto para la capa activa como para cada uno de los estratos existentes (Fig. 43).

Fig. 43. Distribución granulométrica de una capa del sustrato en un modelo simulado con granulometría no uniforme.

Referencias

[1] E. Bladé Castellet, L. Cea, G. Corestein, Numerical modelling of river inundations, Ing. Del Agua. 18 (2014) 68. doi: 10.4995/ia.2014.3144.

[2] E. Bladé, L. Cea, G. Corestein, E. Escolano, J. Puertas, E. Vázquez-Cendón, J. Dolz, A. Coll, Iber: herramienta de simulación numérica del flujo en ríos, Rev. Int. Métodos Numéricos Para Cálculo y Diseño En Ing. 30 (2014) 1–10. doi: 10.1016/j.rimni.2012.07.004.

[3] M. Sanz-Ramos, G. Olivares Cerpa, E. Bladé, Metodología para el análisis de rotura de presas con aterramiento mediante simulación con fondo móvil, Ribagua. 6 (2019) 138–147. doi: 10.1080/23863781.2019.1705198.

[4] J. Anta Álvarez, M. Bermúdez, L. Cea, J. Suárez, P. Ures, J. Puertas, Modelización de los impactos por DSU en el río Miño (Lugo), Ing. Del Agua. 19 (2015) 105. doi: 10.4995/ia.2015.3648.

[5] L. Cea, M. Bermúdez, J. Puertas, E. Bladé, G. Corestein, E. Escolano, A. Conde, B. Bockelmann-Evans, R. Ahmadian, IberWQ: new simulation tool for 2D water quality modelling in rivers and shallow estuaries, J. Hydroinformatics. 18 (2016) 816–830. doi: 10.2166/hydro.2016.235.

[6] M. Sanz-Ramos, E. Bladé Castellet, A. Palau Ibars, D. Vericat Querol, A. Ramos-Fuertes, IberHABITAT: evaluación de la Idoneidad del Hábitat Físico y del Hábitat Potencial Útil para peces. Aplicación en el río Eume, Ribagua. 6 (2019) 158–167. doi: 10.1080/23863781.2019.1664273.

[7] M. Sanz-Ramos, E. Bladé, L. Cea, Iberaula - Dissemination, (2022). https://iberaula.es/56/iber-community/dissemination (accessed July 5, 2022).

[8] A. Coll, R. Ribó, M. Pasenau, E. Escolano, J.S. Perez, A. Melendo, A. Monros, J. Gárate, GiD v.14 User Manual, (2018). https://www.gidsimulation.com/

[9] A. Coll, M. Pasenau, E. Escolano, J.S. Perez, A. Melendo, A. Monros, J. Gárate, www.gidhome.com, (2018). www.gidhome.com.

[10] A. Coll, R. Ribó, M. Pasenau, E. Escolano, J.S. Perez, A. Melendo, A. Monros, J. Gárate, GiD v.14 Reference Manual, (2018). https://www.gidsimulation.com/

[11] J.L. Aragón-Hernández, E. Bladé, Modelación numérica de flujo mixto en conductos cerrados con esquemas en volúmenes finitos, Tecnol. y Ciencias Del Agua. 08 (2017) 127–142. doi: 10.24850/j-tyca-2017-03-08.

[12] E. Bladé, M. Sanz-Ramos, J. Dolz, J.M. Expósito-Pérez, M. Sánchez-Juny, Modelling flood propagation in the service galleries of a nuclear power plant, Nucl. Eng. Des. 352 (2019) 110180. doi: 10.1016/j.nucengdes.2019.110180.

[13] L. Cea, A. López‐Núñez, Extension of the two‐component pressure approach for modeling mixed free‐surface‐pressurized flows with the two‐dimensional shallow water equations, Int. J. Numer. Methods Fluids. 93 (2021) 628–652. doi: 10.1002/fld.4902.

[14] BOE-A-2008-755, Real Decreto 9/2008, de 11 de enero, por el que se modifica el Reglamento del Dominio Público Hidráulico, aprobado por el Real Decreto 849/1986, de 11 de abril, Boletín Of. Del Estado Núm. 14, 16 Enero 2008, Páginas 3141 a 3149. Minist. La Pres. (2008) 9. https://www.boe.es/buscar/doc.php?id=BOE-A-2008-755.

[15] ACA, Recomanacions tècniques per als estudis d’inundabilitat d’àmbit local. Guia Tècnica, Agència Catalana de l’Aigua. Generalitat de Catalunya. Març 2003, 2003. http://www.gencat.net/aca.

[16] E. Bladé, M. Gómez-Valentín, J. Dolz, J.L.L. Aragón-Hernández, G. Corestein, M. Sánchez-Juny, Integration of 1D and 2D finite volume schemes for computations of water flow in natural channels, Adv. Water Resour. 42 (2012) 17–29. doi: 10.1016/j.advwatres.2012.03.021.

[17] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer: Berlin/Heidelberg, Germany, 2009. doi: 10.1007/b79761.

[18] M.H. Chaudhry, Open channel flow, Prentice Hall, New Jersey, USA, 1993.

[19] P.L. Roe, A basis for the upwind differencing of the two-dimensional unsteady Euler equations, Numer. Methods Fluid Dyn. 2. (1986) 55–80.

[20] E. Bladé, M. Gómez-Valentín, M. Sánchez-Juny, J. Dolz, Preserving steady-state in one-dimensional finite-volume computations of river flow, J. Hydraul. Eng. 134 (2008) 1343–1347. doi: 10.1061/(ASCE)0733-9429(2008)134:9(1343).

[21] A. Preissmann, Propagation des intumescences dans les canaux et rivières (pp., in: 1st Congrèss Assoc. Fr. Calc. Grenoble, AFC, Paris, Fr. Sept., 1961: pp. 433–442.

[22] A.G. Barnett, Application of Godunov-type schemes to transient mixed flows, J. Hydraul. Res. 48 (2010) 686–688. doi: 10.1080/00221686.2010.512793.

[23] F. Kerger, P. Archambeau, S. Erpicum, B.J. Dewals, M. Pirotton, A fast universal solver for 1D continuous and discontinuous steady flows in rivers and pipes, Int. J. Numer. Methods Fluids. 66 (2011) 38–48. doi: 10.1002/fld.

[24] J. Dolz, M. Gómez, Problemática del drenaje de aguas pluviales en zonas urbanas y del estudio hidráulico de las redes de colectores, Ing. Del Agua. 1 (1994) 55–66. doi: 10.4995/ia.1994.2631.

[25] J. Despotovic, J. Plavsic, N. Stefanovic, D. Pavlovic, Inefficiency of storm water inlets as a source of urban floods, Water Sci. Technol. 51 (2005) 139–145.

[26] M. Gómez, B. Russo, Methodology to estimate hydraulic efficiency of drain inlets, Water Manag. 164 (2011) 81–90. 10.1680/wama.900070.

[27] J.Á. Aranda, C. Beneyto, M. Sánchez-Juny, E. Bladé, Efficient Design of Road Drainage Systems, Water. 13 (2021) 1661. doi: 10.3390/w13121661.

[28] M. Gómez, B. Russo, Hydraulic Efficiency of Continuous Transverse Grates for Paved Areas, J. Irrig. Drain. Eng. 135 (2009) 225–230. doi: 10.1061/(ASCE)0733-9437(2009)135:2(225).

[29] L. Cea, E. Bladé, A simple and efficient unstructured finite volume scheme for solving the shallow water equations in overland flow applications, Water Resour. Res. 51 (2015) 5464–5486. doi: 10.1002/2014WR016547.

[30] V.T. Chow, D.R. Maidment, L.W. Mays, Applied Hydrology, MCGRAW-HIL, New York, USA, 1988.

[31] M.K. Kenneth, PART 630 HYDROLOGY, in: Natl. Eng. Handb., USDA, Natural resources Conservation Service, 2010: pp. 449–456.

[32] USDA-NRCS, Part 630 Hydrology - Chapter 10, in: Natl. Eng. Handb., US Department of Agriculture, Soil Conservation Service: Washington, DC, USA, 2004: p. 79.

[33] NCRS, Hydrology, in: Natl. Eng. Handb., US Department of Agriculture, Washington DC, USA, 1972.

[34] N. Kannan, C. Santhi, J.R. Williams, J.G. Arnold, Development of a continuous soil moisture accounting procedure for curve number methodology and its behaviour with different evapotranspiration methods, Hydrol. Process. 22 (2008) 2114–2121. doi: 10.1002/hyp.6811.

[35] J. Doorenbos, W.O. Pruitt, Crop water requirements, FAO Irrig Drain. Paper 24. FAO of the United Nations, Rome, 1977.

[36] J.S. Famiglietti, E.F. Wood, Multiscale modeling of spatially variable water and energy balance processes, Water Resour. Res. 30 (1994) 3061–3078. doi: 10.1029/94WR01498.

[37] L.C. van Rijn, Sediment Transport, Part I: Bed Load Transport, J. Hydraul. Eng. 110 (1984) 1431–1456. doi: 10.1061/(ASCE)0733-9429(1984)110:10(1431).

[38] L.C. van Rijn, Mathematical modelling of morphological processes in the case of suspended sediment transport, (1987) 260.

[39] J.D. Smith, Modeling of sediment transport on continental shelves, United States, 1975. https://www.osti.gov/biblio/7156268.

[40] M. Garcia, G. Parker, Entrainment of Bed Sediment into Suspension, J. Hydraul. Eng. 117 (1991) 414–435. doi: 10.1061/(ASCE)0733-9429(1991)117:4(414).

[41] J.A. Zyserman, J. Fredsøe, Data Analysis of Bed Concentration of Suspended Sediment, J. Hydraul. Eng. 120 (1994) 1021–1042. doi: 10.1061/(ASCE)0733-9429(1994)120:9(1021).

[42] G. Parker, G. Seminara, L. Solari, Bed load at low Shields stress on arbitrarily sloping beds: Alternative entrainment formulation, Water Resour. Res. 39 (2003). doi: 10.1029/2001WR001253.

[43] E. Meyer-Peter, MüllerR., Formulas for Bed-Load transport, in: IAHSR 2nd Meet. Stock. Append. 2, IAHR, 1948.

[44] M. Wong, G. Parker, Reanalysis and Correction of Bed-Load Relation of Meyer-Peter and M?ller Using Their Own Database, J. Hydraul. Eng. 132 (2006) 1159–1168. doi: 10.1061/(ASCE)0733-9429(2006)132:11(1159).

[45] M.H. García, Sediment Transport and Morphodynamics, in: Sediment. Eng., American Society of Civil Engineers, Reston, VA, 2008: pp. 21–163. doi: 10.1061/9780784408148.ch02.

[46] G.J.C.M. Hoffmans, H.J. Verheij, Scour Manual, Routledge, 2017. doi: 10.1201/9780203740132.

[47] L. Cea, C. Legout, T. Grangeon, G. Nord, Impact of model simplifications on soil erosion predictions: application of the GLUE methodology to a distributed event-based model at the hillslope scale, Hydrol. Process. 30 (2016) 1096–1113. doi: 10.1002/hyp.10697.

[48] J. Naves, J. Anta, J. Suárez, J. Puertas, Hydraulic, wash-off and sediment transport experiments in a full-scale urban drainage physical model, Sci. Data. 7 (2020) 44. doi: 10.1038/s41597-020-0384-z.

[49] P.B. Hairsine, C.W. Rose, Modeling water erosion due to overland flow using physical principles: 2. Rill flow, Water Resour. Res. 28 (1992) 245–250. doi: 10.1029/91WR02381.

[50] P.B. Hairsine, C.W. Rose, Modeling water erosion due to overland flow using physical principles: 1. Sheet flow, Water Resour. Res. 28 (1992) 237–243. doi: 10.1029/91WR02380.

[51] J. Naves, J. Rieckermann, L. Cea, J. Puertas, J. Anta, Global and local sensitivity analysis to improve the understanding of physically-based urban wash-off models from high-resolution laboratory experiments, Sci. Total Environ. 709 (2020) 136152. doi: 10.1016/j.scitotenv.2019.136152.

[52] M. Hirano, River-bed degradation with armoring, Proc. Japan Soc. Civ. Eng. 1971 (1971) 55–65. doi: 10.2208/jscej1969.1971.195_55.

[53] M. Church, J.K. Haschenburger, What is the “active layer”?, Water Resour. Res. 53 (2017) 5–10. doi: 10.1002/2016WR019675.

[54] M. Spasojevic, F.M. Holly, Two- and Three-Dimensional Numerical Simulation of Mobile-Bed Hydrodynamics and Sedimentation, in: Sediment. Eng., American Society of Civil Engineers, Reston, VA, 2008: pp. 683–761. doi: 10.1061/9780784408148.ch15.

[55] K. Ashida, M. Michiue, Hydraulic Resistance of Flow in an Alluvia Bed and Bed Load Transport Rate, in: Proc. Japan Soc. Civ. Eng., Japan Society of Civil Engineers. No. 206 (in Japanese), 1972: pp. 59–69. https://scirp.org/reference/referencespapers.aspx?referenceid=1312271.

[56] I. V. Egiazaroff, Calculation of Nonuniform Sediment Concentrations, J. Hydraul. Div. 91 (1965) 225–247. doi: 10.1061/JYCEAJ.0001277.

[57] G. Parker, P.C. Klingeman, D.G. McLean, Bedload and size distribution in paved gravel-bed streams., J. Hydraul. Div. - ASCE. (1982). doi: 10.1061/(ASCE)0733-9429(1983)109:5(793).

[58] P.R. Wilcock, J.C. Crowe, Surface-based Transport Model for Mixed-Size Sediment, J. Hydraul. Eng. 129 (2003) 120–128. doi: 10.1061/(ASCE)0733-9429(2003)129:2(120).