Método FOSM

La cuantificación de la incertidumbre se puede expresar a través del coeficiente de variación (ecuación (1)), que relaciona la desviación estándar y el valor esperado de la función de desempeño. En este estudio, la función de desempeño corresponde a la profundidad, velocidad e intensidad del flujo (variables de las ecuaciones de aguas someras), y la variable de entrada para la que se hace el análisis de propagación de incertidumbre es el coeficiente de rugosidad de Manning. Existen diferentes métodos para calcular los momentos estadísticos de la función de desempeño; por ejemplo: estimación puntual, primer orden segundo momento, segundo orden segundo momento, método de Hasofer-Lind, simulaciones de Montercarlo (^{Baecher y Christian, 2003}).

El método empleado en este trabajo es primer orden segundo momento (first order second moment [FOSM]). Para la estimación del primer y segundo momento estadístico (el valor esperado E[g] y la varianza Var[g]) se usan los dos primeros términos de la expansión de la serie de Taylor de una función de desempeño g con n variables aleatorias X _n (ecuación (3)). En FOSM se supone una diferencia pequeña entre la variable X _i y el valor promedio de la misma µX _i , por lo que el cuadrado o las potencias de orden mayor son más pequeñas y pueden ser ignoradas (^{Arévalo-Mendoza et al., 2014}; ^{Baecher y Christian, 2003}).

Para encontrar los momentos estadísticos de la función de desempeño, es necesario conocer información de la función de densidad de probabilidad de las variables analizadas. FOSM permite calcular la esperanza matemática y la varianza de la función de desempeño, sin conocer la función de densidad de probabilidad de las variables de entrada. Estos momentos estadísticos se calculan alrededor del valor medio y la desviación estándar de las variables analizadas. Las ecuaciones (3) y (4) se obtienen de la definición general de E [g] y Var[g], suponiendo que las variables de entrada son independientes entre sí (no hay correlación) y aplicando la expansión de la serie de Taylor para la función de desempeño (^{Baecher y Christian, 2003}; ^{Maskey y Guinot, 2003}).

En la ecuación (4), las derivadas parciales de la función de desempeño son evaluadas en los valores medios de las correspondientes variables de entrada (m = µX ₁ , µX ₂ , . . . , µX _n ) y σ ² es la varianza de cada una de estas variables. Si hay n variables de entrada, es necesario evaluar n derivadas parciales de la función de desempeño.

Modelo numérico de inundación

La función de desempeño corresponde a las variables de salida del modelo numérico de aguas someras implementado en HEC-RAS 2D, solucionadas con el método de volúmenes finitos (^{Brunner, 2016}). Estas ecuaciones son una simplificación de las ecuaciones de Navier-Stokes que describen el movimiento de un fluido en el espacio (x, y, z). Las ecuaciones (5) y (6) corresponden a la ecuación de conservación de masa y las ecuaciones de conservación de momentum, respectivamente, que conforman las ecuaciones de aguas someras bidimensionales.

Donde t es el tiempo; H es la elevación de la superficie de agua; ∇ es el vector de derivadas parciales (∂/∂x, ∂/∂y); h es la profundidad del flujo; V = (u, v) es la velocidad vectorial; a la vez, u es el componente horizontal y v el componente vertical; q son fuentes o sumideros de masa; g corresponde a la aceleración de la gravedad; v _t es el coeficiente horizontal de la viscosidad turbulenta; c _f es el coeficiente de fricción de fondo; f es el parámetro de Coriolis, y k es un vector unitario en la dirección vertical.

HEC-RAS 2D implementa la aproximación de Gauckler-Manning-Strickler (^{Cortés-Zambrano et al., 2022}) para calcular el coeficiente de fricción de fondo c _f (ecuación (7)). Al reemplazarlo en la ecuación de momentum (6), y a partir de las consideraciones para onda difusiva se obtiene una expresión en función del coeficiente empírico de rugosidad de Manning (n) y del radio hidráulico (R) (ecuación (8)).

El índice de intensidad de flujo (I _DF ) no es una variable del modelo numérico de aguas someras, sin embargo, su cálculo depende de la profundidad y la velocidad producto de las simulaciones hidrodinámicas. El I _DF y la magnitud de la velocidad al cuadrado (ecuación (9)); los valores para I _DF se toman como estimaciones de la severidad del flujo (véase^{Altarejos-García et al., 2012}).

Recientemente, ^{Ramos et al. (2021)} plantearon el I _DF para la evaluación de amenaza por avenidas torrenciales en Colombia así: si el I _DF ≥ 50, la probabilidad de colapso es alta y corresponde a amenaza alta. Si 1 ≤ I _DF < 50, la probabilidad de colapso es baja pero la probabilidad de algún daño estructural es considerable, por lo que se asigna una categoría media de amenaza. Por último, si 0 ≤ I _DF < 1, la probabilidad de daño estructural es baja; por tanto, corresponde a amenaza baja.

Incertidumbre del coeficiente de rugosidad de Manning

Para determinar la variabilidad natural del coeficiente de rugosidad de Manning se debe contar con información espacial y temporal detallada de la zona de estudio (^{Lumbroso y Gaume, 2012}; ^{Pappenberger et al., 2005}; ^{Trieste y Jarrett, 1987}). En las simulaciones hidrodinámicas para inundaciones es común que el valor de este coeficiente sea usado para calibrar las manchas de inundación y que el valor supuesto corresponda a un valor representativo de toda la zona de estudio. En la tabla 1 se presenta el valor medio y el coeficiente de variación (COV) para los coeficientes de rugosidad de Manning reportados por ^{Chang et al. (1993)}; ^{Da Costa et al. (2018)}; ^{Dingman y Sharma (1997)}, y ^{Kim et al. (2010)}.

Los valores presentados en la tabla 1 reportados por ^{Dingman y Sharma (1997)} fueron calculados a partir de una muestra total de 621 datos obtenidos de múltiples estaciones en Nueva Zelanda y Estados Unidos. Por consiguiente, los COV de Dingman y Sharma (1997) corresponden a la variación de la muestra total (sin consideración de su ubicación espacial), lo que justifica la diferencia en el orden de magnitud con los otros COV resumidos en la tabla 1. Para los otros autores, el COV de Manning varía entre 1 % y el 10 % aproximadamente.

Tabla 1 Valores reportados en la literatura del valor medio y COV para el Manning 

1Las unidades del coeficiente de rugosidad están dadas por la constante de conversión; para SI k _n = 1 s/m ^1/3

Caso de estudio

^{Prada-Sarmiento et al. (2019)} reportaron que los principales flujos de detritos, flujos de lodos y flujos hiperconcentrados del evento ocurrido el 31 marzo de 2017 se transportaron por el río Sangoyaco y por las quebradas Taruca y Taruquita, afectando la zona nororiental del casco urbano de Mocoa (Putumayo). Luego de estas observaciones y de la delimitación hidrológica para estos cauces, se identificaron catorce subcuencas (microcuencas) que formaron parte de la iniciación y el transporte del flujo. La delimitación del dominio computacional para la simulación hidrodinámica incluyó la zona de depósito, la quebrada San Antonio y parte del río Mocoa (desembocadura de los cauces ya mencionados) (figura 1).

Fuente: Ramos et al. (2020) y ^{Prada-Sarmiento et al. (2019)}

Figura 1 Localización de Mocoa (Colombia) y generalidades de la zona de estudio: microcuencas, dominio computacional y cauces de interés 

Simulación hidrodinámica

Para la simulación en HEC-RAS 2D se usó un modelo digital de elevación DEM. Se generó una malla de 5 m de tamaño de celda para el dominio computacional de la figura 1. En cada punto de cierre de las microcuencas se cargó una hidrógrafa de caudal líquido como condición de borde interna dentro del dominio computacional (fuente de masa). Las hidrógrafas de caudal fueron tomadas del cálculo realizado por ^{Ramos et al. (2021)}. En la figura 2A se muestran las catorce hidrógrafas que se cargaron para un periodo de retorno de cien años. En la figura 2b se ilustra la hidrógrafa con el caudal pico para cada periodo de retorno. Se evidencia la diferencia entre los tres periodos de retorno usados; para el de cien años el caudal de entrada en la microcuenca 2 (79,39 m3/s) es aproximadamente el doble que el del periodo de retorno de 5 años (42,05 m3/s).

Fuente: ^{Ramos et al. (2021)} y Unidad Nacional para la Gestión del Riesgo de Desastres y Pontificia Universidad Javeriana (2018).

Figura 2 Hidrógrafas de entrada del modelo fluidodinámico para una duración de tormenta de tres horas a) Hidrógrafas para las catorce microcuencas de la zona de estudio para un Tr = 100 años. b) Comparación de las hidrógrafas de la microcuenca 2 para los tres periodos de retorno simulados (5, 25 y 100 años). 

Cálculo de la propagación de la incertidumbre

El cálculo de los COV para cada variable de salida de las simulaciones hidrodinámicas requiere del primer y segundo momento estadístico. En el segundo se debe evaluar las derivadas parciales de la función de desempeño para cada variable de entrada (ecuación (4)). Cuando la función de desempeño no tiene forma matemática explícita, es necesario implementar un método numérico para evaluar las derivadas parciales (^{Arévalo-Mendoza et al., 2014}). En este estudio se implementó el método de diferencias finitas centrales para calcular la derivada parcial de la profundidad, la velocidad y la intensidad del flujo respecto al coeficiente de rugosidad de Manning. Teniendo en cuenta que las ecuaciones de aguas someras no tienen una forma matemática explícita, la derivada parcial de la ecuación de la varianza (ecuación (4)) en función del coeficiente de Manning se calcula como se muestra en la ecuación (10); donde g(µn ± ∆ _n ) es cualquiera de las variables de salida evaluadas en µn + ∆ _n . La magnitud de ∆ _n deber ser lo suficientemente pequeña para encontrar la variación en dicho intervalo. Para este estudio se determinó ∆ _n = σ _n (desviación estándar del coeficiente de rugosidad).

Los valores que se suponen en este estudio para µn y σ _n se relacionan con los valores medios y los COV presentados en la tabla 1. Se propone un valor medio del coeficiente de Manning µn representativo (µn = 0.12) para toda la zona de estudio teniendo en cuenta las similitudes del lecho del río analizado por ^{Kim et al. (2010)} y el lecho del río Sangoyaco y las quebradas Taruca y Taruquita. Para considerar el impacto de la variabilidad natural del coeficiente de Manning, se analiza la incertidumbre para varios COV que están en el rango de la tabla 1: COV = 2 %, 5 %, 10 %, 15 % y 20 %.

Las ecuaciones (11), (12) y (14) resumen los cálculos que se realizaron para cada variable de salida del modelo. Estas ecuaciones se aplicaron en cada pixel, resultado de las simulaciones hidrodinámicas mediante una calculadora ráster.

Propagación de la incertidumbre y evaluación de la amenaza en puntos específicos

Se realizó un análisis de la propagación de la incertidumbre en puntos de la mancha de inundación, seleccionados de tal forma que existiera información para todos los periodos de retorno analizados. Los seleccionados fueron: el punto A, en la confluencia de las quebradas Taruca y Taruquita, después del ingreso de la hidrógrafa. El punto B, en la parte baja de la quebrada San Antonio, antes de su desembocadura en el río Mocoa. El punto C, en la intersección con la vía y la quebrada Taruca que se desvía hacia el sector de Pueblo viejo. El punto D, en la parte alta del río Sangoyaco, después de la entrada de la penúltima hidrógrafa de esta quebrada. El punto E, en la quebrada Sangoyaco antes de desembocar en el río Mocoa y de la última hidrógrafa. El punto F, entre los puntos A y C, en la zona de mayor dispersión de la mancha de inundación (figura 5).

Figura 5 Izquierda: ubicación de los puntos específicos en los que se analiza la propagación de la incertidumbre. Derecha: puntos sobre secciones transversales del dominio ST1, ST2 y ST3 

La propagación de la incertidumbre también se analizó a través de tres puntos ubicados perpendicularmente a la dirección del flujo (análogos a secciones transversales [ST]), así: (1) en la zona de mayor dispersión sobre la quebrada Taruquita entre los puntos A y F; (2) en la zona de menor dispersión del flujo sobre la quebrada Taruca y después de la intersección de esta quebrada con la vía, unos metros abajo del punto C; por último (3) unos metros después del punto D en la quebrada Sangoyaco. Las ST se asumieron de acuerdo con tres puntos: margen izquierda (I), centro (C) y margen derecha (D). La imagen de la derecha de la figura 5 muestra los perfiles de cada una de las ST y la ubicación de los puntos I, C, D dentro de cada ST.

La figura 6 muestra la relación entre el COV _h , COV _v , COV _I DF con el COV _n para cada uno de los puntos escogidos (A, B, C, D, E y F), donde se evidencia que a mayor COV _n hay mayor incertidumbre en la profundidad, la velocidad y el I _DF . Sin embargo, esta relación no es proporcional y los aumentos reales de incertidumbre se deben observar a través de la línea de referencia (línea punteada): para todos los puntos hay al menos un resultado que supera dicha línea, excepto el punto B. En la figura 6 también se puede identificar que la velocidad es la única variable que no tiene aumento real de la incertidumbre en ninguno de los puntos analizados, mientras que el comportamiento de la profundidad y del I _DF varían. Tampoco se identifica un comportamiento de la incertidumbre relacionado con los periodos de retorno para cada variable.

Nota: La línea de referencia punteada corresponde a la relación 1:1 entre las ordenadas y las abscisas.

Figura 6 Valor del COV _h , COV _v , COV _I DF en función del COV _n para cada uno de los puntos seleccionados de la mancha de inundación (A, B, C, D, E y F) 

Los resultados mostrados en la figura 6 para los puntos C y E son muy similares, mientras que los resultados de cada variable en el punto D son muy parecidos entre periodos de retorno; en todo caso, estos tres puntos son los que muestran menor dispersión comparados con los puntos A, B y F. Según las figuras 1 y 5, los puntos B y C coinciden con puntos de confluencia en zonas de previa dispersión de la mancha de inundación, pero de diferentes quebradas; mientras que los puntos D y E corresponden a una misma quebrada en una zona uniforme de la mancha de inundación. Por otro lado, los puntos A y F están ubicados en una zona de alta dispersión de la mancha de inundación, la dispersión de los resultados en estos puntos coincide: tanto la profundidad como el I _DF para Tr = 5 y 25 años tienen incertidumbres ubicadas por encima de la línea de referencia.

En cuanto a los resultados para las ST, la ubicación de los puntos se condicionó de tal forma que existiera información en todos los escenarios, i.e. para los tres periodos de retorno y las tres variables analizadas. La figura 7 muestra los valores de los COV _h , COV _v , COV _I DF en función del COV _n para los tres periodos de retorno graficados por cada punto dentro de las ST. De las tres ST de la figura 7, la ST N. º3 es la de menor incertidumbre (p. ej.: no se sobrepasa la línea de referencia); sin embargo, en esta ST no es sencillo identificar un comportamiento particular para cada variable como sí es posible identificar en las otras ST; en la ST 1, por ejemplo, para el margen izquierdo la profundidad tiene una incertidumbre entre 0-8 % con una relación directa para los periodos de retorno (a mayor Tr, mayor COV ), la velocidad entre 0-16 % con una relación inversa para los periodos de retorno y el I _DF entre 1-25 %, también con una relación inversa para los periodos de retorno. La ST 2, margen derecho y centro, coincide con la ST 1 en cuanto a los rangos, pero no para la relación de los periodos de retorno: para la profundidad la relación es inversa, mientras que, para la velocidad y el I _DF la relación es directa. De acuerdo con la ST 1 y ST 2 de la figura 7, la variable con mayor incertidumbre es el I _DF , solo el margen derecho de la ST 1 indica que la profundidad tiene mayor incertidumbre.

Nota: (primera fila), n.º 2 (segunda fila) y n.º 3 (tercera fila); extremo izquierdo (columna izquierda), punto central (columna central) y extremo derecho (columna derecha). La línea de referencia punteada corresponde a la relación 1:1 entre las ordenadas y las abscisas.

Figura 7 Resultados del COV _h , COV _v , COV _I DF en función del COV _n para cada uno de los puntos de las secciones transversales analizadas 

^{Altarejos-García et al. (2012)} exponen mapas de desviación estándar y del valor esperado para la profundidad y la velocidad que son resultados de una simulación bidimensional en el software GUAD-2-D (ecuaciones de aguas someras-volumen finito) para un tramo de corriente con tres caudales picos. En un cálculo del COV al relacionar los valores máximos de cada uno de los mapas presentados, se obtiene un 0 ≤ COV _h ≤ 18 % y 0 ≤ COV _v ≤ 16 %; el rango de la velocidad coincide con los resultados generalizados para la figura 6.

La incertidumbre en la evaluación de la amenaza es menor en las partes altas de las cuencas, mientras que en las zonas donde las pendientes son menores, la incertidumbre en la altura, velocidad e intensidad de flujo es mayor. Lo anterior tiene implicaciones serias en el ordenamiento del territorio, dado que típicamente las comunidades habitan en las partes bajas de las cuencas donde, a pesar de que las velocidades y alturas de flujo pueden ser menores, su comportamiento es menos predecible y por ende su incertidumbre aumenta.

Los resultados evidencian que al evaluar la amenaza por inundación, teniendo en cuenta únicamente la incertidumbre o variabilidad natural temporal asociada a los periodos de retorno, se subestima la cuantificación de la amenaza debida a la variabilidad natural espacial de aspectos relacionadas con las propiedades físicas del flujo.

No se observa un incremento significativo de la incertidumbre de las alturas, velocidades o índice de intensidad de flujo cuando aumenta el caudal representado por escenarios de diferentes periodos de retorno; sin embargo, la evaluación de la magnitud del peligro es mayor debido a que los valores esperados de las alturas, velocidades e índice de intensidad de flujo sí se incrementan con los escenarios de periodos de retorno.

Para la evaluación de la amenaza por cualquier peligro natural (inundaciones, avenidas torrenciales), se usa la curva de amenaza que relaciona en las abscisas los valores de periodo de retorno y en las ordenadas se representa la variable que caracteriza la amenaza (altura, velocidades o índices de intensidad de flujo). En la figura 8A se presenta la curva de amenaza para el I _DF de los puntos F y C de la figura 5. Es claro que la curva del punto C tiene un nivel de amenaza mayor que la curva F y que la incertidumbre natural está dada únicamente por el periodo de retorno, es decir, variabilidad natural temporal. Con el procedimiento acá presentado, se obtiene además la posibilidad de calcular la probabilidad que el I _DF sea mayor a un valor determinado (asociado a un nivel de amenaza alta, media o baja I _{DF −Amenaza} ) dado un periodo de retorno (Tr) P (I _DF ≤ I _{DF −Amenaza} )/Tr, debido a que la relación periodo de retorno e I _DF no es determinística, sino que existe variabilidad natural espacial debido a los parámetros físicos que se encuentran en el modelo de inundación. La incertidumbre asignada a cada valor del I _DF dado un Tr se obtiene de lo reportado en la figura 6. En la figura 8b se presenta la función de distribución de probabilidad del I _DF para Tr 5, 25 y 100 años del punto C. Se observa claramente que la incertidumbre es mayor en la determinación del I _DF a medida que aumenta el periodo de retorno de las hidrógrafas. Para el punto C, se tiene lo siguiente:

Figura 8 a. Curva de amenaza contra inundación de los puntos C y F en términos del IDF. b. Función de distribución de probabilidad del IDF para tres periodos de retorno (5, 25 y 100 años) en el punto C de la figura 5 

P [(I _DF ≤ I _DF = 1 m3/s2)/Tr = 5 años] = 0, 39; P [ (I _DF ≤ I _DF = 1m3/s2)/Tr = 25 años] = 0, 14; P [ (I _DF ≤ I _DF = 1 m3/s2)/Tr = 100 años ] = 0, 0069. Según ^{Ramos et al. (2021)}, un I _DF ≤ 1 m3/s2 se encuentra asociado a un nivel de amenaza baja. Por lo tanto, en el punto C de análisis y a manera de ejemplo, se encuentra que a medida que el periodo de retorno aumenta, la probabilidad de estar en un nivel de amenaza baja disminuye. Consistentemente, la probabilidad de estar en un nivel de amenaza medio o alto es superior con el aumento del periodo de retorno de las hidrógrafas.

En resumen, la figura 8A ilustra la evaluación de la amenaza por inundación en términos del IDF, considerando únicamente la incertidumbre natural temporal dada por la variación de la magnitud de la amenaza en función del periodo de retorno, mientras que la figura 8B muestra la influencia de la incertidumbre natural espacial en los resultados de la evaluación de la magnitud de la amenaza.

La integración de las dos incertidumbres naturales (espacial y temporal) y la evaluación de la amenaza integrada se efectuarían a través de la ecuación (14):

Donde las probabilidades P (I _DF i/TR _i ) se obtienen con el procedimiento presentado en este artículo para los n periodos de retorno y P (TR _i ) están asociados a la probabilidad de excedencia de las hidrógrafas reportadas en la figura 2 .