Análisis hidro-sedimentológico 2d del comportamiento de un campo de espigones para la protección de márgenes en un río meandriforme

ABSTRACT

Alluvial rivers show a highly dynamic morphological patterns, which cause problems in the operation and maintenance of infrastructure located in their environment. In these cases, it is necessary to design structures to control bank morphodynamic processes. This article describes the hydro-sedimentology of the meandering channel called río del Valle (Argentina), for which bank protection was implemented on its outer bank. This study zone was selected for the cutting of an embankment located on the outer bank of the channel due to the excessive lateral migration of the meander. Flow structure and bed morphology patterns are analyzed and discussed using the results of the application of the TELEMAC 2D model coupled to the Sisyphe morphological module. These results allow us to identify vector patterns of velocity and streamlines in the groyne field, showing the development of a recirculation flow between groynes, as well as the influence in the displacement of the core of maximum velocity towards the center of the channel along the bend apex. The bed evolution shows the vertical regrowth of the bed in the groyne field and the reduction of the point bar in the bend entrance.

RESUMEN

En ocasiones, los ríos de llanura exhiben características morfológicas altamente dinámicas, que producen inconvenientes en la operación y mantenimiento de obras de infraestructura ubicadas en su entorno. En esos casos, es necesario disponer de obras de control de los procesos erosivos. En este artículo se describe el comportamiento hidro-sedimentológico de una curva en el rio del Valle (Argentina), sobre la cual se implantó una obra de protección sobre su margen externa. Dicha zona se seleccionó por el corte de un terraplén localizado en la margen externa del cauce debido a la excesiva migración lateral del meandro. Se analizan los patrones de la estructura del flujo y morfología de fondo mediante los resultados de la implementación del modelo TELEMAC 2D acoplado al módulo morfológico Sisyphe. Estos resultados permitieron comparar e identificar los patrones vectoriales de velocidades y las líneas de corriente en el campo de espigones, mostrando el desarrollo de un flujo de recirculación entre espigones, como así también el desplazamiento del núcleo de máximas velocidades hacia la región central de la curva. La evolución del fondo muestra el recrecimiento vertical del fondo en la zona de espigones y la reducción de barra empuntada a lo largo de la curva.

PALABRAS CLAVE:

• protección de márgenes
• simulación numérica
• estructura de flujo
• transporte de sedimento
• morfología de fondo

KEYWORDS:

• bank protection
• numerical modeling
• flow structure
• sediment transport
• bed morphology

1. Introducción

Los ríos de llanura exhiben características morfológicas altamente cambiantes en el tiempo, con una marcada actividad morfodinámica evidenciada por la presencia de meandros, con recurrentes eventos de elongaciones, migraciones longitudinales y laterales y cortes naturales por estrangulamiento y desbordes (o procesos combinados) (Hooke, 2020). Estos cambios morfológicos ocasionan inconvenientes en la operación y mantenimiento de obras de infraestructura ubicadas en su vecindad inmediata, lo cual hace necesario la disposición de estructuras complementarias que permitan la mitigación y control de los procesos erosivos presentes en estos cauces.

El gran desarrollo en los últimos 40 años de los modelos hidro-sedimentológicos bidimensional junto a la capacidad de cómputo de los recientes microprocesadores, ha potenciado la representación precisa de los fenómenos fluviales y por ende ha permitido tomar decisiones a la hora de poner en marcha una obra en el ámbito ingenieril (Villaret et al., 2013; Fernández-Pato et al., 2019). Aun así la aplicabilidad de estos modelos se encuentra limitada a la información con la que se cuenta y a la complejidad del fenómeno que se quiera representar, incrementando lo requerido a medida que aumentamos las dimensiones de resolución. En este sentido, los cauces meandriformes presentan una compleja hidrodinámica y sedimentología a lo largo de sus curvas (Dominguez Ruben et al., 2021). Si bien existen una serie de estudios donde se han representado la influencia de los espigones, algunos de ellos lo han hecho sobre tramos rectos (Sukhodolov, 2014) y en la generalidad de los casos en modelos de laboratorio (Uijttewaal, 2005). Poco se conoce del comportamiento hidro-sedimentológico presente sobre obras de protección en tramos curvos, donde mayormente y dado su dinámica migratoria son necesarias.

El objetivo del artículo es presentar y discutir los resultados de la aplicación del modelo numérico bidimensional TELEMAC 2D, representando la influencia de un campo de espigones, localizado sobre la margen externa de una curva en el río del Valle en la región Norte de Argentina. La particularidad de la zona de estudio se encuentra dada por la complejidad de llevar a cabo relevamientos ya que sus eventos se producen en un periodo relativamente corto (20 días) y en los mismos se movilizan grandes volúmenes de material sólido (en suspensión y fondo). Los resultados abordados permiten optimizar el desempeño del conjunto de obras de estabilización de márgenes fluviales, gracias a la ocurrencia de eventos hidrológicos similares a las condiciones de diseño asumidas inicialmente.

2. Zona de estudio

El río del Valle tiene su naciente en el sistema orográfico de las sierras Cresta de Gallo y del Piquete (Parque Nacional El Rey), escurre hacia el este y atraviesa la traza del FFCCGB (Ferrocarril General Belgrano) a la altura de la localidad de Las Lajitas, Salta, Argentina. Finalmente confluye, a través de un frente de escorrentía difusa, con el río Dorado aproximadamente 30 km aguas abajo, al sureste de la localidad de Apolinario Saravia, sobre la Ruta Provincial N° 5 (R.P. No. 5, Provincia de Salta). Posee una cuenca de 1126 km², mientras que la longitud del curso principal es de 73 km.

El río del Valle en cercanías a la ciudad de Las Lajitas (Figura 1), en su tramo próximo a la traza de la R.P. No. 5 y del Ferrocarril FFCCGB, exhibe un comportamiento fluvio-morfológico altamente dinámico, caracterizado por un patrón de alineamiento en planta meandriforme con fuertes tendencias a desarrollar procesos de migración lateral y longitudinal, acompañados por frecuentes cortes de meandros. El cauce presenta un caudal base de 20 m³/s con una variación diaria que alcanza una relación Q_max/Q_min = 150 (Estación El Ceibal. Fuente: Sistema Nacional de Información Hídrica). El río del Valle se caracteriza por presentar eventos hidrológicos típicos de un subsistema de la cuenca del río Bermejo, donde se movilizan grandes volúmenes de caudales líquidos en no más de 20 días (Farias, 2017). En el estudio de Farias (2017) se indica que para una crecida de recurrencia de ocho años se alcanzan picos de 500 m³/s y generan la movilización de un total de 0.1334 hm³ de material sólido.

Figura 1. Ubicación de zona de estudio (Provincia de Salta, Argentina)

El ancho del rio del Valle activo es moderadamente variable (en el rango de 50 m a 80 m, con algunos sectores de sobreanchos, particularmente sobre el ápice de la curva) de meandros de gran amplitud. Dada la presencia de curvas la sección transversal es de tipo asimétrica con una profundidad media variable entre 3 m y 7 m. La tasa de migración lateral media del cauce es de 60 m/año, aunque en ocasiones los corrimientos alcanzan valores superiores dada la intensidad y duración de las crecidas (Farias, 2017). En el tramo de aguas abajo del Ferrocarril FFCCGB, en un lapso de algo más de 10 años (entre 2005 y 2015) un meandro con forma en planta compuesta de considerable amplitud afectó con procesos erosivos al terraplén ferroviario en el sector de máxima concavidad (Figura 2), produciendo el corte del servicio. Nótese que las márgenes del cauce sobre la planicie aluvial se encuentran confinadas por las actividades de movimiento de suelos vinculadas a los manejos de campos productivos en la región, razón por la cual no se puede afirmar que el desarrollo longitudinal y lateral de los meandros sea totalmente libre.

Figura 2. Evolución morfológica del río del Valle entre los años 2002 y 2016. Digitalización de márgenes

Aproximadamente 800 m aguas abajo del cruce con el FFCCGB, la faja de divagación se re-orienta y tras un giro horario de casi 90º atraviesa a la traza de la R.P. No. 5 con un sentido dominante NW-SE, desarrollando en ese sector un proceso de ‘apilamiento’ (stacking) de meandros inmediatamente aguas abajo de la traza del FFCCGB. La curva de estudio comprende justamente a la desarrollada entre la traza del FFCCGB y la R.P. No. 5.

Las márgenes del río del Valle en el tramo estudiado están constituidas por suelos limo-arcillosos y limoso (CL-ML), según el sistema SUCS. En la zona particular en la que se han producido las erosiones, hay preeminencia de suelos del tipo CL, según los estudios geotécnicos llevados a cabo. El material de fondo se compone por una arena media de D₅₀ = 0.240 mm. Por su parte los perfiles de perforaciones realizadas en la zona indican la presencia de tres estratos, el primero de 1.8 m de profundidad de material SM (limo-arenoso) no plásticas con un ángulo de reposo de 18°, un segundo estrato hasta los 3 m de profundidad del mismo material pero con un ángulo de reposo de 20° y por último hasta los 5.5 m un material fino de tipo CL (arcilla de baja plasticidad).

La Figura 3 muestra las fotografías tomadas para un estado hidrológico bajo (año 2017) en la curva en estudio. Se desprende de la misma la acción erosiva del cauce sobre el terraplén del FFCCGB en la zona de margen izquierda del sector de máxima concavidad del meandro aguas abajo de la traza del ferrocarril.

Figura 3. Procesos de erosión lateral y afectación a la traza ferroviaria (fotografías octubre 2017). (a) Vista (hacia aguas abajo desde el terraplén del FFCCGB) del extremo cóncavo de la curva que ha producido la erosión sobre margen izquierda y el acercamiento al eje de la traza, con erosiones transversales importantes sobre el terraplén. (b) Vista hacia aguas arriba. Se observa estructuras turbulentas del flujo en la zona de la curva, que gobierna la mecánica de los procesos erosivos en el segmento de aproximación del flujo hacia el terraplén, desde donde se realiza la toma fotográfica

3. Diseño y funcionamiento esperables de la obra de protección y control de erosión lateral

Ante la situación de marcada actividad morfodinámica, se dispusieron acciones de mitigación y control de los procesos erosivos en el sector de la margen externa de la curva del río del Valle (Figura 1). Teniendo en cuenta los rasgos morfológicos del río, se definió implantar una obra de estabilización de la margen izquierda, con la finalidad de estabilizar el curso fluvial en ese segmento curvo. La obra principal está conformada por un campo de siete espigones deflectantes materializados en gaviones (numerados progresivamente desde aguas arriba hacia aguas abajo, para facilitar su identificación), en combinación con una cubierta continua en la zona afectada, tal como se indica en el esquema de la Figura 4. El diseño de los espigones se ha realizado de acuerdo a las formulaciones empíricas y lineamientos generales propuestos por Maza Alvarez y García Flores (1994). Cabe resaltar que el diseño final además se basó en experiencias en obras implantadas sobre ambientes fluviales similares, cuyo rendimiento fue positivo.

Figura 4. Planimetría y cortes de la conformación de espigones de obras de control de erosión en el río del Valle

La disposición en planta de los espigones en el río del Valle se definió adoptando un ángulo de expansión variable entre 9º y 14º. El ángulo de orientación de los espigones respecto a la tangente a la línea de defensa se adoptó en 70º. La longitud de empotramiento de los espigones en la margen es de aproximadamente ¼ de la longitud de trabajo. La altura de los espigones es creciente desde la cabeza donde tienen una altura de 0.50 m hacia la raíz donde alcanzan una altura de 3.50 m.

Los espigones se apoyan sobre una colchoneta que en el sector de la cabeza se ensancha y prolonga hacia el cauce en una longitud de 12 m para proteger la erosión.

La Figura 5 presenta un esquema ilustrativo del efecto esperado en cuanto al comportamiento de las obras implantadas, tanto para las etapas de ejecución como de estabilización final. En la situación natural imperturbada (esquema ‘1’), se cuenta con una sección de geometría asimétrica (línea gruesa, marrón oscura), con mayores profundidades recostadas sobre la región de la batiente cóncava de la curva (margen izquierda en este caso), tal como se visualizan las secciones relevadas en campo y modeladas en el estudio. La situación ‘2’ representa la condición esperable, inmediatamente a posteriori de la puesta en operación de la obra de protección. En rigor, estos elementos representan inicialmente una obstrucción al libre escurrimiento, lo cual se traduce en un efecto de remanso que produce una sobreelevación del nivel de la superficie libre del flujo (línea de trazos celeste en el esquema). A partir de este momento, comienzan a desarrollarse los procesos morfológicos de erosión general, erosión local en el extremo de los espigones (aspecto a controlar con los delantales de colchonetas tipo reno proyectadas) y erosiones transversales, en respuesta a la obstrucción del área hidráulica en el sub-sector de margen izquierda.

Figura 5. ’- Esquema con la modalidad de funcionamiento previsto para la obra en tres etapas

El esquema ‘3’ de la Figura 5 ilustra la situación final de equilibrio esperable, luego del desarrollo completo del proceso de erosión lateral descripto en la situación ‘2’. La línea continúa anaranjada describe la nueva geometría hidráulica esperable, respecto de la situación inicial (línea de trazos, color marrón). Es decir, se produce una ‘ganancia de área hidráulica’, que alivia en parte las condiciones de flujo, generando una disminución de los niveles de la superficie libre (línea magenta) tendiente a la convergencia con los niveles correspondientes a la situación inicial (sin perturbación).

Una vez alcanzada esta condición final, la cual se potenciaría con la construcción de un canal piloto para desviar los flujos bajos del río en la instancia de construcción de los espigones, las recurrencias de los caudales de desborde, tenderán a los mismos valores que las actuales (condición ‘natural’ del río en el sector).

3.1 Descripción del modelo matemático y su implementación

3.1.1 TELEMAC 2D

El modelo numérico utilizado para la evaluación del comportamiento hidráulico en el tramo del río bajo análisis fue TELEMAC-MASCARET (Hervouet, 2007). La justificación del uso de este modelo se encuentra asociada, entre otras cosas, al esquema de elementos finitos que implementa, permitiendo controlar las condiciones de borde y logrando una gran adaptabilidad geométrica y versatilidad que posibilitó la modelación numérica bidimensional sobre un ambiente natural y complejo, como es el río del Valle en el tramo seleccionado.

El módulo TELEMAC 2D se utiliza para el cálculo de flujo en superficie libre bidimensional en el plano horizontal. Las principales variables que se obtienen son profundidad de agua y velocidad sobre las dos componentes del plano (x,y). TELEMAC 2D resuelve las ecuaciones de Saint Venant usando el método de elementos finitos sobre una malla triangular, mediante la resolución de las siguientes cuatro ecuaciones:

(1) $\frac{\mathrm{\partial }H}{\mathrm{\partial }t}+\vec{u}\cdot \vec{\mathrm{\nabla }}\left(H\right)+Hdiv\left(\vec{u}\right)=S_h$(1)

(2) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\overrightarrow{u}\bullet \overrightarrow{\mathrm{\nabla }}\left(u\right)=-g\frac{\mathrm{\partial }{Z}_s}{\mathrm{\partial }x}+S_x+\frac{1}{H}div\left(H{\nu }_f\overrightarrow{\mathrm{\nabla }}u\right)$(2)

(3) $\frac{\mathrm{\partial }v}{\mathrm{\partial }t}+\vec{u}\cdot \vec{\mathrm{\nabla }}\left(v\right)=-g\frac{\mathrm{\partial }{Z}_s}{\mathrm{\partial }y}+S_y+\frac{1}{H}div\left(H{v}_f\vec{\mathrm{\nabla }}v\right)$(3)

(4) $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\overrightarrow{u}\bullet \overrightarrow{\mathrm{\nabla }}\left(T\right)={S_T}_r+\frac{1}{H}div\left(H{\nu }_T\overrightarrow{\mathrm{\nabla }}T\right)$(4)

donde H(profundidad de agua), u, v (velocidad promediada en la vertical en las direcciones x e y) y T (trazador pasivo) son nuestras variables desconocidas.

Los términos de las ecuaciones previas son calculados en uno o más pasos (en el caso del término de advección por el método de las características, entre otros) de la siguiente manera: (i) Advección de H, u y v; (ii) Propagación, difusión y términos fuentes de la ecuación de conservación (Ec. 1); (iii) Difusión y términos fuentes de la ecuación de transporte de trazadores.

En cuanto a la viscosidad turbulenta puede ser ingresada por el usuario (viscosidad constante) o determinada a través del modelo de cierre de turbulencia k-ε (entre otros), en donde se evalúa la generación de energía turbulenta con la disipación (Hervouet, 2007).

Diferentes expresiones teóricas para la representación de rugosidad de fondo son posibles implementar en TELEMAC 2D, entre los que se encuentran: Haaland, Chézy, Strickler, Manning, Nikuradse, Ley de la pared (solo en condiciones de borde) y Colebrooke-White (Hervouet, 2007).

3.1.2 Sisyphe

TELEMAC 2D permite acoplarse a un módulo de cálculo de transporte de sedimento denominado Sisyphe (Villaret et al., 2013; Tassi, 2017). De esta manera es posible realizar cálculos del transporte de sedimento para los mecanismos existentes, hablamos de fondo, en suspensión o carga de lavado. Uno de los parámetros más importantes que vincula ambos módulos (TELEMAC 2D y Sisyphe) es la tensión de corte de fondo $({\tau }_b)$, definida:

(5) ${\tau }_b=0.5\rho C_f{\left(u^2+v^2\right)}^{0.5}$(5)

donde u y v son extraídas de TELEMAC 2D, ρ es la densidad del agua y $C_f$ es el coeficiente de fricción relacionado a la rugosidad del fondo (ver Tassi, 2017). La velocidad de corte queda definida de la siguiente manera $u_{\ast }=\sqrt{{\tau }_b/\rho }$.

Para el cálculo del caudal transportado por fondo, Sisyphe ofrece diferentes expresiones, muchas de ellas asumen la existencia de una capacidad transporte iniciada con el desplazamiento de la partícula. Estas fórmulas relacionan la tensión adimensional de corte de fondo debida al flujo y la tensión crítica de corte o parámetro de Shields ($\theta )$ para el cálculo de una tasa de transporte adimensional (${\mathrm{\Phi }}_b$), como:

(6) ${\mathrm{\Phi }}_b=\frac{Q_b}{\sqrt{g\left(s-1\right)D_{50}^3}}$(6)

donde $Q_b$ es la tasa de transporte de fondo por unidad de ancho, ρ_s es la densidad del sedimento, s=ρ_s/ρ la densidad relativa, D₅₀ diámetro medio del sedimento modelado y g aceleración de la gravedad. El parametro ${\mathrm{\Phi }}_b$ se encuentra correlacionado al parámetro adimensional conocido como número Shields $\left(\theta \right)$, definido como:

(7) $\theta =\frac{\mu {\tau }_b}{\left({\rho }_s-\rho \right)g{D}_{50}}$(7)

siendo μ un factor de corrección de la rugosidad del fondo (dependiente de las formas del fondo presentes, ver Tassi, 2017).

Existen diferentes formulaciones semiempíricas para el cálculo de ${\mathrm{\Phi }}_b$, entre las que se encuentran: basadas en conceptos de energía (Engelund y Hansen, 1967) o de derivaciones estadísticas (Einstein, 1950) o basadas en la capacidad de transporte como van Rijn (1984):

(8) ${\mathrm{\Phi }}_b=0.053D_{\ast }^{-0.3}{\left(\frac{\theta -{\theta }_{cr}}{{\theta }_{cr}}\right)}^{2.1}$(8)

donde $D_{\ast }$ es el diámetro adimensionalizado por la viscosidad del agua y ${\theta }_{cr}$ es el número crítico de Shields. La ecuación bidimensional de transporte para una concentración media en la vertical se obtiene de la integración de la ecuación de transporte 3D. Si aplicamos la regla de Leibniz adoptando adecuadas condiciones de borde y asumiendo un fondo suficientemente delgado, la ecuación de transporte de sedimento suspendido resulta:

(9) $\frac{\mathrm{\partial }C}{\mathrm{\partial }t}+u\frac{\mathrm{\partial }C}{\mathrm{\partial }x}+v\frac{\mathrm{\partial }C}{\mathrm{\partial }y}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left({\nu }_T\frac{\mathrm{\partial }C}{\mathrm{\partial }x}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left({\nu }_T\frac{\mathrm{\partial }C}{\mathrm{\partial }y}\right)+\frac{{\left(E_T-D_T\right)}_{z=a}}{H}$(9)

donde $C=C\left(x,y,t\right)$ es la concentración media en la vertical expresada en % del volumen, $a$ es la altura de referencia $a=0.05H$ y $v_T$ es el coeficiente de difusión turbulenta del trazador (en nuestro caso del sedimento). Sisyphe utiliza para el cálculo del termino advectivo el método de las características por defecto, aunque permite optar por otros como SUPG (por siglas en inglés, Streamline Upwind Petrov Galerkin), Esquema Conservativo, entre otros (para mas detalles ver Hervouet, 2007).

En cuanto al término difusivo, el valor de $v_T$ puede ser calculado de manera constante longitudinal y transversalmente (ingresado por el usuario), utilizar el mismo que para flujo (trazador pasivo), $v_T=v_f$ extraído de TELEMAC 2D o puede ser calculado en base al modelo de Elder (Tassi, 2017).

Para el cálculo del último término de la Ec. 9 se acepta la formulación de Celik y Rodi (1988) la cual se basa en el concepto de concentración de equilibrio:

(10) ${\left(E_T-D_T\right)}_{ref}=w_s\left(C_{eq}-C_{ref}\right)\hspace{0.17em}$(10)

E_T y D_T son tasas de erosión y deposición respectivamente, calculadas sobre el nivel de referencia. La tasa de deposición $D_T=w_s{C}_{ref}$, donde $w_s$ es la velocidad de caída de la partículas y $C_{ref}$ es la concentración de referencia en la interfaz fondo y suspensión ($z=z_{ref}$). La tasa de erosión de material no cohesivo es $E_T=w_s{C}_{eq}$, donde $C_{eq}$ es la concentración de equilibrio posible de calcular mediante diferentes ecuaciones empíricas, entre las que se encuentran la siguiente expresión de van Rijn (1984):

(11) $C_{eq}=\frac{0.015D_{50}{\left(\frac{\theta { }^{\mathrm{\prime }}}{{\theta }_{cr}}-1\right)}^{\frac{3}{2}}}{z_{ref}{D}_{\ast }^{0.3}}$(11)

donde $\theta { }^{\mathrm{\prime }}=\mu \theta$ es la tensión generada por la fricción del fondo.

Sisyphe asume que el perfil vertical de concentraciones desarrollado por la partícula es el propuesto por Rouse (1937):

(12) $C\left(z\right)=C_{ref}{\left(\frac{z-H}{z}\frac{a}{a-H}\right)}^{N_R}$(12)

donde $C_{ref}$ es la concentración de referencia y $N_R$ es definido como el número de Rouse de la siguiente manera:

(13) $N_R=\frac{w_s}{\kappa u_{\ast }}$(13)

κ es la constante de von Karman (= 0.4). Mediante la integración de la ecuación de Rouse (Ec. 12) es posible determinar la concentración de referencia $C_{ref}=FC$, donde:

(14) $F^{-1}={\left(\frac{z_{ref}}{H}\right)}^{N_R}\underset{z_{ref}/H}{\overset{1}{\int }}{\left(\frac{1-u}{u}\right)}^{N_R}du$(14)

En Sisyphe, la expresión previa es computada de la siguiente manera:

(15) $\begin{array}{rl}& F^{-1}=\left\{\frac{1}{1-Z}{B}^{N_R\left(1-B\left(1-N_R\right)\right)}\right.siR\ne 1\\ & -BlogBsiR=1\end{array}$(15)

con $B=z_{ref}/H$.

Por último, es relevante destacar que la evolución del fondo es calculada mediante la ecuación de Exner donde se introduce las Ecs. 8 y 10 como así también la porosidad del material (Tassi, 2017).

3.1.3 Modelo Digital del Terreno (MDT)

Para la generación del MDT se utilizaron los datos batimétricos de campo bajo el sistema de referencia Gauss Krugger (Faja 4), medidos durante diciembre del 2017. El MDT se confeccionó en el paquete de software ArcGIS (ESRI, 2011) el cual, entre otros beneficios operativos, permite al usuario modificar el sistema de proyecciones y trabajar en el más adecuado. Para la generación de una primera superficie se recurrió a la función TIN (Triangular Irregular Network), la cual genera una red de triángulos irregulares conectados por los nodos de datos individuales. De esta manera un TIN conserva toda la precisión de los datos de entrada al mismo tiempo que infiere los valores entre los puntos conocidos. Existen diversos métodos de triangulación para conformar esta malla. El método aquí utilizado, y uno de los más populares, fue el de Delaunay.

Una vez generado el TIN, se prosiguió con la modificación o extracción de datos inválidos. Esto permite corregir errores producidos en la toma de datos o puntos incongruentes en la malla. Logrado el TIN filtrado, se lo convirtió a Raster y de éste a Puntos. A continuación se suavizó dichos Puntos con la función de ‘interpolación kernel’ propia del paquete de software ArcGIS (ESRI, 2011). Esta función permitió suavizar la nube de puntos con el método de regresión de crestas (Hoerl y Kennard, 1970), el que evita cualquier inestabilidad de picos en los datos ingresados. El producto de esta función es lo que finalmente se conoce como MDT, el cual es una conjunto de puntos (ordenados en una cuadrícula), y provee una sensación de continuidad visual cuando se realiza una proyección 3D (aunque en rigor, es un modelo discreto de la topografía del terreno, sumergido en éste caso, puesto que representa el lecho del río). Ver Figura 6a

Figura 6. (a) Modelo digital del terreno de toda la región de estudio con espigones implantados; (b) detalle de la zona de espigones; (c) malla de todo el dominio para los casos 0 y 1 y (d) detalle de la malla en la zona de espigones

.

3.1.4 Dominio computacional y malla

El dominio computacional se definió en base a los objetivos de estudio. El dominio se generó con el ente polilínea del paquete de software ArcGIS (ESRI, 2011), con la ayuda de la imagen satelital referente al periodo temporal que se esta analizando (diciembre del año 2017). La malla para la región de estudio (Figura 6b) fue generada con el paquete de software BlueKenue (Bluekenue, 2013). La metodología de construcción de la malla que usa dicho software es la condición de Delaunay. Los archivos de entrada para la generación de mallas son los contornos de la región que se quiere modelar. La practicidad del sistema BlueKenue se basa en que los archivos de entrada pueden ser importados por un gran número de paquetes de software de tipo GIS. Además BlueKenue permite modificar (mover) o incrementar los nodos y elementos de manera manual con las herramientas propias. Esto es sumamente útil en ciertas ocasiones, permitiendo modificar el mallado automático, o densificando nodos en regiones de interés.

Para llevar a cabo las modelaciones se generaron dos mallas; i) mapeando el fondo en su condición natural, e ii) incorporando la presencia de los espigones predimensionados.

Para la determinación de la distancia entre nodos se evaluaron diferentes geometrías. Inicialmente, el espaciamiento entre nodos (Δ_m) de la malla se planteó con base en la siguiente ecuación (Morell et al., 2014):

(16) ${\mathrm{\Delta }}_m\cong \sqrt{\frac{2\mathrm{\Omega }}{N_e}}$(16)

donde Ω es el área del dominio computacional y N_e es el número de elementos (limitado por la capacidad computacional disponible). En función de las dimensiones del dominio y las capacidades de modelación con las que se cuenta, se determinó una separación de nodos de 5 m. Adicionalmente, se generó una submalla de 1 m para la región coincidente con los espigones, de manera de lograr una mejor representación de los mismos. Para la definición de esta configuración no estructurada se realizó un análisis de independencia de malla (Morell et al., 2014).

En cuanto al paso de tiempo, se evaluó la condición de Courant (C_o < 1), siendo ésta necesaria para la convergencia, mientras se resuelven numéricamente las ecuaciones de avance temporal de la solución por el método de las elementos finitos (Hervouet, 2007). De esta manera se determinó un paso de tiempo de Δt = 0.05s.

3.1.5 Condiciones hidrodinámicas y sedimentológicas

El dominio de modelación comprende un tramo de 947 m de longitud, donde se incluyen dos curvas, delimitado aguas arriba por el puente del FFCCGB y aguas abajo por el punto de inflexión entre el final de la curva en cuyo extremo se desarrollan los procesos de erosión lateral que afectaron el terraplén, y el inicio del meandro siguiente, antes del puente vial sobre R.P. No. 5. Se modelaron dos geometrías, la primera que representa la condición natural del cauce sin la presencia de espigones y la segunda con la obra implantada en la margen externa de la curva (Caso0 y Caso1, respectivamente). En ambos casos se simuló la condición hidrológica de una crecida de recurrencia de 10 años igual a 571 m³/s.

Dada la escasa información con la que se contaba de los eventos hidro-sedimentológicos, las condiciones de borde aguas abajo (alturas) se establecieron en base a un modelo unidimensional HEC-RAS, que permitieron validar las profundidades del modelo con las observadas en campo (Farias, 2017). Se usó el modelo TELEMAC 2D acoplado al módulo Sisyphe para la representación de la evolución del fondo en ambos casos de estudio. El tiempo de modelación se definió de modo tal de conseguir la estabilización del transporte de flujo y sedimento en todo el dominio modelado. Se estableció un coeficiente de resistencia del flujo, aplicando un coeficiente de rugosidad de Manning n = 0.030 para el fondo del cauce y sobre la región de los espigones se impuso un valor superior (n = 0.075).

Para la representación de las corrientes secundarias se plantea la corrección del flujo propuesta por Finnie et al. (1999). El modelo de turbulencia utilizado para todos los casos fue k-ε (Hervouet, 2007).

La estimación del transporte de sedimentos y evolución del lecho fueron tratados separando la fracción de carga de fondo y suspensión (material no cohesivo). Ambas fracciones fueron representadas por una sola clase de sedimento (D₅₀ = 0.24 mm). El modelo morfológico de evolución de fondo se aplicó en todo el dominio en el Caso0, mientras que para el Caso1 se fijó el fondo solo sobre la zona de los espigones. Para el cálculo del transporte de fondo se usó la expresión de van Rijn (1984) a través de la Ec. 8. Por su parte, para el transporte en suspensión la estimación de $C_{eq}$ se llevó a cabo con la expresión de van Rijn (1984) Ec. 11 (Tassi, 2017). Nótese que la implementación de las ecuaciones de van Rijn (1984) se encuentra basada en la validez de las mismas para la clase de material sedimentológico modelado. Adicionalmente, se impuso como condición de entrada al dominio, el ingreso de material en suspensión y fondo. La carga de lavado no fue estimada en ambos casos. Por último, se definió una capa activa de material de fondo igual a 3 m, en base a perforaciones de campo con las que se contaba (Farias, 2017).

4. Resultados

A continuación se presentan los resultados de ambos casos de estudio (Tabla 1). Nótese que la elección del caudal a modelar se debe a la ocurrencia de una crecida de tal magnitud entre los meses de enero-febrero del año 2019.

Tabla 1. Resumen de condiciones de bordes, geometría y tiempo de recurrencia modeladas para los 2 casos de estudio

Casos	Condición	TR [años]	Caudal de ingreso al dominio [m^3/s]	Número de elementos
Caso0	Sin espigones	10	571	9030
Caso1	Con espigones	10	571	9030

La Figura 7a,b muestra las velocidades medias en la vertical resultantes de la simulación para ambos casos de estudios (sin y con espigones respectivamente) en un contorno de colores conveniente. En el Caso0 se distingue como el núcleo de máxima velocidad se recuesta sobre la región externa (sobre la primera curva modelada, C1), donde el cauce no migra libremente debido a la densa vegetación de la margen externa y luego de esto ‘salta’ hacia la margen opuesta para mantenerse sobre la misma hasta la finalización de la segunda curva modelada (C2). Este comportamiento de la estructura de flujo concuerda con lo encontrado en mediciones de campo por Dominguez Ruben et al. (2021). Por otra parte, en el Caso1 los espigones desplazan al núcleo de máximas velocidades hacia el interior del cauce. Este desplazamiento reduce notablemente la acción erosiva sobre la margen externa de la curva. En este mismo caso se puede observar patrones de recirculación mayormente desarrollados entre los E4-E5 y aguas abajo del E7 (Figura 7b,c). Posiblemente esto se debe a la combinación de la dirección dominante del flujo a lo largo de la curva y la de los espigones (Figura 7c). Esto genera un gradiente de velocidad localizado en el extremo de los espigones. Cabe resaltar que para esta condición hidrológica el nivel del agua sobrepasa la altura de los espigones, por lo tanto se observan velocidades bajas sobre los mismos (Figura 7c). Nótese además como el núcleo de máximas velocidades inciden en la región externa de la curva C2, en particular en la finalización de la misma impactando (Caso1) sobre la margen vegetada.

Figura 7. Magnitud de velocidades medias en la vertical para la condición de crecida modelada (571 m³/s) en ambas geometrías modeladas (a) sin y (b) con espigones y (c) detalle del campo de velocidades en Caso1 entre E4-E5

Respecto a la región interna de las curvas, el comportamiento de ambos casos de estudio presentan similitudes, en particular sobre C2. Aunque se distingue una leve reducción de velocidad en el Caso0, que se extiende sobre toda la región interna hasta su finalización. No se observa un comportamiento claro de recirculación (en ninguna de las curvas). Esto podría deberse a la geometría de gradual expansión que presenta C2 o al estado hidrológico que evita o confina (a una pequeña región) la recirculación de flujo sobre la barra empuntada (point bar). En este sentido sobre el inicio de C2, para ambos casos modelados, se observa una zona de bajas velocidades. Esto se debe al cambio brusco de la dirección del núcleo de máxima velocidad a la entrada de la curva.

La Figura 8 presenta el comportamiento de la evolución del fondo en el último paso de tiempo simulado para ambos casos de estudio. En concordancia con el patrón de velocidades indicado anteriormente, se observa una leve profundización de la zona central en C2 para el Caso1 (Figura 7b), que limita el desarrollo lateral de la barra empuntada. Por el contrario, el Caso0 muestra una expansión de la barra empuntada, lo cual podría potenciar el desplazamiento del núcleo de máximas velocidades a la región externa (más aún en ausencia de espigones) promoviendo la migración de la curva. Por otro lado para el Caso1, dado el desarrollo de bajas velocidades en la región entre espigones, se produce la deposición de material en la zona de gradientes elevados de velocidad (Figura 7c). Esto permite el desarrollo de una barra longitudinal entre los extremos de E1 a E3, que confina el thalweg entre el extremo de los espigones y la barra empuntada. Por último, se destaca que la presencia de espigones parecería evitar el desarrollo completo de la barra empuntada sobre C1 (ver zona E1 en Figura 8b).

Figura 8. Configuración final del fondo para la condición de crecida modelada (571 m³/s) en ambas geometrías (a) sin y (b) con espigones

Se desprende de las Figuras 7 y 8 que el modelo logra representar de manera razonable los campos de líneas de corriente y la evolución esperable del lecho fluvial, con una disminución considerable de las velocidades del flujo en la zona de espigones (propiciando la futura sedimentación en el sector) y la formación incipiente de una barra empuntada de acreción vertical y horizontal hacia la margen derecha.

5. Discusión

La visualización del comportamiento de las obras se pudo llevar a cabo a través de un monitoreo de las mismas, durante su construcción y su operación inicial. Gracias a la disponibilidad de imágenes de la aplicación Google Earth® se puede resumir la historia reciente del caso presentado, a través de una secuencia de imágenes en los últimos 4 años. En la Figura 9a se muestra el tramo de estudio en febrero de 2016, en donde se aprecia claramente el proceso de erosión de margen izquierda con el corte del terraplén ferroviario que obligó al diseño de las obras de control de erosión. La imagen de noviembre de 2018 (Figura 9b) muestra la obra recién concluida, al inicio de la temporada de lluvias en la región, mientras que la Figura 9c (capturada en enero de 2020) permite apreciar el desempeño inicial de las obras luego de algo más de un año de operación. En particular, durante los meses de enero y febrero de 2019 se presentaron en el río crecidas importantes, muy cercanas a la condición asociada a los 10 años de recurrencia (~500 m³/s), y el comportamiento hidrodinámico y morfológico que exhibió el sistema resultó sorprendentemente similar al simulado con el modelo bidimensional. En efecto, puede observarse la deposición inmediatamente aguas abajo de cada espigón (en particular sobre los primero E1 a E4), también la formación de una vía preferencial de flujo recostado sobre el centro del cauce aguas abajo del Espigón E1, la reducción de la barra empuntada al inicio de la curva y la formación de la barra empuntada de acreción en margen derecha (aguas abajo del ápice de C2), frente al campo de espigones (ver Figura 9c).

Figura 9. (a) Erosión de márgenes antes de construir la protección (24/02/2016); (b) condición del río con obra recién construida (14/11/2018) y (c) condición posterior a la crecida registrada en enero-febrero del 2019 (06/01/2020)

Los patrones de comportamiento observados muestran la eficacia de la herramienta computacional usada (modelo TELEMAC 2D-Sisyphe) en la predicción de la hidrodinámica y procesos morfológicos asociados al esquema de obra propuesto. Resulta evidente que el modelo pudo capturar y representar adecuadamente aspectos que no pueden ser modelados en un paquete 1D, que sólo produce indicios a nivel cualitativo de algunos patrones hidráulicos en escala macroscópica. Aun así, en virtud del dinamismo del cauce, en particular para eventos de crecidas, el modelo actual presenta limitaciones en cuanto a la delimitación de la malla de estudio. Una posible solución sería aplicar el módulo de malla dinámica (Langendoen et al., 2016). Esto permitiría representar de una manera más precisa los cambios morfológicos, librando los límites del dominio computacional impuestos por las condiciones geométricas de una malla fija. Esto podría ser beneficioso para una modelación temporal del orden de décadas en cauces de gran dinamismo.

6. Conclusiones

Los resultados de la modelación numérica han confirmado las principales hipótesis de funcionamiento de las obras proyectadas. En efecto, el análisis detallado de los resultados en lo que se refiere a la estructura del flujo se resume a través de un desplazamiento del núcleo de máximas velocidades y la formación de un patrón de flujo de recirculación entre espigones (marcado en E4-E5 y E7), con la presencia de estructuras en planta en forma (aproximada) de elipses elongadas, con velocidades máximas del orden de 1 m/s y menores para el escenario simulado.

Estas estructuras del flujo, en especial para los caudales ordinarios (de recurrencias menores a las de diseño de las obras) propiciarán los procesos de sedimentación en los recintos comprendidos entre los espigones, fundamentalmente teniendo en cuenta la abundante carga sedimentaria transportada por este río (Farias, 2017). La reducción del ancho de la barra empuntada, para la condición con espigones, permite el desplazamiento hacia la región interna de la curva de la sección hidráulica estable, en especial en espigones E1 a E4. Esto podría deberse a que son los espigones cuya disposición geométrica les permite depositar sedimento en el primer evento de crecidas.

Se puede concluir que la disposición geométrica (‘layout’) del campo de espigones que se ha proyectado resultó adecuada para las condiciones hidráulicas simuladas, que se consideran representativas del comportamiento de este tipo de cauces para el escenario hidrológico analizado.

Nomenclatura

x, y (m)	=	componentes del sistema cartesiano
Z (m)	=	elevación de superficie libre
H (m)	=	profundidad de agua
g (m/s2)	=	aceleración de gravedad
u, v (m/s)	=	componente de velocidad en coordenadas x e y respectivamente
T (ºC, g/l)	=	trazador pasivo
νf (m2/s)	=	coeficiente de difusión turbulenta del flujo
νT (m2/s)	=	coeficiente de difusión turbulenta del trazador
t (s)	=	tiempo
Sh (m/s)	=	fuente o sumidero de fluido
Sx, Sy (m/s2)	=	términos fuentes que representa el viento, fuerza Coriolis, fricción de fondo o una fuente o sumidero de momento en el dominio
STr	=	(unidades del trazador) fuentes o sumideros del trazador

Agradecimientos

El autor Lucas Dominguez Ruben agradece el apoyo a la Universidad Nacional del Litoral por haber permitido el uso de los recursos computacionales para la puesta en marcha de las simulaciones numéricas. Parte de las investigaciones aquí presentadas, se enmarcan en el Proyecto 23/C171, y fueron financiadas por el CICYT-UNSE, Argentina.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.