Simulación de Avenidas: Modelos 1D y 2D
Simulación de Avenidas: Modelos 1D y 2D
Facultad de Ingeniería
Centro Interamericano de Recursos del Agua
Título de la Tesis:
Alumno:
Ing. Juan Antonio Morales Orozco
Director de Tesis
Co-Director de Tesis
Agradezco el apoyo económico otorgado para la realización de esta tesis para obtener el título de
Maestro en Ciencias del Agua: Al Consejo Nacional de Ciencia y Tecnología (CONACyT), dentro
del proyecto de investigación 3952/2015E “Modelación de eventos hidrológicos extremos a partir
de precipitación estimada por percepción remota”, en el marco de la convocatori “Proyectos de
Desarrollo Científico para Atender Problemas Nacionales, 2014”. Así como del proyecto
“Hidrología en Mauritania: Modelación de caudales a partir de precipitación estimaada por
imágenes de satélite” Clave 4192/2016E.
Al Dr. Carlos Diaz Delgado y al Dr. Humberto Salinas Tapia, por su apoyo intelectual y moral para
el avance y mejora del proyecto de tesis que se presenta.
A los Profesores que durante éstos años de estudio de maestría han sido una fuente abundante de
conocimiento y motivación.
A los miembros integrantes del jurado, por sus acertadas recomendaciones y aportaciones:
2
Dedicatoria
Quiero dedicar este trabajo, resultado de años de trabajo y esfuerzo a todas las personas que de alguna
manera han apoyado su desarrollo, bien con aporte de conocimiento, bien con apoyo moral que me
han motivado a continuar.
A mi madre, que siempre tiene tiempo para darme ánimos y un jalón de orejas cuando es necesario.
Y más que nada, que siempre ha confiado en mí, y en todos los proyectos que he emprendido.
A Richard, Karen, Rebeca, Cata y Christian, que en el proceso de esta tesis han permitido que
conserve la cordura.
3
Resumen
La correcta simulación de avenidas es una herramienta importante para analizar el impacto que un
evento de precipitación puede tener sobre una zona . Sobre todo si dicha zona tiene actividad
humana. Por lo tanto se requiere delimitar cuidadosamente la extensión de la zona afectada para
determinar apropiadamente el alcance del impacto de las inundaciones sobre los asentamientos
humanos.
Los modelos de simulación 1D utilizan secciones transversales a lo largo del eje del río como
unidades de cálculo, y resolviendo las ecuaciones de continuidad y de
Para tal labor, contamos con modelos de simulación en una dirección (1D) y dos direcciones (2D),
cada uno con su s ventajas y desventajas. Los modelos 1D por un lado requieren una menor cantidad
de datos y poco poder computacional para funcionar, por lo que requieren una menor cantidad de
tiempo para proporcionar resultados. Los modelos 2D, por otro lado utilizan una mayor cantidad de
datos y poder computacional, a cambio de de poder modelar cambios en las direcciones del flujo que
están fuera del alcance de los modelos 1D.
Existen situaciones en las que ambos modelos tienen poca o ninguna diferencia en sus resultados,
principalmente cuando la avenida presenta una dirección de flujo predominante, es decir, cuando el
flujo está confinado o bien, la dirección del flujo de la planicie de inundación es la misma que la del
rio principal, los resultados no varían sensiblemente entre ambos modelos y por tanto puede optarse
por utilizar el modelo 1D sobre el 2D por su eficiencia computacional.
Si ambos modelos son combinados en un solo esquema de manera adecuada, es posible utilizar el
modelo 1D en las zonas donde el flujo está confinado y el 2D en planicies de inundación con más de
una dirección predominante, optimizando así los tiempos de computación.
Para que las simulasiones puedan llevarse a cabo, se requieren como datos mínimos a) datos
batimétricos y b) datos de uso de suelo. Adicionalmente, para la simulación de eventos específicos,
son necesarios c) datos hidrométricos d) datos fisiográficos de la cuenca.
Los datos batimétricos fueron obtenidos como parte del Continuo de Elevaciones Mexicano Versión
3 (CEM V.3).
i
Los datos batimétricos del CEM V.3 tienen una resulución de 15 m, aunque las elevaciones en
cuerpos de aguas profundas no lograron ser registradas, debido a que en general, los levantamientos
LiDAR no logran penetrar la superficie de agua a demasiada profundidad.
Es por esto que para la correcta simulación hidrodinámica, fue necesario determinar un fondo
sintético para los cuerpos de agua que en el CEM V.3 se marcan como valores nulos. Utilizamos el
método de Corum para determinar un fondo sintético que sea equivalente en cuanto funcionamiento
hidráulico . El método de Corum para obtener lechos sintéticos de rio requiere conocer la elevación
del espejo de agua y el caudal relacionado con esta elevación. Se realizan simulaciones sucesivas en
HEC-RAS usando el caudal conocido, restando en cada una la diferencia de la elevación del espejo
de agua simulado y el espejo de agua conocido al fondo de las secciones del canal, hasta que la
diferencia entre ambos sea menor que una tolerancia predeterminada (0.006 m, en el caso de la
presente investigación).
Los caudales conocidos son los correspondientes a la salida de la cuenca. No es apropiado considerar
estos caudales constantes a lo largo de todo el eje del río, de forma que es distribuido a lo largo del
eje utilizando el modelo y software CEQUEau, por medio de estaciones ficticias. Después que se
determina el caudal aportado por cada estación ficticia en un periodo de tiempo, se utiliza el software
de HEC-RAS para determinar el tiempo de escurrimiento desde cada estación ficticia hacia la salida
de la cuenca, las herramientas del Integrated Water Management son utilizadas para determinar las
propiedades fisiográficas (área y tiempo de concenctración nos serán indispensables) de cada
subcuenca delimitada por las estaciones ficticias, y se usan estos datos para estimar un hidrograma
sintético utilizando el método de la SCS.
Los resultados obtenidos verifican que un esquema combinado de simulación puede incrementar la
velocidad de cálculo hasta por más de veinte veces, siguiendo las recomendaciones siguientes
ii
4. Las zonas de la planicie de inundación cuyo sentido de flujo sea diferente al del eje principal
del rio deben ser discretizados por una malla 2D y unidos lateralmente al modelo 1D del canal
principa.
5. Los meandros que para un periodo de retorno máximo predefinido por el proyectiasta no
presenten desbordamientos que generen varias intersecciones de la mancha de inundación
con el mismo eje del río, se representarán utilizando un modelo 1D cuidando que no exista
traslape entre las secciones trasnversales. De otra manera, se sustituirá el tramo de río 1D por
una malla 2D unida frontalmente aguas arriba y aguas abajo del meandro.
6. En zonas donde haya riesgo de ruptura de diques, es aconsejable unir lateralmente una malla
2D a las secciones transversales en las que se puede presentar la ruptura
7. Cuando existen canales adyacentes al rio principal (naturales o artificiales) es recomendable
modelar tales canales usando una malla 2D unidad lateralmente al modelo 1D.
8. Cuando entre secciones transversales existan zonas que tiendan a fluir en sentido contrario al
del rio principal, se unirá lateralmente una malla 2D al modelo 1D para que el vector de
dirección en esa zona pueda ser representado apropiadamente.
9. Las fronteras entre los modelos 1D y 2D deben tener un tamaño igual o menor al del pixel
del DEM, por un lado para asegurar que la dirección del flujo en la trancisión de los modelos
sea apropiada, y por otro lado para disminuir el número de iteraciones necesarias en la
transferencia de flujo entre los modelos 1D y 2D. Puede modificarse el tamaño de la malla
solo en la frontera mediante el uso de líneas de quiebre. Dado al uso de submalla que hace
HEC-RAS en el modelo 2D, la mancha de inundación no muestra cambios sensibles
En cada caso en el que sea necesario emplear una malla 2D, debe tenerse en consideración que en
cada iteración, todas las celdas presentan carga computacional, incluso si no hay flujo a través de
ellas (celdas secas), por lo es deseable que las mallas no se extiendan innecesariamente. La mejor
manera de determinar las zonas en las que es necesario realizar una combinación de modelos 1D y
2D, es realizar previamente una simulación puramente 1D, transitando los caudales correspondientes
a un periodo de retorno considerado adecuado por el proyectista para verificar las posibles áreas de
inundación.
iii
Contenido
Notaciónes utilizadas y unidades ...................................................................................................... xii
Introducción................................................................................................................................ 1
1.3 Justificación......................................................................................................................... 8
Marco teórico............................................................................................................................ 11
iv
3.2 Tratamiento de datos ......................................................................................................... 42
3.3 Generación de hidrogramas sintéticos por el método de la SCS distribuidos a lo largo del
eje del río principal ....................................................................................................................... 50
Aplicación................................................................................................................................. 63
v
6.1.3 EXP2TIF ................................................................................................................. 147
vi
Lista de figuras
FIGURA 1-1 ESQUEMA DE VECTORES DE FLUJO RESULTANTES BAJO EL EMPLEO DE UN MODELO 1D DONDE EXISTE UN OBSTÁCULO. 6
FIGURA 2-1 MÉTODO DEL PASO ESTÁNDAR PARA LA SOLUCIÓN DE LA ECUACIÓN DE LA ENERGÍA PARA FLUJOS ESTABLES EN 1D ... 12
FIGURA 2-2 ALGORITMO PARA SOLUCIONAR LAS ECUACIONES DE FLUJO NO ESTABLE (1D Y COMBINACIÓN 1D/2D) ................. 13
FIGURA 2-3 EJEMPLO DE DATOS BATIMÉTRICOS DENTRO DE LA MALLA DE CÁLCULO PARA LA SIMULACIÓN 2D .......................... 14
FIGURA 2-4 ESQUEMA DE ELEMENTOS DE CÁLCULO 2D ................................................................................................ 15
FIGURA 2-5 ALGORITMO DE SIMULACIÓN FLUJO NO ESTABLE 2D .................................................................................... 15
FIGURA 2-6 DIVISIÓN DE LA SECCIÓN TRASVERSAL PARA DETERMINAR LA VELOCIDAD PONDERADA (VER ECUACIÓN 2-6) ........... 17
FIGURA 2-7 IMPORTACIÓN DE GEOTIFF A MATLAB PARA EL USO EN TOPOTOOLBOX ......................................................... 24
FIGURA 2-8 ACUMULACIÓN DE UNIDADES DE FLUJO DENTRO DEL DEM (IMAGEN TOMADA DE SCHWANGHART Y KUHN, 2010) . 25
FIGURA 2-9 ALGORITMO DE CORUM PARA OBTENER UN LECHO SINTÉTICO ....................................................................... 27
FIGURA 2-10 DIAGRAMA ADIMENSIONAL DE LA SCS, ADAPTADO POR CHOW, (2000) DEL SOIL CONSERVATION SERVICE
ENGINEERING HANDBOOK, 1972. .................................................................................................................. 37
FIGURA 3-1 ESQUEMA DE MODELO COMBINADO 1D/2D .............................................................................................. 39
FIGURA 3-2 METODOLOGÍA DE INVESTIGACIÓN PROPUESTA .......................................................................................... 40
FIGURA 3-3 ALGORITMO PARA CAMBIAR LOS VALORES NULOS DEL DEM POR LOS CORRESPONDIENTES VALORES DE LA ELEVACIÓN
DE LA SUPERFICIE DE AGUA AL MOMENTO DEL LEVANTAMIENTO LIDAR .................................................................. 44
FIGURA 3-12 ALGORITMO GENERAL PARA PREPARAR LOS DATOS DE USO DE SUELO PARA SU INCORPORACIÓN AL MODELO DE
SIMULACIÓN 1D DE HEC-RAS ....................................................................................................................... 50
FIGURA 3-13 EJEMPLO DE PROCESO DE OBTENCIÓN DE ÁREA DE APORTACIÓN EXCLUSIVA PARA UNA ESTACIÓN FICTICIA PARTICULAR
(REF. UTM-15N) ....................................................................................................................................... 52
FIGURA 3-14 ESQUEMA DE CÁLCULO DE HIDROGRAMAS SINTÉTICOS POR MÉTODO DE LA SCS CON REFERENCIA DE TIEMPO EN LA
DESEMBOCADURA ........................................................................................................................................ 53
vii
FIGURA 3-15 ALGORITMO PARA REALIZAR LAS SIMULACIONES HIDRODINÁMICAS INICIALES DE AVENIDAS EN 1D ...................... 55
FIGURA 3-16 INUNDACIÓN EN SECCIÓN CON CANALES ADYACENTES SIN PUNTOS DE DIQUE SEÑALADOS .................................. 55
FIGURA 3-17 A) DIQUES MAL COLOCADOS (POR INCLUIR LOS CANALES ADYACENTES) B) DIQUES BIEN COLOCADOS (POR AISLAR EL
CANAL PRINCIPAL DE LOS ADYACENTES) C) DIQUES MAL COLOCADOS (POR NO CONSIDERAR LOS CANALES ADYACENTES EN
SITUACIÓN DE INUNDACIÓN) D) DIQUES BIEN COLOCADO (POR CONSIDERAR LOS CANALES ADYACENTES EN SITUACIÓN DE
INUNDACIÓN) .............................................................................................................................................. 56
FIGURA 3-18 EJEMPLO DE MALLA 2D DEMASIADO AMPLIA RESPECTO A LA PLANICIE DE INUNDACIÓN ESPERADA ...................... 60
FIGURA 3-19 EJEMPLO DE MALLA 2D DEMASIADO PEQUEÑA RESPECTO A LA PLANICIE DE INUNDACIÓN ESPERADA .................... 60
FIGURA 3-20 A) EJEMPLO DE PLANICIE DE INUNDACIÓN OBTENIDA CON UNA SIMULACIÓN 1D B) EJEMPLO DE PLANICIE DE
IUNDACIÓN OBTENIDA CON UNA SIMULACIÓN 2D ............................................................................................... 61
FIGURA 4-6 DIRECCIONES POSIBLES DEL SIGUIENTE PÍXEL DEL RÍO PRINCIPAL A PARTIR DEL PÍXEL ACTUAL (PA) ......................... 70
FIGURA 4-7 DIRECCIÓN DE SECCIONAMIENTO SEGÚN LA DIRECCIÓN DEL PÍXEL DEL EJE PRINCIPAL DEL RÍO................................ 70
FIGURA 4-8 ALGORITMO DE EXTRACCIÓN DE DATOS DE LOS PIXELES DEL EJE DEL RÍO ........................................................... 70
FIGURA 4-9 MÉTODO PARA CAMBIAR LOS VALORES NULOS POR LAS ELEVACIONES DE LOS ESPEJOS DE AGUA EN LAS SECCIONES
TRANSVERSALES ........................................................................................................................................... 72
FIGURA 4-10 EJEMPLOS DE RESULTADOS DE LA PRIMERA ITERACIÓN DE CIERRE DEL RÍO A ) RESULTADOS EN TRAMOS RECTOS B)
RESULTADOS EN MEANDROS .......................................................................................................................... 72
FIGURA 4-11 TRATAMIENTO DE LOS PIXELES CON VALOR NULO QUE NO SE PROCESARON DENTRO DE LAS SECCIONES ................ 73
FIGURA 4-12 RESULTADOS PRIMER ITERACIÓN DEL ESPEJO DE AGUA................................................................................ 75
FIGURA 4-13 ALGORITMO SECCIONADO SIN TRASLAPE Y SIN CORTAR ESPEJOS DE AGUA ....................................................... 75
FIGURA 4-14 SITUACIONES POSIBLES DEL ESPEJO DE AGUA DENTRO DE LAS SECCIONES TRANSVERSALES .................................. 76
FIGURA 4-15 VENTANA DE CREACIÓN DE SUPERFICIES EN CIVIL3D ................................................................................. 80
FIGURA 4-16 VENTANA DE IMPORTACIÓN DE PUNTOS HACIA LA SUPERFICIE TIN EN CIVIL3D .............................................. 81
FIGURA 4-17 EJEMPLO DE SCRIPT DE CIVIL3D PARA IMPORTAR LÍNEAS ........................................................................... 81
FIGURA 4-18 APLICACIÓN DE LAS LÍNEAS DE FRONTERA A LOS DATOS DE TRIANGULACIÓN TIN.............................................. 82
FIGURA 4-19 APLICACIÓN DE LAS LÍNEAS DE LÍMITE A LOS DATOS DE TRIANGULACIÓN TIN ................................................... 83
FIGURA 4-20 COMPARACIÓN DE LA PRIMERA ITERACIÓN DEL ESPEJO DE AGUA CON LA SEGUNDA .......................................... 84
FIGURA 4-21 ELEVACIÓN DE LA SUPERFICIE DE AGUA EN LA DESEMBOCADURA .................................................................. 85
viii
FIGURA 4-22 IDEALIZACIÓN DE LA SECCIÓN TRANSVERSAL EN LA DESEMBOCADURA ............................................................ 86
FIGURA 4-23 DISTRIBUCIÓN ESPACIAL DE LA CUENCA EN CEQUEA (ADAPTADA DE MERCADO, V., 2010) ............................ 88
FIGURA 4-24 POSICIÓN INICIAL DE LOS PUNTOS DE BANCO Y DIQUES ............................................................................... 91
FIGURA 4-25 EJEMPLO DE SECCIÓN CON LECHO SINTÉTICO ............................................................................................ 91
FIGURA 4-26 HISTOGRAMA DE CAUDALES MÁXIMOS ANUALES DE LA ESTACIÓN 30016 ...................................................... 93
FIGURA 4-27 DISTRIBUCIÓN DE LOS DATOS OBSERVADOS EN COMPARACIÓN DE LOS VALORES ESPERADOS (PTE. 1) .................. 98
FIGURA 4-28 DISTRIBUCIÓN DE LOS DATOS OBSERVADOS EN COMPARACIÓN DE LOS VALORES ESPERADOS (PTE. 2) .................. 99
FIGURA 4-29 SUBCUENCAS POR ESTACIÓN FICTICIA .................................................................................................... 103
FIGURA 4-30 SUBCUENCA DE LA ESTACIÓN FICTICIA 8, CON DETALLE DEL RÍO PRINCIPAL DE LA SUBCUENCA Y DE LA CUENCA
GENERAL .................................................................................................................................................. 109
FIGURA 4-31 HIDROGRAMA SINTÉTICO DE LA ESTACIÓN FICTICIA 1, PARA TR=5 AÑOS ...................................................... 111
FIGURA 4-32 HIDROGRAMAS POR ESTACIÓN FICTICIA Y SU CONVOLUCIÓN PARA TR=5, CORRESPONDIENTES A UN DÍA DE
PRECIPITACIÓN. ......................................................................................................................................... 111
FIGURA 4-33 CONVOLUCIONES DE LOS HIDROGRAMAS DE LAS ESTACIONES FICTICIAS POR PERIODO DE RETORNO ................... 112
FIGURA 4-34 EJEMPLO DE HIDROGRAMA PARA TORMENTAS AISLADAS OBSERVADAS......................................................... 113
FIGURA 4-35 EJEMPLO DE ARCHIVO DE TEXTO PARA LA CONVERSIÓN DE LA MATRIZ DE CLAVES DE USO DE SUELO POR SU
RESPECTIVO VALOR N .................................................................................................................................. 114
FIGURA 4-44 CAUDAL A TRAVÉS DE LA DESEMBOCADURA DE LA CUENCA, CALCULADOS CON MODELOS 1D Y 2D PARA TR=5 AÑOS
.............................................................................................................................................................. 124
FIGURA 4-45 CAUDAL A TRAVÉS DE LA DESEMBOCADURA DE LA CUENCA, CALCULADOS CON MODELOS 1D Y 2D PARA TR=10 AÑOS
.............................................................................................................................................................. 125
FIGURA 4-46 CAUDAL A TRAVÉS DE LA DESEMBOCADURA DE LA CUENCA, CALCULADOS CON MODELOS 1D Y 2D PARA TR=20 AÑOS
.............................................................................................................................................................. 125
FIGURA 4-47 CAUDAL A TRAVÉS DE LA DESEMBOCADURA DE LA CUENCA, CALCULADOS CON MODELOS 1D Y 2D PARA TR=25 AÑOS
.............................................................................................................................................................. 126
ix
FIGURA 4-48 CAUDAL A TRAVÉS DE LA DESEMBOCADURA DE LA CUENCA, CALCULADOS CON MODELOS 1D Y 2D PARA TR=50 AÑOS
.............................................................................................................................................................. 126
FIGURA 4-49 CAUDAL A TRAVÉS DE LA DESEMBOCADURA DE LA CUENCA, CALCULADOS CON MODELOS 1D Y 2D PARA TR=100
AÑOS ....................................................................................................................................................... 127
FIGURA 4-50 COMPARACIÓN DE LAS ELEVACIONES DEL ESPEJO DE AGUA EN EL EXUTORIO PARA MODELOS 1D Y 2D................ 128
FIGURA 4-51 EJEMPLO DE PLANICIE DE INUNDACIÓN DEBAJO DEL NIVEL DE DESBORDE DEL RÍO ........................................... 129
FIGURA 4-52 ERROR EN LA SIMULACIÓN PARA TERRENOS CON GRANDES PLANICIES DE INUNDACIÓN POR DEBAJO DEL NIVEL DE
DESBORDE DEL RÍO. A) PLANICIE DE INUNDACIÓN OBTENIDA POR SIMULACIÓN 1D B) PLANICIE DE INUNDACIÓN OBTENIDA
FIGURA 4-53 MEANDROS Y CANALES ADYACENTES A) SIMULACIÓN 1D A.1) SECCIÓN CON DOS CANALES INUNDADOS A.2)
SECCIÓN CON TRES CANALES INUNDADOS B) SIMULACIÓN 2D (SE UTILIZÓ EN LA SECCIÓN 62+675 UN QP=1000 M3/S)
.............................................................................................................................................................. 132
FIGURA 4-54 DESBORDAMIENTO DE UNA SOLA SECCIÓN TRANSVERSAL A) MODELACIÓN 1D B) MODELACIÓN 2D ................ 133
FIGURA 4-55 COMPARACIÓN DE INUNDACIONES PEQUEÑAS CON REINGRESO DE AGUA AL SISTEMA SIMULADOS CON MÉTODOS 1D Y
2D. A) SIMULACIÓN 1D B) SIMULACIÓN 2D .................................................................................................. 134
Lista de tablas
TABLA 1-1 VENTAJAS Y DESVENTAJAS DEL MODELADO 1D Y 2D (ADAPTADA DE PAREDES ET AL, 2014) .................................. 5
TABLA 2-1 VALORES TÍPICOS DE D (BRUNNER, G., 2016). ........................................................................................... 19
TABLA 2-2 VALORES TABULADOS DEL DIAGRAMA ADIMENSIONAL DE LA SCS ..................................................................... 38
TABLA 3-1 DATOS REQUERIDOS POR CADA ACTIVIDAD DE LA INVESTIGACIÓN ..................................................................... 41
TABLA 3-2 FUENTE DE OBTENCIÓN DE DATOS REQUERIDOS ............................................................................................ 41
TABLA 3-3 ELEMENTOS DE ENTRADA DE DATOS AL SISTEMA 1D DE SIMULACIÓN Y DATOS REQUERIDOS POR ÉSTOS ................... 58
TABLA 3-4 PARÁMETROS PRINCIPALES PARA LOS DIFERENTES TIPOS DE SIMULACIÓN ........................................................... 59
TABLA 4-1 PARÁMETROS FISIOGRÁFICOS DE LA CUENCA DE LA ESTACIÓN HIDROMÉTRICA PUEBLO NUEVO (30016) ................. 64
TABLA 4-2 MARCADORES DE VALORES ANÓMALOS, DESDE AGUAS ABAJO HASTA AGUAS ARRIBA DE LA ZONA DE LA CUENCA ....... 68
TABLA 4-3 VARIABLES USADAS EN EL CÓMPUTO DE LA BATIMETRÍA (PRIMER ITERACIÓN) ..................................................... 70
TABLA 4-4 VARIABLES USADOS EN EL CÓMPUTO DE LA BATIMETRÍA (SEGUNDA ITERACIÓN) .................................................. 78
TABLA 4-5 EQUIVALENCIA DE CUADROS CEQUEAU EN COORDENADAS UTM Y CADENAMIENTO DEL EJE DEL RÍO ..................... 89
TABLA 4-6 DISTRIBUCIÓN DEL ESCURRIMIENTO DEL DÍA 02/01/2008 A LO LARGO DEL EJE .................................................. 90
TABLA 4-7 CAUDALES MÁXIMOS ANUALES OBSERVADOS EN LA ESTACIÓN 30016 .............................................................. 92
TABLA 4-8 PROPIEDADES GENERALES DE LA SERIE DE CAUDALES MÁXIMOS ANUALES DE LA ESTACIÓN 30016 .......................... 93
x
TABLA 4-9 RESULTADOS DE PRUEBAS DE CALIDAD SOBRE LA SERIE DE QMAX ANUALES DE LA ESTACIÓN 30016........................ 94
TABLA 4-10 NÚMERO DE OBSERVACIONES POR RANGO ................................................................................................ 94
TABLA 4-11 PARÁMETROS DE LA MUESTRA DE QMAX ANUALES PARA LA DISTRIBUCIÓN GAMMA .......................................... 95
TABLA 4-12 PARÁMETROS DE LA MUESTRA DE QMAX ANUALES PARA LA DISTRIBUCIÓN EXPONENCIAL ................................... 95
TABLA 4-13 PARÁMETROS DE LA MUESTRA DE QMAX ANUALES PARA LA DISTRIBUCIÓN PEARSON III ..................................... 95
TABLA 4-14 PARÁMETROS DE LA MUESTRA DE QMAX ANUALES PARA LA DISTRIBUCIÓN LOGNORMAL .................................... 95
TABLA 4-15 PARÁMETROS MUESTRA DE QMAX ANUALES PARA LA DISTRIBUCIÓN LOGPEARSON III........................................ 96
TABLA 4-16 PARÁMETROS MUESTRA DE QMAX ANUALES PARA LA DISTRIBUCIÓN GUMBEL TIPO I ......................................... 96
TABLA 4-17 COMPARACIÓN DE QMAX OBSERVADO CONTRA QMAX ESTIMADO POR DIFERENTES FUNCIONES DE AJUSTE ............ 97
TABLA 4-18 VALORES 𝒙𝒊 − 𝒙𝒊𝟐 PARA LA PRUEBA DE EFICIENCIA DE AJUSTE POR RAÍZ DE ERROR CUADRÁTICO MEDIO ............ 100
TABLA 4-19 VALORES 𝒙𝒊 − 𝒙𝒊𝟐 PARA LA PRUEBA DE EFICIENCIA DE AJUSTE POR RAÍZ DE ERROR CUADRÁTICO MEDIO ............ 101
TABLA 4-20 APROXIMACIONES DE CAUDAL PARA DIFERENTES TIEMPOS DE RETORNO UTILIZANDO LA FUNCIÓN GAMMA COMO
MODELO PROBABILÍSTICO DE AJUSTE .............................................................................................................. 102
TABLA 4-21 CARACTERÍSTICAS DE LAS SUBCUENCAS CORRESPONDIENTES A LAS ESTACIONES FICTICIAS .................................. 104
TABLA 4-22 LÁMINA POR EXCESO DE PRECIPITACIÓN SOBRE LA CUENCA POR PERIODO DE RETORNO .................................... 105
TABLA 4-23 DISTRIBUCIÓN ACUMULADA DE LAS APORTACIONES DE LAS ESTACIONES FICTICIAS AL QMAX MEDIO DIARIO ESPERADO
POR PERIODO DE RETORNO .......................................................................................................................... 105
TABLA 4-24 CAUDALES MEDIOS DIARIOS APORTADOS Y LÁMINA ESCURRIDA POR CADA ESTACIÓN FICTICIA POR PERIODO DE
RETORNO (RESTANDO EL FLUJO BASE) (PTE. 1) ................................................................................................ 106
TABLA 4-25 CAUDALES MEDIOS DIARIOS APORTADOS Y LÁMINA ESCURRIDA POR CADA ESTACIÓN FICTICIA POR PERIODO DE
RETORNO (RESTANDO EL FLUJO BASE) (PTE. 2) ................................................................................................ 107
TABLA 4-26 TIEMPO DE CONCENTRACIÓN TOTAL PARA LAS ESTACIONES FICTICIAS ............................................................ 108
TABLA 4-27 FACTORES DE FORMA PARA EL HIDROGRAMA SINTÉTICO DE LA SCS POR ESTACIÓN FICTICIA ............................... 109
TABLA 4-28 HIDROGRAMA SINTÉTICO PARA UN TR=5 AÑOS DE LA ESTACIÓN FICTICIA 1 .................................................... 110
TABLA 4-29 CAUDALES PICO POR HIDROGRAMA COMPUESTO ...................................................................................... 112
TABLA 4-30 COEFICIENTES DE RUGOSIDAD DE MANNING PARA LAS CLAVES DE USO DE SUELO ............................................ 115
TABLA 4-31 CAUDALES PICO (EN M3/S) DISTRIBUIDOS POR ESTACIÓN FICTICIA Y PERIODO DE RETORNO ................................ 118
TABLA 4-32 DIFERENCIA DE FLUJO PROMEDIO EN LA DESEMBOCADURA ENTRE SIMULACIONES 1D Y 2D ............................... 127
TABLA 4-33 TIEMPOS DE PROCESAMIENTO DE LAS SIMULACIONES ................................................................................ 135
xi
Nomcenclatura
Notaciónes utilizadas y unidades
xii
tp: tiempo de retardo de escurrimiento, en horas
tr: Duración efectiva de la lluvia, en horas
u,v Velocidades del flujo en las direcciones cartesianas del plano de la malla
V Velocidades promedio (Gasto/Área)
vt Coeficiente de viscosidad de torbellinos
X,Y: Coordenada en X e Y del elemento de la matriz de matlab (Ver Figura 2-7).
XUTM,YUTM: Coordenadas X e Y en formato UTM
Y Tirante del agua en la sección transversal
Z Elevación en el fondo del canal principal
Z1 Pendiente de la margen izquierda
Z2 Pendiente de la margen derecha
Zi, Zj: Elevación de los pixeles i y j, respectivamente
β Coeficiente del momento tomado en cuenta para una distribución variante de la
velocidad en canales irregulares
Γ: Función gamma
ΔA Cambio de área hidráulica dentro del canal principal
ΔA Cambio de área hidráulica de la planicie de inundación
Δx Longitud del canal principal entre dos secciones transversales
Δx Distancia de flujo equivalente
Δx La longitud de la planicie de inundación entre dos secciones transversales
𝜇 : Media de la muestra
𝜎 : Desviación estándar de la muestra
Acrónimos
xiii
Introducción
Paredes, Pedrozo y Mejía (2014) mencionan que a nivel global, la amenaza natural con mayor número
de eventos incidentes en poblaciones humanas son las inundaciones. La tendencia en aumento de los
eventos de inundación puede inferirse de las estadísticas del Consejo Internacional de la Cruz Roja,
ya que en el último siglo (1919-2004) la Cruz Roja participó en 500 eventos de inundaciones a nivel
mundial, mientras que en el periodo comprendido entre el año 2004-2014, participó en 1,752 eventos
de inundación (Cannon y Schipper, 2014).
Tan solo en México, entre los años de 1980 a 1999 las inundaciones han causado alrededor de 10,000
muertes (Salas y Jimenez, 2004) y en términos económicos, en el año 2002 los daños por fenómenos
hidrometeorológicos (inundaciones, lluvias y ciclones) ascienden a los 10,544 mdp, (Bitrán et al,
2002). Para el año 2014, esta cifra sube a 27,962 mdp, correspondiendo al 84.9% del total de daños
por desastres (naturales y antrópicos) ocurridos en la república durante ese año (García, Méndez y
Sánchez, 2015). Hay que señalar de cualquier forma que de ese total, 24,133 mdp en afectaciones
fueron originados por el huracán Odile. En consecuencia, el conocimiento del comportamiento de
éstos fenómenos es información poderosa e indispensable al momento de tomar decisiones acerca de
las acciones a ejecutar para intentar controlarlas y/o mitigar sus efectos negativos.
De acuerdo con los datos que se tienen sobre el impacto económico de los fenómenos
hidrometeorológicos ocurridos en México en los últimos quince años puede apreciarse una clara
tendencia al alza (García, Méndez y Sánchez, 2015). Esto en el marco del cambio climático toma
mayor importancia, Olcina (2008) menciona la incertidumbre que existe acerca de la evolución del
comportamiento de las precipitaciones, puesto que los modelos actuales fallan en hacer predicciones
a nivel global. Al respecto Fenoglio et al (2015) señalan que la pérdida de fiabilidad de los resultados
de los modelos de pronóstico hidrológico, se limita a las series temporales de los fenómenos
climatológicos. Por otro lado, Olcina (2008) menciona que debido a la incertidumbre climatológica,
también pueden crearse nuevas zonas de riesgo, de forma que hay que considerar el concepto de
“escenarios futuros de riesgo”.
Los métodos de teledetección para determinar regímenes de inundación utilizan imágenes satelitales
disponibles a lo largo del tiempo, teniendo como desventaja la baja disponibilidad temporal de éstas
imágenes cuyo uso además está condicionado por la cobertura de nubes al momento de la toma de la
imagen. Los modelos hidrológicos por otro lado requieren de datos de la región de interés como
entrada. El modelo consiste en vincular el balance hídrico, los volúmenes de agua, topografía y la
batimetría de la región de estudio (Gibbs, Clark y Taylor, 2015).
Una vez que los parámetros del modelo de simulación hidrodinámica de la zona de interés han sido
determinados, es posible realizar una modelación hidrodinámica que permita simular el alcance de
una inundación en periodos de tiempo en los que no se tiene información geoespacial de la inundación
(Karim et al, 2011) o en escenarios de predicción de eventos.
Los modelos 1D tienen como ventaja la velocidad de cálculo, pero dejan de lado algunos aspectos
importantes de la variabilidad espacial, como la velocidad que es analizada en un sólo sentido
(Paredes, Pedrozo y Mejía, 2014) y el tratamiento de los flujos de inundación es muy simplificada
(Kuiry, Sen y Bates, 2010).
El análisis en dos dimensiones, comparado con el análisis en una dimensión requiere de una elevada
cantidad de información y poder de procesamiento, pero simula un comportamiento más fiel del flujo
en un río a través del tiempo (Morales, Petaccia y García, 2015).
Entre los modelos de simulación hidráulica de tránsito de avenidas se destaca el modelo HEC-RAS
(Brunner, 2016).
2
Hec-RAS es un software ampliamente usado a nivel mundial debido a que es gratuito y que cuenta
con una amplia documentación de aplicaciones. Corum (2015), utilizó el software de HEC-RAS
para modelar diferentes ríos del Noroeste del Pacífico (Los ríos Green River, Skykomish River,
Dungeness River en Whashington y Clark Fork entre Idaho y Montana, todos en Estados Unidos).
En algunas situaciones, donde el flujo permanece confinado al canal principal del río, los resultados
de las simulaciones usando modelos 1D y 2D son similares en cuanto a la extensión de la inundación
(Werneret al, 2005, Paredes et al, 2014). Cuando el flujo no está confinado, caso que se presenta
durante el desbordamiento del canal de un río, los modelos de sumulación hidrodinámica 2D ofrecen
mejores resultados. Es por esto, que el análisis más eficiente sería uno que combine el modelado de
flujo en 1D y 2D (Morales, Petaccia y García, 2015). De ahí, la importancia de desarrollar
estrategias y diseñar criterios que permitan seccionar tramos de un río para que sean analizados con
uno u otro modelo según corresponda. El presente proyecto de tesis busca proponer un conjunto de
éstos criterios.
1.1 Antecedentes
Los modelos climatológicos han comenzado a perder fiabilidad en cuanto a los análisis probabilísticos
de ocurrencia y magnitud de los eventos de precipitación debido, entre otras cosas, a que el cambio
climático que se está presentando en el planeta afecta las características de estacionalidad de los datos
hidroclimáticos que alimentan tales modelos (Olcina, 2008).
Entonces es posible señalar dos factores que aumentan la importancia de la simulación hidrodinámica
sobre el área de estudio. Por un lado, la posibilidad de que los escenarios de riesgo cambien (Olcina,
2008) y por otro, que la vulnerabilidad debida a la batimetría de la zona varía muy lentamente en el
tiempo (Fenoglio et al, 2015).
Así que una vez que el modelo hidrodinámico ha sido calibrado en una zona de interés, pueden
simularse una gran cantidad de situaciones de flujo sobre ésta, ya sean simulaciones por evento
(predicciones por periodo de retorno), simulaciones usando datos en tiempo real (percepción remota
y/o estaciones de monitoreo automáticas) o incluso determinar indicadores que señalen el “umbral de
inundación” sobre una zona de interés (Gibbs, Clarke y Taylor, 2015) .
3
Con la simulación de éstos escenarios de riesgo por inundación, es posible tomar mejores decisiones
sobre las medidas de prevención de inundaciones en zonas donde tradicionalmente éstos problemas
se tenían controlados (Fenoglio et al, 2015). Cabe señalar que el esfuerzo y costos requeridos en las
medidas de prevención y mitigación de inundaciones son demasiado elevados como para llevarse a
cabo sin un estudio previo de costo-beneficio (Baró et al, 2012).
Según Paredes et al (2014) los modelos presentan mayor certidumbre cuando el modelo se calibra
con el gasto medio, así que entre mayor es la diferencia entre el gasto pico calibrado y el gasto pico
simulado, los modelos pierden calidad. Igualmente este autor resume las ventajas y desventajas bajo
el empleo de modelos de tránsito de avenidas de una o dos dimensiones (Ver Tabla 1-1)
Paredes, Pedrozo y Mejía (2014) recomiendan para los modelos 2D utilizar una densidad de datos
elevada en los primeros 5 m respecto a la cota más baja de la llanura de inundación, de forma que el
cauce establezca las direcciones del flujo adecuadas. Para los modelos 1D se recomienda reducir el
intervalo entre secciones transversales, y tener especial cuidado al modelar las condiciones de frontera
de los modelos de las propiedades físicas (modelos de elevación, modelos de uso de suelo, modelos
de precipitaciones, etc.). También se hace notar que en casos donde los flujos bi-direccionales son
limitados, los resultados de las simulaciones usando modelos 1D y 2D son similares.
4
Tabla 1-1 Ventajas y desventajas del modelado 1D y 2D (Adaptada de Paredes et al, 2014)
Modelos Hidrodinámicos
Resuelve las ecuaciones de Saint-Venant en un sentido (1D) por
medio de una serie de secciones transversales del río y la llanura de
inundación al flujo del río.
Ventajas
Capacidad del sistema de No necesita altos requerimientos
cómputo de cómputo
Simulación de niveles y gastos Permiten una evaluación rápida
del río de la distribución del nivel de agua
y caudales del río
Desventajas
Una dimensión (1D) Simulación de Son menos precisos para modelar
desbordamiento del río el flujo de desbordamiento en ríos
Simulación de planicie o Son menos exactos para modelar
humedal con baja pendiente planicies de inundación
Simulación de estuarios y Son inciertos para modelar
flujos costeros estuarios y flujos costeros
Velocidad del flujo principal No se obtiene un dato aceptable
en el espacio de la velocidad, porque sólo se
tiene un sentido en la dirección del
flujo (X) (Ver Figura 1-1)
Resuelven las ecuaciones de aguas someras en 2D, en algunos casos
con modelos de cierre de turbulencia. Es posible discretizar la llanura
de inundación por medio de mallas regulares o adaptables
(rectángulos y triángulos, respectivamente)
Ventajas
Simulación de Modela de manera adecuada un
desbordamiento del río desbordamiento de río
Simulación de planicie o Modela de manera adecuada un
humedal con baja pendiente flujo sobre una planicie o humedal
con baja pendiente en el terreno
Simulación de estuarios y Se aplican con éxito en la
flujos costeros modelación de estuarios y flujos
Dos dimensiones (2D) costeros
Velocidad del flujo principal La velocidad en un punto se
en el espacio considera aceptable, ya que es
variable en cada dirección (X,Y)
por ser un vector
Desventajas
Capacidad del sistema de Su aplicación está limitada por los
cómputo altos requerimientos de datos,
hardware y software
Simulación de niveles y El tiempo que realiza una
gastos del río simulación de flujos o niveles de
agua en un río es considerable
Una comparación de ambos modelos hecha por Cook (2008), identifica las siguientes diferencias
significativas de resultados entre modelados efectuados bajo el empleo de modelos 1D y 2D.
5
b) Los modelos 1D son más susceptibles al cambio de precisión en los datos. Los modelos 2D
son menos susceptibles al cambio de precisión en la configuración de la malla que a los
cambios en la precisión del DEM.
c) Los modelos 1D interpolan los límites de la planicie de inundación entre cada par de
secciones transversales de terreno que tiene disponible, lo que regularmente deriva en áreas
de inundación discontinuas.
d) Los modelos 1D utilizan el coeficiente de fricción n de Manning para la calibración del
sistema de ecuaciones. Los modelos 2D utilizan el coeficiente de fricción n de Manning y el
coeficiente de viscosidad de remolino para calibrar las ecuaciones del sistema.
e) Los modelos 1D no tienen vectores afectados por la topografía del terreno, debido a que el
análisis se hace únicamente en la dirección del eje del río, de forma que el modelo sólo calcula
la magnitud de los vectores, pero la dirección y el sentido son definidos por la dirección del
eje del río.
Elevación de
terreno
Figura 1-1 Esquema de vectores de flujo resultantes bajo el empleo de un modelo 1D donde existe un
obstáculo.
De acuerdo con van Dijk et al (2014) los modelos 1D son suficientemente precisos para generar
simulaciones de flujo urbano detallado, esto debido a que las calles son estructuras artificiales que
tienen más en común con canales confinados que con el lecho natural de un río. De forma que modelar
las calles como una red de flujo 1D resulta adecuado y computacionalmente más eficiente. Aunque
en donde existen transiciones fuertes de pendiente, es mejor utilizar un modelo combinado 1D-2D.
Como explica Llamas (1993), los modelos tienen que ser calibrados con criterios deductivos o con
observaciones anteriores de solicitaciones y respuestas. Paredes, Pedrozo y Mejía (2014) proponen
calibrar los modelos hidrodinámicos de simulación de flujos utilizando el coeficiente de rugosidad de
6
Manning (n). Éstos autores señalan que para valores más elevados del coeficiente de rugosidad, se
genera una mayor área afectada.
Bajo este principio, se comparan los resultados de la simulación hidrodinámica, que ha sido realizada
utilizando los datos hidrométricos de un periodo en particular, con la llanura de inundación
correspondiente a ese mismo periodo. Entonces, es posible afinar el coeficiente de rugosidad en cada
elemento de análisis (secciones transversales para las simulaciones con modelos 1D, elementos de la
malla para las simulaciones con modelos 2D) como parte del proceso de calibración (Karim et al,
2011).
Sin embargo, se ha demostrado que la mayor parte de la información requerida para la delimitación
de inundaciones puede ser obtenida de bandas satelitales que detecten la reflexión infraroja de rango
medio (TM: 1.55–1.75 mm) del terreno natural, por lo que, con el fin de reducir tiempos y costos, las
imágenes de la banda térmica 5 de LANDSAT ha sido elegida como criterio único para determinar
la presencia de agua sobre el terreno (Overton, 2005).
7
Una vez calibrado el modelo, necesita ser validado. Autores como Paredes et al (2014) y van Dijk et
al (2014) utilizan una simulación enteramente realizada con modelos 2D para validar los modelos de
simulación combinados.
Incluso en el levantamiento LiDAR (Laser Detection and Ranging), los cuerpos de agua regresan
valores nulos. Es por esta razón que la información de los fondos de los cuerpos de agua deben ser
completados. Algunos autores, como Coleman (2006) y Wang (2012) simularon los datos faltantes
usando una combinación de sonares e interpolaciones ponderadas (tipo “vecino natural”).
Ambos métodos generan buenos resultados al reconstruir los datos batimétricos, pero requieren de un
estudio por sonar independiente de los levantamientos LiDAR. Corum (2015) propone un método
para crear una batimetría sintética, cuya metodología se aborda en 2.7 Método de obtención de
batimetría sintética propuesto por Corum.
1.3 Justificación
Se han descrito las ventajas y desventajas de cada tipo de modelo. De ahí se desprende la importancia
de desarrollar estrategias y diseñar criterios que puedan ser fácilmente empleados durante el análisis
del tránsito de avenidas y de inundaciones, para definir qué modelo o combinación de modelos es
apropiada para optimizar los tiempos de cómputo a la vez que se logra la calidad del resultado.
8
Por combinación de modelos, debe entenderse que, para un mismo tramo de río, existirán subtramos
que por sus características resulte más conveniente simular usando modelos 1D o 2D. Esto resulta en
que, en un mismo tramo, coexistan simulaciones en serie, donde los datos de salida de una simulación,
son los datos de entrada de otra. Tal intercambio de información se da a través de zonas denominadas
fronteras, que son las líneas de unión entre subtramos con simulaciones de distinto modelo.
Si se busca además una automatización de este proceso, habrá que tomarse en cuenta las fronteras de
éstos criterios para dividir la zona de estudio en segmentos independientes que tendrán que analizarse
con un modelo en particular, y los cuales deben integrarse al final en un sólo análisis hidrodinámico.
Los resultados del análisis, particularmente la cartografía de la inundación, pueden ser exportados
para utilizarse en otros análisis, tales como la estimación de costos.
Debe destacarse que el presente protocolo de tesis de Maestría en Ciencias del Agua forma parte del
proyecto de investigación “Modelación de eventos hidrológicos extremos a partir de
precipitación estimada por percepción remota” financiado por CONACyT, con clave UAEM:
3952/2015E en el marco de la convocatoria “Proyectos de Desarrollo Científico para Atender
Problemas Nacionales, 2014” (Bâ, K.M., Díaz-Delgado, C., Franco Plata, R., Gómez, M., Fonseca,
C.R, Mastachi Loza, C.A., Torres Maya, A., Vilchis-Francés, A.Y.). Otro proyecto de investigación
que complementa este proyecto es: Diseño de herramienta hidrológica como apoyo a los sistemas de
alerta temprana ante inundaciones (López García, A., Díaz-Delgado, C., Bâ, K.M. y Gómez, M.,).
Igualmente, este trabajo forma parte del proyecto UAEM 4192/2016E “Hidrología en Mauritania:
Modelación de caudales a partir de precipitación estimada por imágenes de satélite”.
1.4 Hipótesis
Obtener secciones transversales sintéticas en las áreas faltantes del modelo de elevación que
no afecten los resultados de la simulación hidrodinámica.
Simular el tránsito de avenidas con modelos 1D y 2D (HEC-RAS) para comparar los
resultados en función de la calidad de los datos topo-batimétricos disponibles y características
morfológicas del terreno inundable.
Deducir estrategias de simulación hidrodinámica del tránsito de avenidas para definir en qué
casos es más conveniente utilizar cada tipo de simulación.
Deducir criterios de aplicación para el uso de cada tipo de modelo (1D y 2D) para simular el
tránsito de avenida con la calidad adecuada.
10
Marco teórico
El modelo HEC-RAS (River Analysis System) apoyádo en el software del mismo nombre, en su
versión 5.0, es una herramienta de software de uso libre, que permite realizar simulaciones hidráulicas
de flujo constante usando modelos 1D y de flujo no constante de agua usando modelos 1D, 2D y la
combinación 1D-2D (Lluén, W., 2015, Brunner, G., 2016). y su desarrollo es por parte del
Hidrologic Engineering Center del U.S. Army Corps of Engineers,
El análisis de flujo estable (1D) en una dirección tiene la finalidad de determinar la elevación del
agua en flujos estables gradualmente variados. El principio básico del análisis de flujo estable dentro
de este software es la resolución de la ecuación de la energía en una dirección. Las pérdidas son
evaluadas de acuerdo con el coeficiente de fricción de Manning y coeficientes de expansión y
contracción. Los saltos hidráulicos, las confluencias de ríos y el análisis de la hidráulica sobre
estructuras hacen uso de la ecuación del momento. El software tiene capacidad de evaluar un canal o
una red de canales naturales o artificiales (Brunner, G., 2016).
El procedimiento de análisis del flujo se da por el método de pasos estandar (ver Figura 2-1)
El análisis de flujo no estable en 1D se hace mediante el modelo UNET desarrollado por el Dr.
Robert L. Barkau (Barkau, 1992), basando la simulación de flujo entre secciones en la resolución de
las ecuaciones de continuidad de masa (Ecuación 2-8) y la de continuidad de momento (Ecuación
2-9), siguiendo el algoritmo descrito en la Figura 2-2
11
Figura 2-1 Método del paso estándar para la solución de la ecuación de la energía para flujos estables en
1D
12
Inicio
Se asume el tirante
inicial y se calculan los
elementos de flujo (Ecs.
8 y 9)
Los elementos de
cálculo (secciones
transversales y/o celdas Se realiza una nueva
de malla 2D) se estimación del tirante
comparan con la No
tolerancia numerica
Si
Se muestra una ¿Algún valor de error
Se procede con el
advertencia con el es menor a la
siguiente paso de Si
error que supera la tolerancia de
tiempo
tolerancia inestabilidad?
No
Fin
Figura 2-2 Algoritmo para solucionar las ecuaciones de flujo no estable (1D y combinación 1D/2D)
13
El modelo de análisis de flujo no estable en 2D fue desarrollado dentro del HEC y basa la simulación
en la solución de las ecuaciones de conservacion de masa y momento (Ecs. 10 a 14) en dos
dimensiones.
Los elementos de cálculo del flujo 2D son las celdas de la malla trazada en la zona de estudio. Una
contribución importante del cálculo efectuado por el modelo 2D de HEC-RAS es el hecho de que
además de utilizar éstas celdas, si se cuenta con una batimetría com mayor resolución, utilizará
algunos de los datos de esta batimetría que estén dentro de la celda a modo de una “sub-malla” (Ver
Figura 2-3) para mejorar los resultados de la simulación
Malla modelo 2D
Sub-malla a partir de batimetría
Figura 2-3 Ejemplo de datos batimétricos dentro de la malla de cálculo para la simulación 2D
Si bien esta sub-malla no es utilizada directamente como elemento de cálculo, algunos datos como la
elevación y el coeficiente de fricción de manning se conservan para obtener un coeficiente de fricción
compuesto y un área de canal mejores que con el simple uso de la malla 2D.
El cómputo de la simulación de la malla, se realiza creando “canales” entre los centroides de celdas
con aristas contiguas de la malla y entre cada uno de éstos resolviendo las ecuaciones de continuidad
de masa y momento (Ver Figura 2-4)
14
Malla de cálculo 2D
Canal entre centroides
Centroides de celdas
Hi Elevación en punto i
Hj Elevación en punto j
n Dirección del flujo inicial
n’ Dirección del flujo final
T’ Cambio de dirección del flujo
El algoritmo de simulación para flujo no estable en dos direcciones (Ver Figura 2-5)
Inicio
Se computan las
velocidades de
acuerdo a las Ecs. 9 y
Se pre-computan la 10
geometría, la
ortogonalidad de las
celdas y la sub-malla por
batimetría
Mostrar mensaje
de error mayor a la
tolerancia
A partir de las elevaciones
propuestas se calculan las Si
propiedades hidráulicas del No
canal (Áreas, Radios, etc.)
Si
Se procede al
siguiente paso
de tiempo
Se determina la continuidad
de momento y volumen en
todas las celdas (para
conocer la dirección en la
que se realizará la Fin
simulación)
15
2.2 Principios matemáticos de las simulaciones hidrodinámica de flujo
Los modelos hidrodinámicos unidimensionales funcionan resolviendo las ecuaciones de <0o de onda
dinámica, que corresponden a las ecuaciones de momento y de continuidad en una dirección (Timbe,
2011).
𝑎 𝑉 𝑎 1
𝑍 +𝑌 + =𝑍 +𝑌 + +ℎ Ecuación 2-2
2𝑔 2𝑔
𝑎 𝑉 𝑎 𝑉
ℎ = 𝐿𝑆 + 𝐶 − Ecuación 2-3
2𝑔 2𝑔
𝐿 𝑄 +𝐿 𝑄 +𝐿 𝑄 Ecuación 2-4
𝐿=
𝑄 +𝑄 +𝑄
𝑆 +𝑆
𝑆 = Ecuación 2-5
2
Donde:
Para tener una mejor aproximación de los datos de cada sección, ésta es dividida en varias secciones,
generalmente en cada punto del perímetro mojado (Ver Figura 2-6). Para obtener el coeficiente de
16
ponderación, determinamos la velocidad, el gasto y la carga por velocidad en cada una de las
subdivisiones, y se emplea la ecuación 6.
Fondo de la S.T.
Espejo de agua
Figura 2-6 División de la sección trasversal para determinar la velocidad ponderada (Ver Ecuación 2-6)
𝑄 𝑉 𝑄 𝑉 𝑄 𝑉
+ +⋯
2𝑔 2𝑔 2𝑔 Ecuación 2-6
𝑎=
(𝑄 + 𝑄 + ⋯ 𝑄 ) ∗ 𝑉
La Ecuación 2-7 describe la conservación de momento entre secciones transversales para el flujo
estable.
𝑄 𝛽 𝐴 +𝐴 𝐴 +𝐴 𝑄 𝛽
+𝐴 𝑌 + 𝐿𝑆 − 𝐿𝑆̅ = +𝐴 𝑌 Ecuación 2-7
𝑔𝐴 2 2 𝑔𝐴
Donde:
17
La Ecuación 2-8 describe la continuidad (conservación de masa en 1D).
ΔA ΔA ΔS
Δ𝑄 + Δx + + Δx − Q = 0 Ecuación 2-8
Δt Δt Δt
Donde:
Δ 𝑄 Δx + Q Δx
+ Δ(βVQ) + gAΔz + gAS Δx = 0 Ecuación 2-9
Δt
𝑔𝐴 𝑆 Δx + 𝑔𝐴 𝑆 Δx
Δx = Ecuación 2-10
𝑔𝐴̅𝑆
Donde:
Para flujo no estable, la Ecuación 2-11 describe la conservación de momento en 2D para velocidad
en sentido u, mientras que la Ecuación 2-12 describe la ecuación de momento en 2D para la
velocidad en sentido v.
18
𝛿𝑣 𝑢𝛿𝑣 𝑣𝛿𝑣 𝑔𝛿𝐻 𝑢𝛿 𝑣 𝛿 𝑣
+ + =− +𝑣 + − 𝑐 𝑣 + 𝑓𝑢 Ecuación 2-12
𝛿𝑡 𝛿𝑥 𝑑𝑦 𝛿𝑥 𝛿𝑥 𝛿𝑦
𝑛𝑔 /
𝑢∗ = |𝑉| Ecuación 2-14
𝑅 /
𝑔|𝑉|
𝑐 = Ecuación 2-15
𝐶 𝑅
Donde
u,v: Velocidades del flujo en las direcciones cartesianas del plano de la malla
vt: Coeficiente de viscosidad de torbellinos
f: Parámetro de Coriolis
cf: Fricción del fondo del canal
t: Es el tiempo de simulación
H: Elevación de la superficie de agua
R: Es el radio hidráulico
𝑢∗ : Velocidad de corte
g: Aceleración de la gravedad
|𝑉| : Magnitud del vector de velocidad
D: Constante empírica
h: Tirante hidráulico
C: Coeficiente de Chézy
19
2.3 Modelo CEQUEau
El modelo CEQUEAU, que es empleado para la simulación de escurrimiento dados los datos
hidrometeorológicos obtenidos en tiempo real, es un modelo hidrofisiográfico capaz de simular
caudales simultáneamente en tiempo y espacio (CHARBONNEAU, 1977).
Este modelo basa su funcionamiento en la discretización de la cuenca de estudio por medio de una
cuadrícula, de la cual se obtendrá una síntesis del modelo de escurrimiento por medio de una
subdivisión de esta cuadrícula. Luego se asignan estaciones meteorológicas y estaciones
hidrométricas a los elementos de la cuadrícula correspondiente.
El programa cuenta con dos funciones principales, la función de producción, que simula el flujo
vertical del agua en el suelo, esquematizado por varios reservorios intercomunicados entre ellos y la
segunda, llamada función de transferencia, enfocada en el escurrimiento usando factores de drenaje
y factores fisiográficos de la cuenca (Morín. y Paquet. , 2007).
El modelo discretiza las caracterísitcas físicas de la cuenca en cuadros que a su vez pueden ser
divididos hasta en 4 sub-cuadros (denominados A,B,C y D) y utilizando como datos de entrada la
elevación del terreno (m) mediante un archivo que contenta la información fisiográfica del terreno,
las precipitaciones (mm) de las estaciones climatológicas pertinentes cercanas a la cuenca de estudio,
las temperaturas (°C) tomadas también de las estaciones climatológicas, y los caudales (m3/s)
obtenidos de las estaciones hidrométricas cercanas la desembocadura (Diaz., 2010).
20
Dentro de los cuadros en los que se ha discretizado la zona de estudio también es necesario indicar el
uso de suelo. Para CEQUEau esto se hace indicando para cada cuadro el porcentaje que contiene de
los cuatro usos de suelo que puede diferenciar el programa:
1. Cuerpos de Agua
2. Bosques
3. Ciénegas
4. Otros
Como Morín y Paquet (2007) señalan, la calibración del modelo puede realizarse por prueba y error,
o por optimización, utilizando el módulo que el programa incluye para tal fin. La calibración por
prueba y error consiste en hacer variar los parámetros y ver el efecto que tienen tales modificaciones
sobre el resultado de la simulación. Hay algunas prácticas que facilitan el proceso de calibración,
tales como:
Si bien no existen reglas para los demás parámetros, la experiencia y la experimentación permiten
realizar un buen ajuste.
Una vez que se han definido los datos fisiográficos y se ha calibrado el modelo CEQUEau, es posible
posicionar estaciones ficticias, para las cuales de acuerdo a la interpolación combinada de datos
fisiográficos y meteorológicos, pueden obtenerse estimaciones de gastos máximos diarios que
corresponden a la posición de éstas estaciones ficticias.
2.4 Matlab
Matlab es un lenguaje de scripting científico, con un entorno de desarrollo integrado (un conjunto de
editor de código fuente, un depurador, compilador y librerías para automatizar el desarrollo), que
trabaja en lenguaje M (lenguaje propietario de Mathworks, empresa desarrolladora del software),
teniendo Matlab la característica de contar con una amplia variedad de herramientas incluidas
orientadas a matemáticas, pudiendo expandir su funcionalidad con paquetes de herramientas
21
adicionales (toolboxes), que por la naturaleza de MATLAB, vienen en forma de código fuente
(Borrel., 2008).
Principalmente, el trabajo de datos en Matlab es por medio de vectores (incluso las matrices, no son
sino vectores con formateo suficiente para poder ser utilizados como matrices, aunque como
reminiscencia de su naturaleza vectorial, es posible acceder a los datos por medio de las coordenadas
matriciales o bien, con la posición escalar de los elementos dentro del vector).
Como se ha mencionado, más herramientas pueden anexarse al acervo de funciones de Matlab por
medio de las toolboxes, aunque la versión base del programa ya cuenta con una gran cantidad de
éstas. De éstas, las principales que se utilizarán son las siguientes:
Geotiffread: Esta herramienta permite importar la matriz de datos de una imagen Geotiff,
así como las etiquetas correspondientes a la información de referencia
geoespacial de la imagen (Límites X e Y del DEM, Resoluciones en X e Y,
etc.)
Geotiffinfo: Esta herramienta importa las etiquetas correspondientes a las propiedades de la
imagen Geottif (Resolución, Cantidad de Bits, tipo de compresión, nombres
de las unidades, nombres de las referencias, etc.)
Geotiffwrite: Usando como datos de entrada la matriz de datos de elevaciones (Mat), el
arreglo de propiedades de georreferencia (Ref), y la información de las
propiedades de la imagen (Inf) puede escribir un archivo en formato Geottif.
Textread: Crea una matriz de datos a partir de un archivo de texto, donde cada fila se
delimita por los espacios en el archivo de texto y cada columna se delimita por
un espaciador especificado por el usuario (regularmente, una coma, un
espacio o un tabulador)
Fprintf: Escribe datos especificados por el usuario en una línea de un archivo de texto
22
El conjunto de herramientas externas que es utilizado dentro de este proyecto es Topotoolbox (ver
Herramientas Topotoolbox, pag. 23).
Las matrices bidimensionales de Matlab se numeran de izquierda a derecha (Eje X) y de arriba abajo
(eje Y) (ver Figura 2-7 Importación de Geotiff a MATLAB para el uso en Topotoolbox Es importante
tener en cuenta esta peculariedad de que los valores en Y incrementan hacia abajo, puesto que en los
archivos GeoTiff este valor se incrementa hacia arriba.
Es por esto, que para una coordenada en formato matricial (Y,X) cualquiera, su correspondiente
coordenada (XUTM,YUTM) puede obtenerse mediante la Ecuación 2-16, o bien, se puede estimar la
coordenada matricial mas cercana (Y,X) a partir de coordenadas (X UTM,YUTM) con la Ecuación 2-17.
Ecuación 2-17
𝑌 − 𝑅𝑒𝑓𝑌 𝑋 − 𝑅𝑒𝑓𝑋
𝑓(𝑋 ,𝑌 ) = 𝑟𝑒𝑑𝑜𝑛𝑑𝑒𝑎𝑟 , 0 , 𝑟𝑒𝑑𝑜𝑛𝑑𝑒𝑎𝑟 ,0
𝑅𝑒𝑠𝑌 𝑅𝑒𝑠𝑋
Donde:
Las herramientas que componen topotoolbox tienen capacidad de análisis de dirección de flujo,
sedimentos, químicos y nutrientes.
23
En el presente trabajo, las herramientas utilizadas fueron las relacionadas con el análisis de dirección
de flujo, acumulación de unidades de flujo y limpieza de valores anómalos que podrían generar ciclos
(PitRemoval).
La importación de información geoespacial se hace por medio de archivos tipo Geotiff. Por medio de
Matlab, la información de la matriz numérica de elevaciones es guardada dentro de una matriz,
mientras que los datos geográficos son almacenados dentro de diferentes arreglos, para utilizarse
posteriormente (Ver Figura 2-7)
La dirección del flujo tradicionalmente es calculada desde un píxel a otro, calculando el gradiente de
elevación (ver Ecuación 2-18) en los pixeles adyacentes y dirigiendo el flujo hacia el píxel adyacente
que presente el mayor gradiente.
𝑍𝑖 − 𝑍𝑗
𝑆 = Ecuación 2-18
𝑑
Donde:
Sij: Gradiente de elevación entre los pixeles i y j
Zi, Zj: Elevación de los pixeles i y j, respectivamente
dij: Distancia entre los pixeles i y j
Inicialmente, cada píxel tiene una unidad de flujo. Una vez que se ha calculado la dirección del flujo
de un pixel a otro, la cantidad de unidades de flujo que se encuentren acumulados en el pixel inicial,
se suman a la cantidad de unidades de flujo acumuladas en el pixel final (Ver Figura 2-8).
24
Sin embargo, según señalan Schwanghart y Kuhn (2010), este método termina ignorando los flujos
cruzados y propicia sobre concentraciones de flujo en las zonas bajas de la topografía. Es por esto por
lo que en TopoToolBox se implementó un algoritmo de escurrimiento múltiple, donde para cada píxel
se crea una matriz M, denominada matriz de transferencia, que señala la transferencia de unidades de
flujo hacia hasta ocho pixeles adyacentes, proporcionalmente a la intensidad de la pendiente hacia
tales pixeles.
Figura 2-8 Acumulación de unidades de flujo dentro del DEM (Imagen tomada de Schwanghart y Kuhn, 2010)
Para mantener la continuidad de la transferencia de flujo a través de zonas planas (secciones del
modelo de elevación en el que el gradiente de elevación es igual a cero), topotoolbox tiene dos
algoritmos diferentes, Routeflats y Crossflats.
Ambos algoritmos comienzan las iteraciones identificando valores nulos (NaN) y agujeros (pixeles
que no se encuentran en las orillas del DEM, que alrededor no tienen ningún píxel con una elevación
igual o menor) para ser evitados durante el cálculo.
El algoritmo Routeflats consiste en delimitar el borde de la zona plana, siendo la “salida” del flujo la
elevación más baja del borde, y luego determina las rutas dentro de la zona plana tomando como
inicio los pixeles que transfieren su flujo a los pixeles que conforman el borde, y después se encuentra
la ruta a lo largo de gradientes interpolados del borde, mejorando las rutas subsecuentes con cada ruta
calculada.
25
El algoritmo Crossflats de igual manera delimita el borde de la zona plana, determinando la “salida”
del flujo. Después, simplemente une cada entrada del flujo con la salida por medio de una línea recta.
Contemplan un amplio espectro de funciones sobre el manejo del agua. En particular, el presente
trabajo hace uso de las herramientas que posee acerca de las propiedades hidrofisiográficas del
terreno, particularmente, las propiedades de área de la cuenca, tiempo de concentración por el método
de Kirpich y longitud del cauce principal.
El levantamiento por medio de LiDAR no puede devolver datos, o devuelve datos insuficientes en
cuerpos profundos de agua debido a la baja reflectancia (Schenk y Csatho, 2002). Además, los datos
devueltos por el dispositivo LiDAR son puntos ubicados al azar sobre el terreno, por lo que no pueden
ubicar correctamente las orillas y bordes del terreno en elementos como ríos, crestas o construcciones.
Es por eso que, en cuerpos de agua, como ríos, si solo se dispone de levantamientos por medio de
LiDAR, es probable que éstas zonas no cuenten con cuerpos de agua.
Es por esto, que Corum (2015) propone un método para determinar una batimetría sintética si se
conoce la elevación de la superficie de agua y el gasto correspondiente a tal elevación.
26
Figura 2-9 Algoritmo de Corum para obtener un lecho sintético
Aún queda el problema de la distribución del gasto conocido a lo largo del tramo, principalmente
cuando una parte importante de la aportación del gasto proviene de canales adyacentes tributarios. En
tales casos, puede elegirse entre aislar el flujo completo dentro del canal principal o elaborar un
esquema con un tramo de río secundario conectado al principal (Corum, 2015).
En el presente trabajo, se eligió realizar la distribución del gasto conocido a lo largo del río, mediante
la simulación de gasto en estaciones ficticias utilizando el modelo CEQUEau (ver Modelo CEQUeau,
pág. 20) para obtener una estimación de los porcentajes del gasto conocido correspondientes a cada
tramo del río (la discretización de éstos tramos está limitada por las restricciones de discretización
27
fisiográfica del software CEQUEau (ver Modelo CEQUeau, pag. 20)), lo cual incluye las
aportaciones de los ríos tributarios.
Tanto para la simulación de caudales usando CEQUEau (Pág. 20) como para el tránsito de avenidas
utilizando HEC-RAS (Pag. 11,) es necesario realizar análisis estadísticos de los datos
hidrometeorológicos e hidrométricos.
Primeramente, para comprobar la calidad de las series de datos obtenidos se realizarán pruebas de
aleatoriedad, de homogeneidad y de independencia.
Una muestra aleatoria, desde el punto de vista estadístico, es una en la que cada una de las muestras
tiene la misma probabilidad de ser escogida. En hidrología, la aleatoriedad asegura que los datos son
debidos a causas naturales (Cano, 2005).
Cada observación, se evalúa respecto a cada observación que le precede. Es decir, para cada
elemento Xj, se evalúan los elementos X1, X2,… Xj-1.
NS=0
Ni=0
28
Si el elemento Xj es el máximo dentro del grupo evaluado (es decir Xj = max(X1, X2,… Xj)) entonces
se agrega una unidad al contador NS ( NS = NS +1), si por el contrario es el mínimo (es decir Xj =
min(X1, X2,… Xj)) se agrega un valor al contador Ni (Ni= Ni+1). Si no nos encontramos en ninguno
de estos dos casos, no se toma ninguna acción y se procede a evaluar el siguiente elemento.
Se repite el proceso hasta terminar con la serie, y después se obtiene el estadístico ε mediante la
Ecuación 2-19.
𝑁 −𝑁
𝜀= Ecuación 2-19
2 ∗ ln(𝑛) − 0.845
Se supone que la muestra tiene una distribución Normal. Si se considera un valor de significancia α,
entonces la muestra se considera aleatoria si se cumple
−𝑍 ≤ 𝜀 ≤ 𝑍
Las pruebas de homogeneidad determinan indican si los datos probados pertenecen a una misma
muestra. Particularmente, indica si es que los datos probados se desenvuelven alrededor de una misma
media. Esto tiene significancia hidrológica principalmente debido a las modificaciones en el ciclo
hidrológico que sufre una zona debido a la deforestación, la urbanización y el cambio climático. Es
decir, si a lo largo del periodo en el que fueron tomados los datos, ha cambiado la tendencia de los
datos (Por ejemplo, debido a la urbanización y deforestación de una cuenca, con unas décadas el
escurrimiento esperado para un periodo de retorno puede aumentar drásticamente). De ser así, se
podría llegar al caso en que datos más antiguos tengan que descartarse del análisis estadístico, dado
que ya no son representativos para el ciclo hidrológico de la cuenca.
29
si Xj<𝑥 , se señala con un 0.
Ni0=0
Ni1=0
N0= Número de ceros señalados en la muestra
N1= Número de unos señalados en la muestra
Entonces, se hace un recorrido sobre el vector de 0’s y 1’s obtenidos de la muestra, en orden
cronológico. Por cada cambio que haya de 0 a 1, se añade una unidad a Ni1 (Ni1=Ni1+1), por cada
cambio que haya de 1 a 0, se añade una unidad a Ni0 (Ni0=Ni0+1). Se repite el proceso hasta terminar
con la serie.
(𝑁 + 𝑁 ) − 𝑥
𝜀= Ecuación 2-20
𝑆
2𝑁 𝑁 Ecuación 2-21
𝑥 = +1
𝑛
Donde:
30
Se considera que la prueba sigue una distribución Normal, por lo que para un nivel de significancia
α, la muestra se considera homogénea si se cumple que
−𝑍 ≤ 𝜀 ≤ 𝑍
La prueba empleada dentro de la investigación fue la de Wald y Wolfowitz (recopilada por Cano,
2005).
Para una muestra de tamaño n esta prueba considera la media y la varianza con los siguientes valores
𝜓 −𝜓 Ecuación 2-23
𝑥 =
𝑛−1
Ecuación 2-25
𝜓 −𝜓 𝜓 − 4𝜓 𝜓 + 4𝜓 𝜓 + 𝜓 − 2𝜓
𝑆 = + −𝑥
𝑛−1 (𝑛 − 1)(𝑛 − 2)
𝑤= 𝑥𝑥 +𝑥 𝑥 Ecuación 2-24
Donde
𝜓 : 𝑛𝑚′
31
𝑛: Tamaño de la muestra
𝑚′ : j-ésimo momento respecto al origen
Entonces, el estadístico de prueba para la prueba de Wald y Wolfowitz está dado por la Ecuación
2-26.
𝑤−𝑥
𝜀= / Ecuación 2-26
𝑆
Se supone que la muestra tiene una distribución Normal. Si se considera un valor de significancia α,
entonces la muestra se considera independiente si se cumple que
−𝑍 ≤ 𝜀 ≤ 𝑍
Una vez que la calidad de los datos ha sido evaluada, se procede a realizar un análisis estadístico de
los datos.
Para el análisis, se realizan diferentes ajustes probabilísticos, que después son comparados con los
datos observados para determinar la bondad de ajuste que presenta cada función y elegir la que se
apegue de mejor manera a los datos observados.
García (2010) ofrece una recopilación de las principales distribuciones estadísticas empleadas en
modelos hidrológicos, mismas que son listadas abajo.
1 ( )
𝑓(𝑥) = 𝑒 Ecuación 2-27
√2𝜋𝜎
32
Donde
(𝑎 𝑥 𝑒 )
𝑓(𝑥) = Ecuación 2-29
Γ(𝑏)
𝑥̅
𝑎= Ecuación 2-28
𝑠𝟐
𝑥
b= Ecuación 2-30
𝑆
La distribución de densidad de la función Pearson III está dada por la Ecuación 2-32.
( )
(|𝑐|(𝑥 − 𝑎𝑝) 𝑒 )
𝑓(𝑥) = Ecuación 2-32
Γ(𝑏𝑝)
2𝑆
𝑎𝑝 = 𝑥̅ − Ecuación 2-33
𝐶
4
𝑏𝑝 = Ecuación 2-34
𝐶
2
𝑐𝑝 = Ecuación 2-35
𝑆𝐶
33
/
𝐶 = 𝑚 /𝑚 Ecuación 2-36
𝐶 𝑛(𝑛 − 1)
𝐶 = Ecuación 2-37
𝑛−2
Donde:
La función de densidad de la distribución Gumbel Tipo I está dada por la ecuación Ecuación 2-38
1 Ecuación 2-38
𝑓(𝑥) = 𝑒 ( )
𝑏
√6𝑆
𝑏= Ecuación 2-40
𝜋
Donde:
La función de densidad de la distribución log-Normal está dada por la ecuación Ecuación 2-42
34
1
𝑓(𝑥) = 𝑒 Ecuación 2-42
√2𝜋𝜎
Donde
y=ln(x)
𝜎 : Desviación estándar de la transformada “y” de la muestra
La distribución de densidad de la función log-Pearson III está dada por la Ecuación 2-43.
( )
(|𝑐𝑙𝑝|(𝑦 − 𝑎𝑙𝑝) 𝑒 )
𝑓(𝑥) = Ecuación 2-43
Γ(𝑏𝑙𝑝)
2𝑆
𝑎𝑙𝑝 = 𝑦 − Ecuación 2-44
𝐶
4
𝑏𝑙𝑝 = Ecuación 2-45
𝐶
2
𝑐𝑙𝑝 = Ecuación 2-46
𝑆𝐶
/
𝐶 = 𝑚 /𝑚 Ecuación 2-47
𝐶 𝑛(𝑛 − 1)
𝐶 = Ecuación 2-48
𝑛−2
35
Donde:
y=log(x)
Cs: Coeficiente de asimetría de la transformación “y” de la muestra
Csc: Coeficiente de asimetría corregido de la transformación “y” de la muestra
alp, blp y clp: Parámetros de la función de densidad de la distribución log-Pearson III
n: Tamaño de la muestra
Pero mientras que los hidrogramas unitarios aplican únicamente para la cuenca y el punto donde
fueron medidos los caudales, los procedimientos de desarrollo de hidrogramas unitarios sintéticos
permiten desarrollar hidrogramas unitarios para otros puntos en la corriente de la misma cuenca
(Chow, 2000).
El hidrograma unitario adimensional SCS, es un H.U. sintético que puede ser estimado a partir del
tiempo de ocurrencia del pico y el caudal pico (Chow, 2000), cuya relación puede apreciarse en la
Ecuación 2-49
𝐶𝐴𝑐
𝑞 = Ecuación 2-49
𝑇
𝑡
𝑇 = +𝑡 Ecuación 2-50
2
36
𝑡𝑟 = 0.133𝑇 Ecuación 2-52
Donde:
Una vez calculados los parámetros de forma, éstos son multiplicados por las coordenadas del
hidrograma adimensional del SCS mostrado en la Figura 2-10 Diagrama adimensional de la SCS,
adaptado por Chow, (2000) del Soil Conservation Service Engineering Handbook, 1972. tabulado en
la Tabla 2-2.
0.5
0.4
0.3
0.2
0.1
0
0 1 2 3 4 5
t/Tp
Figura 2-10 Diagrama adimensional de la SCS, adaptado por Chow, (2000) del Soil Conservation
Service Engineering Handbook, 1972.
37
Tabla 2-2 Valores tabulados del diagrama adimensional de la SCS
t/tp q/qp
0 0
0.2 0.075
0.4 0.28
0.6 0.6
0.8 0.89
1 1
1.2 0.92
1.4 0.75
1.6 0.56
1.8 0.42
2 0.32
2.2 0.24
2.4 0.18
2.6 0.13
2.8 0.098
3 0.075
3.5 0.036
4 0.018
4.5 0.009
5 0.004
Una vez que se ha obtenido el hidrograma unitario, para adaptarlo cualquier otro volumen, se
multiplica el caudal pico unitario, qp (Ecuación 2-49) por la lámina correspondiente en centímetros.
38
Materiales y métodos
La presente investigación pretende determinar las situaciones en las que es apropiado utilizar
secciones de modelo 1D y 2D, de forma que se obtenga un esquema combinado de simulación tan
eficiente como se pueda (Ver Figura 3-1).
1D
1D
2D
1D
Para este fin, es necesario llevar a cabo un modelo de investigación que permita alcanzar los objetivos
planteados (Ver Pág. 10). Es necesario entonces considerar dentro de tal metodología la investigación,
los procesos de simulación y el desarrollo de herramientas informáticas necesarios para cumplir con
éstos.
39
Recopilación de datos
•Datos meteorológicos
•Datos hidrométricos
•Datos batimétricos
•Datos de reflectancia (LANDSAT 5)
Tratamiento de datos
•Tratamientos estadísticos
•Tratamientos geoespaciales
•Corrección de valores anómalos
•Desarrollo de herramientas para automatizar el
tratamiento de datos
40
3.1 Recopilación de datos
Los datos de elevación de alta resolución (LiDAR a cada 5 m) correspondientes a la zona de estudio,
pueden ser obtenidos en la página de INEGI, en la sección de topografía.
Por su parte, las imágenes raster de uso de suelo, que serán utilizadas para indicar el coeficiente de
fricción de Manning dentro de las simulaciones de flujo, pueden obtenerse de la página de CONABIO
(http://www.conabio.gob.mx/informacion/gis/).
41
La base de datos y el archivo vectorial con las ubicaciones de las estaciones hidrométricas fueron
descargadas desde la página de CONABIO (http://www.conabio.gob.mx/informacion/gis/).
Los datos hidrométricos para esta estación, pueden ser obtenidos de la página de BANDAS
(ftp://ftp.conagua.gob.mx/Bandas/Bases_Datos_Bandas).
Los datos meteorológicos para esta estación fueron descargados de la página de CLICOM
(http://clicom-mex.cicese.mx/, en la pestaña descarga de datos).
Para la calibración del modelo, como se mencionó en la sección de antecedentes, es necesario el uso
de las imágenes de la banda térmica 5 de LANDSAT. Éstas imágenes es posible conseguirlas de la
página EARTEXPLORER (www.earthexplorer.org).
Una vez obtenidos los datos necesarios para las actividades de la investigación (Ver Recopilación de
datos, pág. 41), éstos serán sometidos a un proceso que garantice su calidad, que los ordene en el
formato adecuado y que extraiga de éstos la información requerida para la investigación (Ver
Herramientas Topotoolbox, Pag. 23, Método de elaboración de batimetría sintética de Corum, pág.
26 y Métodos estadísticos de evaluación de datos, Pag. 28).
De acuerdo con lo que se aprecia en la Tabla 3-1 Datos requeridos por cada actividad de la
investigación la primera actividad consiste en la delimitación de la cuenca de estudio, para lo que es
necesaria la información batimétrica de la zona.
El tratamiento que recibirán los datos batimétricos obtenidos como se indica en la sección “3.1
Recopilación de datos”, consiste la remoción de datos que impidan el flujo hidráulico a través del
DEM, que se consigue usando la herramienta TopoToolbox (2.5 Herramientas TopoToolbox, pág.
23), la delimitación de la cuenca y la extracción de los datos fisiográficos correspondientes de esta
fueron conseguidos utilizando las herramientas Integrated Water Management.
42
Por otro lado, para generar el archivo DEM en un formato de datos fisiográficos que pueda utilizar
CEQUEau, fue necesario el uso de las herramientas de Idrisi-CEQUEau.
Si bien es posible utilizar los datos batimétricos con correcciones relativamente pequeñas para su uso
dentro de CEQUEau y la delimitación de la zona de estudio, para las simulaciones hidrodinámicas,
es necesario además corregir el DEM generando un lecho sintético como se menciona en Método de
obtención de batimetría sintética propuesto por Corum, Pág. 26.
Aunque para llevar a cabo este método, ha sido necesario primero pasar el DEM de datos batimétricos
por un proceso de corrección de datos nulos, de forma que se cumpla con el requisito establecido por
el método de Corum de establecer la superficie conocida del espejo de agua.
Dado que la superficie del espejo de agua es representada como valores nulos (Ver. Método de
obtención de batimetría sintética propuesto por Corum, pág. 26), es necesario corregir el DEM de
forma que los valores nulos representen el espejo de agua. En términos generales, el algoritmo de este
proceso se resume en la Figura 3-3.
Una vez que la elevación de la superficie de agua inicial ha sido determinada, se sigue el algoritmo
para la determinación del lecho sintético (Ver Figura 2-9 Algoritmo de Corum para obtener un lecho
sintético).
Los datos hidrometeorológicos, luego de ser obtenidos (Tabla 3-2) se someten a un análisis de calidad
estadística, como se describe en Métodos estadísticos para la evaluación de la calidad de los datos,
pág. 28.
Una vez que han sido evaluados y en su caso corregidos, los datos hidrometeorológicos (Temperatura
máxima, temperatura mínima y precipitación) son colocados en los formatos correspondientes para
su análisis en CEQUEau.
43
Figura 3-3 Algoritmo para cambiar los valores nulos del DEM por los correspondientes valores de la
elevación de la superficie de agua al momento del levantamiento LiDAR
44
archivo para los datos hidrometeorológicos Temperatura máxima, Temperatura mínima y
Precipitación, en el formato mostrado en la Figura 3-5.
45
Para obtener éstos archivos, se elaboró una herramienta que pide la dirección de una carpeta que
contenga tres subcarpetas de nombres Pres, Tmax y Tmin (Ver Figura 3-6) que contengan los datos
obtenidos desde CLICOM, y relacionados entre ellos con el nombre de la estación al que pertenecen
(P.Ej. éstos tres datos obtenidos para la estación 7005, deben tener en cada carpeta, un nombre que
inicie por 7005). Entonces los datos son entrelazados en un mismo archivo con formato de CEQUEau
mediante el algoritmo mostrado en la Figura 3-7.
Figura 3-7 Algoritmo para transformar datos en formato CLICOM a formato CEQUEau
46
3.2.3 Datos hidrométricos
Los datos hidrométricos, luego de ser obtenidos (Tabla 3-2) se someten a un análisis de calidad
estadística, como se describe en Métodos estadísticos para la evaluación de la calidad de los datos,
pág. 28.
Por otro lado, los archivos en CEQUEau que contienen la información hidrométrica, tienen el mismo
formato que se muestra en Figura 3-5, excepto que todos los datos de la matriz corresponden a los
datos hidrométricos diarios, mientras que los datos hidrométricos diarios obtenidos de BANDAS
tienen el formato mostrado en la Figura 3-8. Para colocar los datos de BANDAS en el formato de
CEQUEau, se sigue el algoritmo mostrado en la figura
Figura 3-9 Algoritmo para cambiar los datos de formato BANDAS a formato CEQUEau
47
Además, para el análisis hidrométrico de avenidas, los datos hidrométricos máximos anuales son
sometidos a un ajuste estadístico como se describe en Análisis Estadístico de Datos, pág. 32 para
obtener los caudales picos esperados de acuerdo con un periodo de retorno.
Los datos de uso de suelo (Ver 3.1 Recopilación de datos) son empleados tres procesos importantes
dentro de la investigación.
1. En la creación del archivo de datos fisiográficos para usarse en CEQUEau (Ver 2.3 Modelo
CEQUEau).
2. En la determinación del coeficiente de fricción de Manning para cada punto de las
secciones transversales que se usarán en la simulación de avenidas en 1D usando HEC-
RAS.
3. En la determinación del coeficiente de fricción de Manning sobre la zona de estudio, que
podrán ser accedidos desde cada celda de la malla generada.
Los datos de uso de suelo obtenidos de CONABIO están en formato vectorial y abarcan toda la
república mexicana, lo que hace necesario realizar una transformación de estos datos a formato ráster,
recortar los datos a la zona de estudio y transformar los valores de las claves de uso de suelo a su
correspondiente valores n de Manning. En general, el tratamiento que reciben los datos de uso de
suelo para incorporarlos dentro de las diferentes simulaciones se describe en las Figuras 22, 23 y 24.
48
Figura 3-10 Algoritmo general para incorporar los datos de uso de suelo a la simulación de caudales de
CEQUEau
Figura 3-11 Algoritmo general para preparar los datos de uso de suelo para su incorporación al modelo
de simulación de avenidas 2D de HEC-RAS
Como se ha mencionado en 2.1 Modelo y software HEC-RAS (V. 5.x), mientras que para un modelo
puramente 1D es posible realizar una simulación de flujo estable, donde el tiempo no es una variable
y que permite conocer una llanura de inundación correspondiente a un caudal único (usualmente, el
caudal pico de un evento de avenida, pues se corresponde con la llanura de inundación más
desfavorable).
49
Figura 3-12 Algoritmo general para preparar los datos de uso de suelo para su incorporación al modelo
de simulación 1D de HEC-RAS
3.3 Generación de hidrogramas sintéticos por el método de la SCS distribuidos a lo largo del
eje del río principal
Sin embargo, tanto los modelos de simulación 2D, como en las simulaciones combinadas 1D/2D, es
forzoso que sean simulaciones de flujo no estable, es por esto por lo que no basta con conocer el
caudal pico del evento, sino que es necesario conocer el hidrograma correspondiente.
Los puntos donde se elaborarán hidrogramas sintéticos se distribuyen a lo largo del eje del río en las
mismas posiciones que las estaciones ficticias que se simulan en CEQUEau (Ver 3.4.1 Simulación
de caudales medios usando CEQUEau).
50
El procedimiento que se utiliza para la elaboración de estos hidrogramas es el siguiente
1. Se simulan los caudales medios en una serie de estaciones ficticias a lo largo del eje del río
(Ver 3.4.1 Simulación de caudales medios usando CEQUEau).
2. De acuerdo con el estudio estadístico realizado sobre los datos hidrométricos (Ver 2.8.2
Análisis estadístico de los datos) se determinan los caudales para cada periodo de retorno.
3. Se buscan dentro de los resultados de las simulaciones de caudales medios realizadas con
CEQUEau caudales totales (correspondientes a la estación hidrométrica en la desembocadura
de la cuenca) similares al caudal de diseño para el periodo de retorno que se está utilizando.
4. Usando el registro de resultados obtenido en el punto anterior, determinar cuál es la
distribución del gasto total (en porcentaje) que hay entre todas las estaciones ficticias.
5. Basados en los porcentajes del gasto total que corresponden a cada estación ficticia (como se
enunció en el punto anterior), se divide el escurrimiento correspondiente al periodo de retorno
entre los puntos que espacialmente corresponden a la ubicación de las estaciones ficticias
utilizadas en CEQUEau. Se denominará en adelante estos puntos simplemente como
estaciones ficticias, dado que sobre éstos se realizarán varios cálculos más.
6. Apoyados con las herramientas de Integrated Water Management (Ver 2.6 Herramientas
Integrated Water Management), para cada estación ficticia se determina su área de aportación
correspondiente (solo la propia de cada estación ficticia, no se usa el área de aportación
acumulada hasta las estaciones ficticias, ver Figura 3-13) y los datos hidro-fisiográficos de
éstas (área, tiempo de concentración, longitud del cauce principal, pendiente media, etc.).
7. Se realiza una simulación en HEC-RAS del gasto medio correspondiente al periodo de
retorno utilizado, para determinar en cada sección del río principal (cada tramo entre dos
estaciones ficticias) el tiempo de escurrimiento sobre este.
8. Con los datos obtenidos en los puntos 5, 6 y 7 de esta lista, se elaboran hidrogramas de
acuerdo con lo establecido en 2.9 Hidrograma sintético de la SCS.
9. Utilizando la desembocadura de la cuenca de estudio como marco de referencia de tiempo,
se agrega a cada hidrograma obtenido en el punto anterior un tiempo de retraso,
correspondiente a al tiempo de escurrimiento sobre el río principal desde la estación ficticia
correspondiente, hasta la desembocadura de la cuenca de estudio (Ver Figura 3-14).
51
a) Área de aportación acumulada para la estación 2
Estación 1
Estación 2
Área de aportación de estación 1
Área de aportación de estación 2
Figura 3-13 Ejemplo de proceso de obtención de área de aportación exclusiva para una estación ficticia particular (Ref.
UTM-15N)
52
Estación Hidrométrica
Estación Ficticia
Río Principal
Río Principal Sub Cuenca Est. Fic.
2) Se determina el tiempo Cuenca de estudio
de escurrimiento hasta la
Subcuenca de estación ficticia
desembocadura.
1) Se calcula el hidrograma
3) Se agrega el tiempo de
de acuerdo a las secciones
escurrimiento como
2.9 y 3.3.
retardo al inicio del
hidrograma.
2 2
Tiempo de
1.5 1.5 escurrimiento
1 1
0.5 0.5
0 0
0 100 200 300 400 500 600 0 200 400 600 800
-0.5
Figura 3-14 Esquema de cálculo de hidrogramas sintéticos por método de la SCS con referencia de tiempo en la
desembocadura
53
3.4 Simulaciones iniciales
De acuerdo con lo establecido en 3.2 Tratamiento de datos, una vez que la información ha sido
revisada, corregida y convertida al formato adecuado, es posible hacer que el programa CEQUEau
opere y nos permita realizar simulaciones sobre los caudales máximos.
Ya que CEQUEau tiene la capacidad de simular estaciones ficticias en cualquiera de los cuadros en
que se ha discretizado el terreno (Morin y Paquet., 2007), se hace un contraste entre las direcciones
generadas por IDIRIS-CEQUEau y el eje del río principal (ambos obtenidos como se describe en
3.2.1 Datos batimétricos), para determinar en qué cuadros de la discretización es factible colocar una
estación ficticia.
Los resultados obtenidos de estas simulaciones servirán más adelante para determinar la distribución
del gasto pico a lo largo del cauce principal, utilizando los porcentajes de gasto que en cada estación
corresponden para situaciones similares de escurrimiento.
Las simulaciones iniciales en 1D utilizando HEC-RAS tienen por prerrequisitos los siguientes:
Una vez que se tiene esta información, se procede conforme al algoritmo presentado en la Figura 3-15
54
Figura 3-15 Algoritmo para realizar las simulaciones hidrodinámicas iniciales de avenidas en 1D
Figura 3-16 Inundación en sección con canales adyacentes sin puntos de dique señalados
55
a) b)
c) d)
Figura 3-17 a) Diques mal colocados (por incluir los canales adyacentes) b) Diques bien colocados (por aislar el
canal principal de los adyacentes) c) Diques mal colocados (por no considerar los canales adyacentes en situación
de inundación) d) Diques bien colocado (por considerar los canales adyacentes en situación de inundación)
Cuando existen canales adyacentes al canal principal (ya sean naturales, artificiales, u estructuras
como carreteras que debido a la resolución del DEM utilizado (Tamaño de pixel de 15 m, Ver 3.1
Recopilación de datos) se ven reflejadas en la batimetría), es necesario señalar el canal principal por
medio de puntos de bancos y de ser necesarios, puntos de dique. Esto se debe a que, si no son
colocados, el software de HEC-RAS asignará flujo a todas las zonas de la sección que estén por
debajo de la elevación de la superficie de agua (Ver Figura 3-16).
56
Los puntos de dique señalan coordenadas dentro de la sección, que indican que HEC-RAS no deberá
distribuir flujo al resto de la sección mientras que la elevación de la superficie de agua no sobrepase
la elevación de esos puntos.
Sin embargo, cuando hay varios canales inundados, es necesario delimitar con estos puntos el
conjunto de canales que pueden ser afectados por la inundación. Así que, en tales zonas, puede ser
necesario cambiar la ubicación de estos puntos para que las zonas inundadas sean representativas.
Como se aprecia en la Figura 3-17, para una misma sección, la posición adecuada de los puntos de
dique puede variar con un cambio relativamente pequeño en el tirante.
La información de entrada de caudales al sistema se hace por medio de los hidrogramas sintéticos de
la SCS (Ver 2.9 Hidrograma sintético de la SCS), ya sea con los caudales pico para el análisis 1D de
flujo estable, o bien el hidrograma completo distribuido a lo largo del eje (Ver 3.3 Generación de
hidrogramas sintéticos por el método de la SCS distribuidos a lo largo del eje del río principal) para
el flujo no estable en 1D. Cada uno de estos datos es aplicado a las secciones transversales que le
corresponde.
Todas las secciones transversales del proyecto son candidatas para recibir datos, sin embargo, la
primera y la última sección son las únicas que los requieren obligatoriamente. En general, los
elementos de entrada de datos y los datos que pueden y/o deben contener se resumen en la Tabla 3-3.
Finalmente, antes de realizar la simulación es importante determinar los parámetros que serán
utilizados para la simulación (Ver Tabla 3-4).
Como se puede apreciar en el algoritmo general que utiliza HEC-RAS para llevar a cabo las
simulaciones de avenidas en 2D (Ver Figura 2-5) para determinar el sentido del flujo de una celda a
otra, el programa calcula la continuidad de momento y de volumen en todos los elementos de cálculo
(celdas).
57
Esto ocasiona que, si se proyecta una malla demasiado amplia sobre el área de estudio, es probable
que se eleve innecesariamente la cantidad de operaciones matemáticas para cada paso de tiempo, pues
aún si las celdas no llegan a contener volumen de agua (Ver Figura 3-18), de cualquier forma tienen
que ser evaluados en ese paso del algoritmo.
Tabla 3-3 Elementos de entrada de datos al sistema 1D de simulación y datos requeridos por éstos
Simulación Elemento de Entrada Datos de entrada
Flujo inicial de escurrimiento (m3/s)
- Elevación de superficie de agua conocida
Primera sección (Aguas
- Profundidad crítica
arriba)
- Pendiente línea de energía
- Curva de aforo
- Elevación de superficie de agua conocida
Última sección (Aguas - Profundidad crítica
abajo) - Pendiente línea de energía
Flujo estable 1D - Curva de aforo
Nuevo gasto que se empleará a partir de esta
sección hasta una nueva sección de cambio de
flujo o hasta el final del tramo del río (en
Sección de cambio de
dirección aguas abajo, m3/s)
flujo (opcional,
- Elevación de superficie de agua conocida
Intermedio)
- Profundidad crítica
- Pendiente línea de energía
- Curva de aforo
-Hidrograma de nivel de agua
Primera sección (Aguas
-Hidrograma de caudal
arriba)
- Hidrograma nivel/caudal
-Hidrograma de nivel de agua
-Hidrograma de caudal
Última sección (aguas
- Hidrograma nivel/caudal
abajo)
-Curvas de aforo
Flujo no estable 1D
-Pendiente de la línea de energía
-Hidrograma de entrada de flujo lateral
-Hidrograma de entrada de flujo lateral
Sección de condiciones
distribuido
de frontera internas
- Hidrograma de elevación de flujo subterráneo
(Intermedio)
-Hidrograma nivel/caudal de condiciones
internas
Entrada obligatoria
Dato obligatorio si se señaló una entrada
Obligatorio ingresar uno de los datos listados si se señala
una entrada
58
Tabla 3-4 Parámetros principales para los diferentes tipos de simulación
Tipo de simulación Parámetro Opciones
Régimen de flujo -Subcrítico
Flujo estable 1D -Crítico
-Combinado
Intervalo de cómputo - Se puede elegir entre
0.1 segundo a 1 día
Periodo de simulación - Fecha y hora de inicio de la
simulación
- Fecha y hora de fin de la
simulación
Pasos de calentamiento 1 a 100,000
Máximo número de iteraciones 0 a 40
Flujo no estable 1D Tolerancias Elevación S.A.: 0 a 0.06 m
Tol. Aborto: 0 a 30 m
59
Figura 3-18 Ejemplo de malla 2D demasiado amplia respecto a la planicie de inundación esperada
Por otro lado, una malla demasiado pequeña podría no llegar a cubrir la totalidad de la inundación
(Ver Figura 3-19).
Figura 3-19 Ejemplo de malla 2D demasiado pequeña respecto a la planicie de inundación esperada
60
Por esa razón, es deseable contar con una simulación de avenidas en 1D previa al establecimiento del
área de terreno que será cubierta con una malla para la simulación 2D, para determinar las zonas en
condición de inundación y la extensión de ésta (Ver Figura 3-20).
a) b)
Figura 3-20 a) Ejemplo de planicie de inundación obtenida con una simulación 1D b) Ejemplo de planicie
de iundación obtenida con una simulación 2D
Una vez que se ha elegido el contorno que debe llevar la malla (auxiliados por los resultados de las
simulaciones 1D previas), el tamaño de celda y parámetros como el coeficiente de Manning, lo que
prosigue es el ingreso de los datos de flujo dentro de la malla.
Los datos son introducidos a la malla por medio de líneas de Condiciones de Frontera (CF, Boundary
Conditions Line) colocadas a lo largo de las celdas que deben ser afectadas (Ver Figura 3-21). Éstos
datos corresponden a los hidrogramas sintéticos calculados por el método del diagrama adimensional
de la SCS, que fueron obtenidos de acuerdo a lo discutido en 2.9 Hidrograma sintético de la SCS y
3.3 Generación de hidrogramas sintéticos por el método de la SCS
61
Figura 3-21 Líneas de CF colocadas en las celdas pertinentes de la malla
Finalmente, antes de realizar la simulación es importante determinar los parámetros que serán
utilizados para la simulación (Ver Tabla 3-4).
Una vez realizadas las simulaciones, la capa de resultados de tirante hidráulico en cada celda se
exporta a formato geotiff. Como se mencionó en 2.4 Matlab, los archivos en dicho formato pueden
ser importados a manera de matrices dentro de Matlab.
Ya que éstas matrices representan las manchas de inundación, una comparación entre ellas es posible
mediante el simple proceso de encontrar el valor absoluto de la resta de ambas matrices.
El archivo geotiff correspondiente a los tirantes de la mancha de inundación que se obtiene a partir
de HEC-RAS marca los pixeles “secos” como valores nulos (NaN). Así que también es posible
encontrar zonas que en una de las matrices de inundación contiene celdas mojadas que corresponden
a celdas secas en la otra, señalando las partes de la matriz donde hay diferencia en el área que abarca
la mancha de inundación.
62
Aplicación
Para la presente investigación, se utilizará el río La Sierra como caso de estudio, en la cuenca
correspondiente a la estación hidrométrica “Pueblo Nuevo”, en el estado de Tabasco, mientras que
la cuenca está distribuida entre los estados de Tabasco y Oaxaca (Ver
Figura 4-1)
Esta estación se encuentra en las coordenadas geográficas, Long. -92.879, Lat. 17.854, con número
de identificación de estación hidrométrica a nivel nacional 30016.
(-92.24°,17.86°)
(-93.07°,17.86°)
(-92.24°,16.77°)
(-93.07°,16.77°)
Estación hidrométrica
Escurrimientos principales
63
La cuenca del Río la Sierra está ubicada entre las longitudes -93.07° E y -92.24° E, y las latitudes
16.77° N y 17.86° N.
Las propiedades fisiográficas correspondientes a la cuenca en estudio fueron obtenidas con las
herramientas “Integrated Water Management” (Ver 2.6 Herramientas Integrated Water Management)
a partir de la información que se detalla en 3.2.1 Datos batimétricos y se resumen en la Tabla 4-1.
Tabla 4-1 Parámetros fisiográficos de la cuenca de la estación Hidrométrica Pueblo Nuevo (30016)
Parámetro Valor
Tiempo de Concentración (Kirpich) 28.18 hrs
Área 4,366.13 km2
Perímetro 610 km
Pendiente media 1.99%
Lon. Cauce Principal 244.04 km
En la zona de estudio se encuentran las comunidades de: Pueblo Nuevo de las raíces, La Agraria, La
Huasteca, Hueso de Puerco, Nuevo Progreso, San Antonio, Puyacatengo, Payacatengo, El Perú,
Ceibita, Morelia, Tacotalpa, Teapa, Las Nieves, Graciano Sánchez, Tapijulapa, Solosuchiapa, La
Gloria, Francisco I. Madero, San Lorenzo, Amatán, Oxolotán, Ixhuatan, Chapayal Grande, Tapilula,
Rayón, Zacatonal de Juarez, Los Naranjos, Huitiupan, Simojovel, Carmen Zacatal, El Bosque, San
Pedro Nishtalucum, Zizim, Pueblo Nuevo Sitalá, Esquipulas, Tzajalá, Sitalá, Guaquitepec, San Juan
Cancuc, Larrainzar, Chalchihuitán, San Cayetano, Larrainzar, Tenejapa, Tenango, Pantheló, la mayor
parte de éstas teniendo menos de 2 Km2, y en general tienen distribuciones amplias entre unas y otra
(Ver Figura 4-2).
64
Estación hidrométrica
En este proyecto se utilizan las cartas E15D11, E15D21, E15D31, E15D41, E15D51, E15D12,
E15D22, E15D32, E15D42 y E15D52 (cada una con su subdivisión A1-A4, B1-B4, C1-C4, D1-D4)
con una resolución de 1:10,000 (Ver Figura 4-3). Se hace notar, sin embargo, para la zona de estudio,
solo se cuenta con datos de alta resolución para las cartas E15D11, E15D12, E15D21 y E15D22.
Para el resto de la zona de estudio, se cuenta con una resolución de pixel de 15 m (Ver Figura 4-4),
que fue obtenida del CEM (Continuo de Elevación Mexicano) V.3.0 proveniente de la misma página
de INEGI.
65
Figura 4-3 Cuadrícula de correspondencia geográfica de las cartas topográficas disponibles en INEGI
66
Elevación (m)
Parteaguas
Res. 5 x 5 m
Res 15 x 15 m
Por las limitaciones que tiene el método de levantamiento LiDAR, que se han mencionado en 2.7
Método de obtención de batimetría sintética propuesto por Corum, en los cuerpos de agua no hay
retorno de datos batimétricos. Éstos valores son señalados con un marcador de valores nulos. Sin
embargo, cada carta batimétrica unida para formar la totalidad del terreno de estudio tiene su propio
marcador de valores nulos (regularmente un valor que no aparece dentro de la carta), cuyo valor puede
ser un valor legítimo en otra de las cartas.
Es por esto que la delimitación de valores nulos se tiene que hacer carta por carta, para después unirse
en una sola discontinuidad completa a lo largo del modelo de elevación completo usado en la zona
(Ver Tabla 4-2 y Figura 4-5).
67
Tabla 4-2 Marcadores de valores anómalos, desde aguas abajo hasta aguas arriba de la zona de la
cuenca
Marcadores de valores anómalos
-8
6
151
197
213
399
440
444-49
452
490
525
Elevación (m)
Parteaguas
Valores nulos
Figura 4-5 Acercamiento a los valores nulos correspondientes a los datos faltantes del levantamiento
batimétrico
68
Una vez que han sido detectados los valores nulos de todas las cartas, se unifica el valor numérico
para datos faltantes con el valor de -10 (Puede realizarse tal operación desde algún SIG o
matricialmente utilizando MATLAB. Por facilidad de uso, se utilizó IDRISI).
Tales valores nulos dentro de la batimetría son un faltante crítico dentro de las simulaciones debido
representan justamente las áreas principales de simulación. Es por esto, que de acuerdo con el método
de Corum (2.7 Método de obtención de batimetría sintética propuesto por Corum) se debe obtener un
lecho sintético equivalente para poder realizar sobre éstas las simulaciones.
El primer paso de dicho método es conocer la elevación del espejo de agua. Para lograr esto se siguió
el algoritmo descrito en la Figura 3-3.
Para “cerrar” el DEM, sustituyendo los valores nulos por las elevaciones del espejo de agua al
momento del levantamiento, se inició importando los archivos Geotiff del DEM y de uso de suelo
haciendo uso de la herramienta IMPTIF (Ver 6.1.1 IMPTIF) dentro de MATLAB.
Haciendo uso de las herramientas TopoToolBox obtenemos un eje del río principal temporal ya que
en esta etapa del proceso de “cierre” del DEM solo necesitamos un eje continuo que esté dentro de
los valores nulos. En el caso de estudio no fue necesario, pero podría forzarse al eje a pasar dentro de
los valores nulos señalando los valores nulos como un número más negativo (hemos mencionado que
se utiliza el valor de -10 para señalar los valores nulos, pero nada impide usar -9999 o cualquier otro
valor mientras que sea lo suficientemente más bajo que la elevación del terreno).
En general, el proceso de obtención del eje del río se realiza mediante una función elaborada para tal
propósito denominada CALRIO (Ver 6.1.2 CALRIO), que sigue el algoritmo descrito en la Figura
4-8.
La dirección de cada pixel del río, que se menciona en la Figura 4-8, se obtiene mediante una función
elaborada para tal propósito denominada DIRECCION (Ver 6.1.12 DIRECCIÓN), que obtiene las
coordenadas de uno de los pixeles del eje del río principal (Arreglo RÍO, variable RIOCAL, Ver Tabla
4-3) y del pixel siguiente y determina su dirección (Ver Figura 4-6).
69
8 1 2
7 PA 3
6 5 4
Figura 4-6 Direcciones posibles del siguiente píxel del río principal a partir del píxel actual (PA)
Se procede entonces para cada pixel del eje principal del río a obtener una sección transversal, que
tendrá un ángulo acorde a las direcciones obtenidas de los pixeles (Ver Figura 4-7) usando la función
SECCIONAR (Ver 6.1.16 SECCIONAR). Procedemos a crear una matriz que contenga los datos de
las secciones transversales obtenidas de cada píxel (Arreglo RÍO, variable SECRIO, Ver Tabla 4-3)
8 1 2 8 1 2 8 1 2
7 PA 3 7 PA 3 7 PA 3
6 5 4 6 5 4 6 5 4
Figura 4-7 Dirección de seccionamiento según la dirección del píxel del eje principal del río
Figura 4-8 Algoritmo de extracción de datos de los pixeles del eje del río
Tabla 4-3 Variables usadas en el cómputo de la batimetría (Primer iteración)
Variable o arreglo Variables dentro del Descripción
arreglo
MAT Matriz con los valores numéricos de las
DEM
elevaciones
70
REF Arreglo con las referencias geográficas
del archivo Geo TIFF que se importó
INFO Arreglo con las características de la
imagen Geottif que se importó
MATPR Es la matriz MAT que ha sido sujeta a un
proceso de limpieza de valores anómalos
con la herramienta TopoToolBox.
SECERRADO Matriz de secciones transversales con el
espejo de agua sustituyendo los valores
nulos (Ver Figura 4-9)
DEMCERRADO La matriz de elevaciones con la primera
iteración del espejo de agua
implementado
DEMCERRADOPR La matriz DEMCERRADO que ha sido
sujeta a un proceso de limpieza de
valores anómalos con la herramienta
TopoToolBox
DFLUJO Matriz de datos de dimensiones
correspondientes a la matriz del DEM
que contiene las unidades de flujo en los
pixeles del río principal y 0 en los demás
pixeles (Ver Figura 2-8)
RIOPRIN Matriz de datos de dimensiones
correspondientes a la matriz del DEM
que tiene valores de 1 en los pixeles del
RÍO
río principal y 0 en los demás pixeles
RIOCAL Matriz con número de filas igual al
número de pixeles en el río principal, con
datos de coordenadas, direcciones, y
jerarquía de cada píxel
SECRIO Contiene los datos de las secciones
transversales de cada píxel tomadas
conforme a la dirección de éste.
La primera iteración del espejo de agua se realiza uniendo los valores válidos que limitan los valores
nulos de cada sección transversal. Ya que es posible que los valores en las orillas sean diferentes, se
toma el menor de los dos (Ver Figura 4-9). Sin embargo, no se puede tomar los resultados de este
proceso como el espejo definitivo, principalmente porque, aunque los resultados son buenos en
tramos rectos, en los meandros debido a la curvatura existen pixeles que no pueden ser incluidos
dentro de las secciones, lo que genera agujeros (Ver Figura 4-10).
71
Sección Transversal con valores nulos
14
9
Elevación (m)
4
Valores nulos
-1
0 100 200 300 400 500 600 700
-6
-11
Distancia (m)
Se cambian los valores
nulos por la menor de las
elevaciones de las orillas
del espejo de agua
6.5
6
5.5
5
4.5
4
0 100 200 300 400 500 600 700
Distancia (m)
Figura 4-9 Método para cambiar los valores nulos por las elevaciones de los espejos
de agua en las secciones transversales
Figura 4-11 Tratamiento de los pixeles con valor nulo que no se procesaron dentro de las secciones
Puede apreciarse en la Figura 4-12, que mientras en los tramos más rectos del río y curvas suaves de
los meandros los resultados son buenos, en los meandros más pronunciados comienzan a notarse los
efectos del tratamiento descrito en la Figura 4-11 y se muestran discontinuidades en el fondo del
río.
73
Elevación (m)
Parteaguas
Corte transversal
a) b)
74
Figura 4-12 Resultados primer iteración del espejo de agua
Para la segunda iteración del espejo de agua, tomamos el espejo obtenido en la primera iteración y se
repite el proceso de limpieza de valores anómalos y trazado del eje del río principal con
TopoToolBox. El trazado de este eje de río mejor concebido podría considerarse la función principal
de haber obtenido la primera iteración del espejo del río y se usará por el resto del proceso, así que
será llamado eje definitivo.
Sobre el eje definitivo se realiza un nuevo seccionamiento sobre la matriz de elevaciones original
(Arreglo DEM, variable MAT, Ver Tabla 4-3), aunque en esta ocasión no se realizará pixel por pixel
porque queremos evitar que las secciones se traslapen unas con otras.
Este seccionamiento es un paso previo para la interpolación vectorial, por lo que las secciones que se
obtendrán necesitan ser tan pequeñas como se pueda y que abarquen la totalidad del espejo de agua
y que estén tan juntas como se pueda sin que se traslapen entre ellas, el algoritmo se detalla en la
Figura 4-13.
Figura 4-13 Algoritmo seccionado sin traslape y sin cortar espejos de agua
Se hace notar que a diferencia del método descrito en la Figura 4-7, que se realiza en su totalidad en
coordenadas matriciales, puesto que al hacerse pixel por pixel solo existen 8 posibles direcciones y
cada punto de la sección corresponde a un pixel de la matriz numérica del DEM, este método requiere
75
obtener los puntos de la sección a través de la pendiente del vector unitario, lo que se hace en
coordenadas UTM, que después pueden transformarse a coordenadas matriciales (Ver ecuaciones
Ecuación 2-16 y Ecuación 2-17).
b) Un solo espejo de
agua (Sección
habitual. No hay
batimetría en el
c) Dos espejos de
agua en la misma
dirección
d) Dos espejos de
agua en direcciones
diferentes
e) Dos espejos de
agua en cuerpos de
agua diferentes
Figura 4-14 Situaciones posibles del espejo de agua dentro de las secciones transversales
76
Una vez que se han obtenido estas secciones, se procede a sustituir los valores nulos de las secciones
con el menor de los valores de elevación de las orillas del espejo de agua (Ver Figura 4-9)
Ahora bien, en cuanto al espejo de agua dentro de una sección pueden presentarse cinco situaciones
(Ver Figura 4-14). Que la sección no contenga espejo de agua (como en las zonas altas, donde el
tirante es somero y el levantamiento LiDAR alcanzó a muestrear el fondo), que contenga un espejo
de agua, que contenga más de un espejo de agua que tengan el mismo sentido de flujo (Islas dentro
del río), que contenga más de un espejo de agua que tengan diferente sentido de flujo (La sección
interseca más de una vez el eje del río) o que contenga más de un espejo de agua y que pertenezcan a
cuerpos de agua diferentes (Canales adyacentes).
Existe gran importancia en el hecho de diferenciar cada caso, puesto que se podría incurrir en errores
tales como nivelar en igualdad de condiciones dos espejos de agua en una sección que no tienen el
mismo flujo o no pertenecen al mismo cuerpo, o bien, ignorar uno de los dos lados del mismo espejo
de agua dividido por una isla fluvial.
Se desarrolló una herramienta que permite identificar estos casos y delimitar el espejo de agua
principal (aún si está separado, se considera uno) dentro de una sección transversal, denominada
ESPBOOL (Ver 6.1.13 ESPBOOL)
Esta herramienta primeramente crea una lista de las secciones que no fueron omitidas por distancia
mínima o por traslape (Segundo seccionamiento sin traslape, algoritmo mostrado en la Figura 4-13,
Arreglo SECCESPEJO, variable LISTSECCTRASLA ), después, a partir de la matriz de secciones
transversales creada por ese mismo algoritmo, se crea un vector a partir de la columna de elevaciones,
en el que las filas correspondientes a valores nulos se marcan con un 1 y todas las demás se marcan
con un 0.
Por cada elemento de la lista de secciones mencionada, se aíslan las filas de datos de secciones
transversales que le corresponden dentro de la matriz (SECCSESPEJO.SECCSNULAS, Ver Tabla
4-4) después se cuentan los “Inicios de espejos de agua” todos los valores que su elemento actual sea
1 y el anterior 0 y se cuentan los “Finales de espejos de agua”, todos los valores que su elemento
actual sea 0 y el anterior 1. El número de espejos de agua dentro de la sección se toma como el mayor
de estos dos.
77
Si no hay espejos de agua (Figura 4-14, a), no se realiza ninguna acción y se continúa a la siguiente
sección de la lista SECCSESPEJO.LISTSECCTRASLA, en el caso de haber un espejo de agua
(Figura 4-14,b), se copia la información del vector booleano SECCSESPEJO.SECCSNULAS en la
columna del espejo principal.
Si hay más de un espejo de agua, se realiza una búsqueda de colisiones con el eje del río. Si hay más
de una colisión (Figura 4-14,d) , se toma como espejo principal únicamente el espejo más cercano al
píxel del eje en la sección transversal. Si sólo hay una colisión, se buscan valores nulos (indicadores
de espejo de agua) hacia arriba y hacia abajo del píxel del eje del río en la sección transversal. Los
espejos que se alejan más de la tolerancia ingresada por el usuario no se consideran parte del mismo
cuerpo de agua (Figura 4-14, e), de los restantes se copia su información booleana dentro del vector
booleano (SECCSESPEJO.SECCSNULAS, Tabla 4-4) del espejo principal (Figura 4-14, c).
“El objetivo de convertir DEM’s a TIN’s es obtener el menor número de puntos de elevación como
sea posible mientras se obtenga la mayor cantidad de información que se pueda acerca de las
estructuras topográficas. Desafortunadamente, es usualmente imposible alcanzar ambas metas
simultáneamente. Ya sea que más información estructural sea representada por más puntos de
elevación o que una red más pequeña provea una representación estructural más generalizada.
Cuando se convierten DEM’s a TIN’s, este intercambio tiene que llevarse a cabo, encontrando un
punto intermedio, pero ya que éste difiere entre un método y otro y entre un tipo de superficie y otra,
experimentos de prueba y error son usualmente necesarios para obtener mejores resultados”
En el presente trabajo se pudo constatar esa afirmación y es por esto qué, en este paso, fue necesario
complementar la inspección y corrección de datos (mediante el arreglo de las fronteras) manualmente.
Para tal propósito, se importaron los datos vectoriales de las orillas del espejo de agua al programa
CIVIL3D y se ajustó la triangulación del TIN por medio de fronteras y límites.
Por un lado las líneas de frontera (Boundary lines), que son polígonos cerrados, limitan la
triangulación al interior de este polígono (Ver Figura 4-18), mientras que las líneas de límite
(Breaklines) indican que los vértices de los triángulos no deben cruzar a través de ellas (Ver Figura
4-19), es decir, evitan que un punto correspondiente al terreno genere un vértice con un punto
correspondiente al espejo de agua.
Las líneas de frontera evitan que secciones cercanas, pero no continuas, interpolen entre ellas. Las
líneas límite ayudan a corregir la pérdida de información debido a la discretización de las curvas.
79
La importación de los puntos hacia CIVIL3D se hace en coordenadas (X,Y,Z). El procedimiento de
importación de coordenadas es como sigue:
80
Figura 4-16 Ventana de importación de puntos hacia la superficie TIN en CIVIL3D
Se procede a la importación de las líneas de frontera y las líneas de límite de los espejos de agua.
Ya que dentro del programa estos elementos deben ingresarse como líneas, se ha creado un script de
AutoCAD que importe las coordenadas. Es por esto, que se crearon archivos *.scr con la instrucción
_pline y las coordenadas de los vértices de las líneas en orden (Ver Figura 4-17). Para terminar con
la importación de las líneas, se ejecuta el script correspondiente a cada una (consola run script de
CIVIL3D).
81
Eje del río
Orillas del límite de las secciones
Elevación (m)
Parteaguas
Orilla de la superficie TIN
Aristas de la triangulación
0km 20 km
Elevación (m)
Figura 4-18 Aplicación de las líneas de frontera a los datos de triangulación TIN
82
Eje del Río
Orillas del límite de las secciones
Parteaguas
Elevación (m) Orilla de la superficie TIN
Curvas de nivel
0km 20 km
Elevación (m)
Figura 4-19 Aplicación de las líneas de límite a los datos de triangulación TIN
83
Parteaguas
Elevación (m) Corte transversal
Elevación (m)
Figura 4-20 Comparación de la primera iteración del espejo de agua con la segunda
Finalmente, la información de la malla TIN que ha sido corregida es incorporada dentro de la matriz
de datos numéricos dentro de MATLAB, el cual se exporta como un DEM usando la función
EXP2TIFF (Ver 6.1.3) . Los datos obtenidos para el espejo de agua se consideran adecuados para
seguir adelante con el proceso de elaboración de lecho sintético (Ver Figura 4-20).
84
El método de Corum para determinar un lecho sintético requiere de dos elementos a) La elevación
del espejo de agua y b) El caudal que causó que existiera ese espejo de agua.
Con la elevación del espejo de agua, el caudal se determina a partir de los metadatos de las cartas
batimétricas. Los registros del limnígrafo en la estación 30016 (Contenido en el mismo archivo de
caudales que se obtuvo de BANDAS, CONAGUA, Ver Tabla 3-2 Fuente de obtención de datos
requeridos).
De acuerdo con los metadatos que acompañan las cartas batimétricas obtenidas de INEGI, la mayor
parte de la exploración LiDAR fue realizada en el año 2008 (Solamente las cartas E15D11 y E15D12
fueron levantadas el año 2007). Sin embargo, no se cuenta con una fecha exacta del levantamiento,
por lo que se deberá que determinarla.
De acuerdo con la elevación del espejo de agua que ha sido determinada, en la desembocadura se
tiene una elevación de 5.52 m (Figura 4-21).
Sección en 0+000
8.3
7.8
7.3
6.8
6.3
5.8
5.3
Elev. Sup. Agua: 5.52
4.8
0 100 200 300 400 500 600
De acuerdo con las mediciones del limnígrafo en la estación 30016 para el año 2008 y dado que
tampoco se conoce la elevación del fondo en la sección, se recurrió al encuentro de una solución
geométrica y funcional (se considera adecuado usar este año en los registros, dado que no hay datos
de escala del limnígrafo para esa estación en el año 2007, porque la mayor parte del levantamiento
85
ocurrió en el 2008 y porque el periodo es lo suficientemente cercano como para que no sea probable
que hayan existido grandes cambios en el caudal).
De la Figura 4-22 se deducen las propiedades geométricas de la sección mojada con las ecuaciones
Ecuación 4-1 y Ecuación 4-2, combinada con la ecuación de Manning en la
para determinar el caudal.
1 𝐴
𝑄= √𝑆 Ecuación 4-3
𝑛 𝑃
Donde:
A: Área en la sección en la desembocadura
P: Perímetro mojado
Q: Caudal en la desembocadura
b: ancho del fondo de canal
Y: Profundidad del canal hasta el espejo de agua conocido
DY: Cambio en el tirante por encima del espejo de agua conocido
86
Z1: Pendiente de la margen izquierda
Z2: Pendiente de la margen derecha
S: Pendiente del canal
n: coeficiente de rugosidad de Manning
Z1=12.93 m/m
Z2= 3.30 m/m
b = 98.56 m
S=0.000022
Como coeficiente de Manning, tomamos el sugerido por Chow (1994) para canales dragados rectos
y uniformes (por las imágenes satelitales se puede ver que, ya que hay comunidades en esa zona, el
canal del río ha sido alterado).
n=0.02
Ahora, dentro del archivo de información hidrométrica obtenido desde CONABIO se encuentran
curvas elevación-gasto.
Para cada valor de elevación se hace una aproximación del gasto mediante la
, usando el valor de la lectura como elevación “Y” y la diferencia con la lectura siguiente como DY.
Se compara el valor del gasto calculado con el gasto medido para cada tirante, y se toma el valor que
logre una mejor aproximación con el gasto de avenida.
Q=523.969 m3/s.
De la curva de aforo correspondiente al primer tercio del año (Transición), es decir, del registro de
lecturas del limnígrafo podemos determinar que la fecha en que se tiene este caudal es el 2 de enero
del 2008.
Figura 4-23 Distribución Espacial de la cuenca en CEQUEA (Adaptada de Mercado, V., 2010)
Los valores mostrados en las Tabla 4-6 son los valores de las estaciones y los caudales que se
distribuirán en ellas. El proceso de encontrar el lecho sintético corresponde a un lecho “erosionado”
dado a la naturaleza del proceso de ir bajando el nivel del espejo principal la diferencia entre el tirante
observado y el calculado.
Antes de comenzar con las simulaciones necesarias para la elaboración del lecho sintético, es
necesario colocar los puntos de bancos y diques. Como se ha mencionado antes, en 3.4.2 Simulación
hidrodinámica de avenidas en 1D usando HEC-RAS , la adecuada colocación de éstos puntos puede
causar enormes diferencias en los resultados.
88
Se elaboró una herramienta denominada BANYDIQS (6.1.15 BANYDIQ) que se encarga de hacer
una primera sugerencia de la posición de los puntos de bancos y diques (al final, aún se recomienda
hacer una inspección visual de los resultados para detectar si existen errores en el posicionamiento.
Tabla 4-5 Equivalencia de cuadros CEQUEau en coordenadas UTM y cadenamiento del eje del río
Iniciales X 485000 Y 1845000
UTM15
X Y X Y Cadenamiento
23 12 A 23 12 550000 1855000 244+621.31
24 12 A 24 12 555000 1855000 237+898.71
24 13 C 24 13 555000 1860000 231+026.54
24 14 B 24 14 555000 1865000 224+363.31
25 15 D 25 15 560000 1870000 215+546.74
25 17 A 25 17 560000 1880000 188+647.09
24 18 A 24 18 555000 1885000 179+055.74
23 19 D 23 19 550000 1890000 167+729.79
22 19 D 22 19 545000 1890000 157+827.08
21 19 A 21 19 540000 1890000 147+485.55
20 21 A 20 21 535000 1900000 127+334.94
19 23 A 19 23 530000 1910000 113+312.82
18 24 B 18 24 525000 1915000 104+918.12
18 25 A 18 25 525000 1920000 98+712.97
17 27 A 17 27 520000 1930000 82+791.77
17 28 A 17 28 520000 1935000 73+830.53
17 29 D 17 29 520000 1940000 68+267.11
17 30 A 17 30 520000 1945000 55+340.06
16 34 D 16 34 515000 1965000 17+262.23
15 35 B 15 35 510000 1970000 6+693.05
Para este posicionamiento inicial, se aplica una media móvil de tamaño definido por el usuario para
suavizar los datos de la sección y evitar que pequeñas protuberancias obliguen al programa a
declararlas un punto de banco o dique. Sobre esa información de elevaciones suavizadas, se determina
el punto que corresponde al eje del río y se busca hacia la izquierda y la derecha valores máximos,
señalando un valor máximo cuando el valor evaluado sea mayor que el valor siguiente y el anterior.
Para el fin del erosionado, se elaboró una herramienta que con el fin de ejecutar el proceso de Corum
para obtener lechos sintéticos. La herramienta EROSIONADOR que hace uso de otra herramienta
desarrollada denominada TIRANTES (Ver 6.1.11 TIRANTES). TIRANTES ejecuta una simulación
89
en HEC-RAS y lee las elevaciones de los fondos y la superficie de agua, para determinar el tirante en
cada sección.
Como indica el algoritmo descrito en la Figura 2-9, la diferencia entre la superficie calculada y la
observada se considera como el error que deberá erosionarse de la elevación del lecho (Ver
Figura 4-25).
El proceso se repite hasta que la superficie de agua calculada coincida con la superficie de agua
conocida. En este caso, se consideró un umbral de error de 0.006 m (Ya que es el mismo umbral de
error que tiene HEC-RAS para tirantes).
Tabla 4-6 Distribución del escurrimiento del día 02/01/2008 a lo largo del eje
No. Est. Ficticia Sección %Esc. Acum. % Esc Q (m3/s)
1 244+621.31 0.009 0.009 4.695
2 237+898.71 0.020 0.011 10.262
3 231+026.54 0.027 0.007 14.059
4 224+363.31 0.048 0.021 25.180
5 215+546.74 0.048 0.000 25.277
6 188+647.09 0.100 0.052 52.413
7 179+055.74 0.140 0.040 73.182
8 167+729.79 0.184 0.045 96.590
9 157+827.08 0.216 0.032 113.275
10 147+485.55 0.228 0.012 119.688
11 127+334.94 0.467 0.238 244.604
12 113+312.82 0.596 0.129 312.135
13 104+918.12 0.645 0.050 338.167
14 98+712.97 0.711 0.066 372.588
15 82+791.77 0.736 0.025 385.470
16 73+830.53 0.737 0.001 386.219
17 68+267.11 0.739 0.002 387.348
18 55+340.06 0.744 0.004 389.691
19 17+262.23 0.747 0.003 391.148
20 6+693.05 1 0.253 523.969
90
Elevación
(msnm)
Sección erosionada
Sección original
Elevación
(msnm)
91
4.3 Hidrogramas sintéticos por el método de la SCS
Se procede primeramente a realizar un análisis estadístico de los datos obtenidos para la estacón
30016, con los procedimientos que se indican en 2.8.2 Análisis estadístico de los datos.
92
Tabla 4-8 Propiedades generales de la serie de caudales máximos anuales de la estación 30016
Datos Generales
Longitud 58
Q1 (m3/s) 768 RCI 139.000
Q2 (m3/s) 862.5 LimInf 559.500
Q3 (m3/s) 907.00 LimSup 1115.499
Media 847.22
Desv. Est. 92.57
Coef.Var. 0.11
Coef. Asim.
Fisher 0.15
Curtosis 2.53
14
HISTOGRAMA QMAX ANUALES EST. 30016
12
10
OBSERVACIONES
0
677 - 720 720 - 763 763 - 806 806 - 849 849 - 892 892 - 935 935 - 978 978 - 1021 - >1064
1021 1064
RANGO DE GASTOS (M3/S)
Las pruebas de calidad de los datos (aleatoriedad, homogeneidad e independencia) descritas en 2.8.1
Calidad de las series de datos, fueron realizada sobre la serie y se obtuvieron los resultados
presentados en la Tabla 4-9.
La serie de datos no mostró valores anómalos y pasó las pruebas de calidad, por lo que no se realizarán
cambios adicionales sobre esta serie de datos cronológica.
93
Tabla 4-9 Resultados de pruebas de calidad sobre la serie de Qmax anuales de la estación 30016
Wald y
Prueba Hatanaca Prueba Iteración Prueba Wolfowitz
Alfa 5 Alfa 5 Alfa 5
Ns 4 Ni0 8 W 41662279.857
Ni 2 Ni1 8 Xw 41622713.155
w 2 N0 27 Sw^2 3991357503.250
e 0.304 N1 14 e 0.626
Media 20.806 Xw 19.439 Za/2 1.960
Sigma 7.866 S2w 8.039 Resultado La prueba se
Za/2 1.960 w 16.000 considera
Resultado La muestra se e -1.213 independiente
considera aleatoria Za/2 1.960
Resultado
La muestra se
considera
homogénea
Las pruebas de ajuste descritas en 2.8.2 Análisis estadístico de los datos, se realizan sobre la serie
de Qmax anuales (Tabla 4-7), con los siguientes resultados
La serie se dividió en diez rangos, con límite inferior 677 y un tamaño de rango de 43, obteniendo
la Tabla 4-10 con las observaciones por rango para la elaboración del histograma.
94
Los parámetros obtenidos para las pruebas de ajuste tabulan de la ¡Error! No se encuentra el origen
de la referencia. a la Tabla 4-16
α 83.760
β 10.115
Tabla 4-13 Parámetros de la muestra de Qmax anuales para la distribución Pearson III
Variable Valor
μ3 116647.605
C.As.Fish 0.147
α 6.806
β 185.000
λ 0.147
E -411.889
Xgog -411.889
95
Tabla 4-15 Parámetros muestra de Qmax anuales para la distribución LogPearson III
Variable Valor
μ3 116647.605
C.As.Fish 0.147
α 6.806
β 185.000
λ 0.147
E -411.889
Xgog -411.889
Tabla 4-16 Parámetros muestra de Qmax anuales para la distribución Gumbel Tipo I
Variable Valor
a 805.558
b 72.178
Los datos generales de la muestra que requieren algunas de las distribuciones pueden verse en la
Tabla 4-8.
96
Tabla 4-17 Comparación de Qmax observado contra Qmax estimado por diferentes funciones de ajuste
Qmax Obs. ProbEx D. Gamma D. Exp D. Pearson III D. Gumbel I D. LogNorm D. LogPear. III
677.00 0.0169 662.80 14.48 658.86 704.11 667.42 671.21
692.50 0.0339 686.22 29.22 683.54 717.56 689.36 691.98
696.00 0.0508 701.58 44.21 699.61 726.77 703.85 705.78
708.00 0.0678 713.42 59.48 711.95 734.10 715.08 716.52
717.11 0.0847 723.25 75.02 722.15 740.35 724.44 725.50
718.50 0.1017 731.77 90.86 730.97 745.89 732.58 733.34
723.00 0.1186 739.36 107.00 738.81 750.93 739.86 740.36
737.00 0.1356 746.27 123.45 745.93 755.60 746.50 746.78
739.00 0.1525 752.64 140.23 752.48 759.98 752.64 752.73
744.00 0.1695 758.60 157.34 758.59 764.14 758.39 758.32
747.00 0.1864 764.21 174.81 764.34 768.13 763.83 763.60
750.00 0.2203 774.63 210.87 774.99 775.69 773.95 773.47
750.00 0.2203 774.63 210.87 774.99 775.69 773.95 773.47
754.10 0.2373 779.52 229.49 779.98 779.32 778.72 778.13
767.00 0.2542 784.25 248.53 784.79 782.86 783.33 782.64
771.00 0.2712 788.83 268.01 789.45 786.35 787.81 787.03
787.12 0.2881 793.28 287.94 793.98 789.78 792.17 791.32
788.00 0.3051 797.64 308.36 798.39 793.18 796.44 795.51
795.00 0.3220 801.90 329.28 802.71 796.54 800.63 799.63
805.00 0.3390 806.08 350.73 806.94 799.88 804.75 803.69
806.00 0.3559 810.20 372.74 811.10 803.21 808.81 807.70
820.00 0.3729 814.26 395.33 815.21 806.54 812.82 811.66
831.00 0.3898 818.28 418.54 819.26 809.87 816.80 815.60
831.80 0.4068 822.26 442.41 823.27 813.20 820.74 819.50
834.90 0.4237 826.21 466.97 827.25 816.56 824.66 823.39
837.60 0.4407 830.14 492.26 831.20 819.93 828.57 827.28
840.06 0.4576 834.06 518.33 835.14 823.33 832.47 831.15
842.00 0.4746 837.97 545.23 839.06 826.77 836.38 835.04
857.00 0.4915 841.89 573.01 842.99 830.25 840.29 838.94
865.51 0.5085 845.81 601.73 846.92 833.78 844.21 842.85
868.00 0.5254 849.75 631.46 850.86 837.37 848.16 846.80
874.00 0.5424 853.72 662.27 854.82 841.02 852.13 850.78
881.00 0.5593 857.71 694.25 858.80 844.75 856.15 854.80
883.54 0.5763 861.74 727.47 862.82 848.55 860.20 858.87
884.87 0.5932 865.82 762.06 866.89 852.45 864.31 863.00
885.90 0.6102 869.95 798.12 871.00 856.46 868.49 867.19
888.13 0.6271 874.15 835.78 875.17 860.58 872.73 871.47
889.00 0.6610 882.78 916.53 883.74 869.22 881.49 880.30
889.00 0.6610 882.78 916.53 883.74 869.22 881.49 880.30
890.00 0.6780 887.24 959.98 888.16 873.77 886.02 884.89
894.87 0.6949 891.82 1005.79 892.69 878.51 890.68 889.60
900.50 0.7119 896.52 1054.22 897.34 883.45 895.48 894.47
904.00 0.7288 901.37 1105.58 902.13 888.63 900.44 899.51
906.00 0.7458 906.40 1160.26 907.09 894.08 905.59 904.75
907.33 0.7627 911.62 1218.71 912.24 899.83 910.95 910.21
920.00 0.7797 917.08 1281.49 917.61 905.94 916.57 915.94
922.00 0.7966 922.80 1349.31 923.23 912.46 922.47 921.98
922.40 0.8136 928.85 1423.03 929.17 919.47 928.71 928.38
929.00 0.8305 935.27 1503.77 935.46 927.07 935.37 935.21
932.00 0.8475 942.15 1593.04 942.19 935.38 942.51 942.56
944.00 0.8644 949.59 1692.83 949.47 944.58 950.27 950.56
950.00 0.8814 957.75 1805.96 957.43 954.90 958.80 959.37
960.00 0.8983 966.84 1936.56 966.27 966.70 968.33 969.25
979.89 0.9153 977.17 2091.02 976.31 980.53 979.21 980.57
988.57 0.9322 989.29 2280.07 988.06 997.29 992.02 993.94
990.75 0.9492 1004.17 2523.80 1002.44 1018.70 1007.85 1010.53
1021.00 0.9661 1023.95 2867.32 1021.49 1048.59 1029.03 1032.86
1100.80 0.9831 1055.22 3454.57 1051.45 1099.25 1062.86 1068.82
97
Distribuición Gamma 1) Distribuición Pearson III 2)
1200 1200
1100 1100
1000 1000
Qmax
Qmax
900 900
800 800
700 700
600 600
-2 -1 0 1 2 3 4 -2 -1 0 1 2
Z Z
Qmax
900 900
800 800
700 700
600 600
-2 -1 0 1 2 3 4 5 -2 -1 0 1 2 3 4
(X-a)/b Z
Figura 4-27 Distribución de los datos observados en comparación de los valores esperados (Pte. 1)
98
Distribuición Log-Pearson III 5)
1200.00
1100.00
1000.00
Qmax
900.00
800.00
700.00
600.00
-2 -1 0 1 2 3 4
Z
Figura 4-28 Distribución de los datos observados en comparación de los valores esperados (Pte. 2)
99
A continuación, se procede a realizar una prueba de eficiencia de ajuste por medio de la Raíz del error
cuadrático (Ecuación 4-4).
1
𝑅𝐸𝐶𝑀 = ( (𝑥 − 𝑥 ) ) Ecuación 4-4
𝑛
Donde:
RECM: Raíz del error cuadrático medio
n: número de elementos de la muestra
xi: Valor i observado
𝑥 : Valor i esperado
Tabla 4-18 Valores (𝒙𝒊 − 𝒙𝒊 )𝟐 para la prueba de eficiencia de ajuste por raíz de error cuadrático medio
Normal Gamma Exp P. III Gum. I Log Norm Log P.III
683.85 201.63 438929.13 329.08 735.12 91.72 33.55
205.90 39.38 439943.86 80.27 628.07 9.84 0.27
0.09 31.16 424826.58 13.06 946.91 61.67 95.73
1.13 29.37 420580.25 15.59 681.33 50.13 72.62
8.57 37.63 412278.08 25.39 539.94 53.67 70.39
120.40 175.98 393931.35 155.49 750.14 198.24 220.13
219.61 267.64 379457.85 249.99 779.98 284.24 301.37
69.89 85.87 376443.76 79.68 345.86 90.23 95.64
176.33 186.10 358529.05 181.75 440.27 186.14 188.60
216.23 213.03 344166.23 212.91 405.81 207.21 204.96
314.35 296.07 327398.84 300.67 446.41 283.18 275.56
667.82 606.56 290661.63 624.58 659.96 573.62 550.79
667.82 606.56 290661.63 624.58 659.96 573.62 550.79
725.22 646.50 275210.91 670.10 636.03 606.24 577.54
361.37 297.47 268811.18 316.62 251.64 266.66 244.66
392.92 317.84 253001.83 340.49 235.54 282.53 257.01
70.07 38.06 249174.06 47.10 7.10 25.59 17.65
144.57 92.85 230055.94 108.01 26.78 71.29 56.45
89.30 47.56 216896.43 59.44 2.37 31.71 21.48
14.29 1.16 206362.75 3.77 26.18 0.06 1.71
49.39 17.61 187718.33 26.06 7.76 7.89 2.89
7.81 32.96 180345.35 22.98 181.18 51.53 69.48
93.67 161.91 170121.55 137.87 446.57 201.76 237.28
41.11 91.06 151624.59 72.77 345.75 122.31 151.19
30.12 75.54 135374.11 58.57 336.47 104.81 132.39
17.58 55.58 119257.09 40.91 312.04 81.45 106.52
7.24 36.00 103510.10 24.24 279.74 57.57 79.30
100
Tabla 4-19 Valores (𝒙𝒊 − 𝒙𝒊 )𝟐 para la prueba de eficiencia de ajuste por raíz de error cuadrático medio
Normal Gamma Exp P. III Gum. I Log Norm Log P.III
0.47 16.21 88073.47 8.62 231.85 31.64 48.44
138.00 228.32 80651.21 196.34 715.35 279.39 326.26
266.31 387.74 69577.00 345.51 1006.21 453.49 513.06
221.33 332.93 55950.92 293.85 938.13 393.77 449.53
286.62 411.43 44828.79 367.96 1087.50 478.19 539.36
398.57 542.44 34876.97 492.67 1314.24 617.76 686.66
342.78 475.23 24356.66 429.19 1224.04 544.67 608.85
250.14 362.82 15081.56 323.29 1050.60 422.43 478.38
163.35 254.26 7705.45 221.95 866.78 303.12 349.90
118.61 195.29 2740.24 167.74 758.99 236.89 277.44
11.16 38.64 757.70 27.64 391.43 56.42 75.65
11.16 38.64 757.70 27.64 391.43 56.42 75.65
0.00 7.60 4897.62 3.38 263.40 15.82 26.14
0.20 9.31 12303.65 4.75 267.64 17.53 27.71
2.40 15.84 23628.47 9.98 290.62 25.19 36.33
0.15 6.89 40633.65 3.48 236.20 12.65 20.14
5.90 0.16 64646.20 1.19 142.16 0.17 1.56
36.97 18.42 96954.56 24.08 56.30 13.12 8.31
1.97 8.53 130677.91 5.72 197.76 11.79 16.46
4.06 0.65 182591.89 1.52 91.01 0.22 0.00
53.42 41.55 250625.85 45.77 8.57 39.85 35.71
45.35 39.28 330365.25 41.72 3.72 40.53 38.51
103.22 102.98 436970.66 103.89 11.45 110.55 111.49
25.80 31.30 560739.84 29.91 0.34 39.34 43.01
43.82 60.12 732660.72 55.17 24.05 77.43 87.86
24.66 46.77 953660.60 39.34 44.96 69.35 85.57
30.18 7.39 1234611.74 12.80 0.41 0.47 0.45
10.22 0.52 1667978.42 0.26 75.98 11.93 28.84
63.73 180.18 2350261.39 136.72 781.11 292.43 391.29
22.20 8.71 3408901.73 0.24 761.48 64.51 140.66
3273.05 2077.91 5540227.95 2435.20 2.40 1439.46 1022.70
RECM:
13.99 13.54 670.74 13.57 20.49 13.60 13.88
Podemos apreciar que varias de la mayoría de las pruebas tienen resultados similares, sin embargo,
la función Gamma fue la que en general tuvo un mejor desempeño, así que esta función es
seleccionada, y los resultados pueden verse en la Tabla 4-20.
101
Tabla 4-20 Aproximaciones de caudal para diferentes tiempos de retorno utilizando la función Gamma
como modelo probabilístico de ajuste
Prob. De Ex. Tr (años) Qmax Esp.
(m3/s)
0.0001 1.0001 545.32
0.0005 1.0005 575.29
0.001 1.0010 589.68
0.005 1.0050 627.77
0.01 1.0101 646.84
0.02 1.0204 668.12
0.05 1.0526 700.92
0.1 1.1111 730.96
0.2 1.2500 768.49
0.3 1.4286 796.34
0.5 2.0000 843.85
0.7 3.3333 893.21
0.8 5 923.98
0.9 10 967.81
0.95 20 1005.02
0.96 25 1016.03
0.98 50 1048.00
0.99 100 1077.32
0.995 200 1104.61
0.999 1000 1162.31
0.9995 2000 1185.29
0.9999 10000 1235.44
Ahora bien, como se mencionó en 3.3 Generación de hidrogramas sintéticos por el método de la
SCS distribuidos a lo largo del eje del río principal, el proceso de creación de hidrogramas sintéticos
que puedan ser simulados dentro de HEC-RAS implica conocer datos de las subcuencas que le
corresponden a cada estación ficticia. Utilizando las herramientas Integrated Water Management y
las coordenadas mencionadas en la Tabla 4-5 se determinó el área tributaria de que corresponde a
102
cada estación ficticia (Ver Figura 3-13) junto con sus principales atributos, mismos que se resumen
en la Tabla 4-21.
Estación ficticia
Estación Hidrométrica
103
Tabla 4-21 Características de las subcuencas correspondientes a las estaciones ficticias
SUBCUENCA
PARÁMETRO 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Área (km2) 9.27 39.55 43.52 57.31 202.22 469.32 276.89 28.79 28.79 131.02 1247.36 129.28 31.92 40.77 266.41 84.91 16.06 40.23 102.24 978.41
Perímetro (km) 18.93 49.08 45.36 54.3 108.66 144.15 117.81 31.5 31.5 87.3 308.7 72.93 32.19 42.09 135.87 69.75 26.91 58.71 100.26 330.78
Elevación media de 2399.43 2260.39 2110.35 1878.74 1656.18 1273.47 1104.49 620.45 620.45 984.63 1337.75 740.2 751.85 394.94 586.41 273.85 96.68 34.02 14.17 450.04
la cuenca (msnm)
Pendiente media de 19.6 20.33 18.73 16.57 18.77 17.38 20.1 22.45 22.45 19.73 20.56 21.09 26.69 21.33 22.36 20.17 9.67 3.9 1.02 14.13
la cuenca (°)
Pendiente media de
38.06 39.55 36.66 31.53 35.78 32.97 38.67 46.63 46.63 37.58 58.64 41.82 54.52 42.46 44.27 40.23 19.07 7.38 1.84 28.11
la cuenca (%)
Coeficiente de 1.75 2.2 1.94 2.02 2.15 1.88 2 1.66 1.66 2.15 2.46 1.81 1.61 1.86 2.35 2.13 1.89 2.61 2.8 2.98
compacidad
Relación de 0.32 0.21 0.27 0.24 0.22 0.28 0.25 0.36 0.36 0.22 0.16 0.31 0.39 0.29 0.18 0.22 0.28 0.15 0.13 0.11
circularidad
Relación
0.95 2.2 0.84 0.83 1.49 1.21 1.55 1.57 1.57 1.52 1.35 1.62 1.36 3.46 2.17 -1 11.64 -1 0.14 -1
hipsométrica
Longitud del eje del 5.74 13.54 10.44 12.34 27.17 41.88 36.2 9.96 9.96 23.99 66.73 21.06 9.96 14.95 47.2 18.72 9.52 22.01 39.96 124.99
río principal (km)
Longitud directa del 4.71 10.07 4.31 9.27 17.31 23.32 13.07 5.8 5.8 17.67 41.76 15.69 6.54 11.74 22.11 9.61 6.42 10.94 21.86 68.92
río principal (km)
Coeficiente de
sinuosidad 1.22 1.35 2.42 1.33 1.57 1.8 2.77 1.72 1.72 1.36 1.6 1.34 1.52 1.27 2.14 1.95 1.48 2.01 1.83 1.81
hidráulico
Elevación máxima
del río principal 2608 2714 2334.87 2199.31 2397.53 1630 1598 234.25 234.25 1711 2372 1177 1264 1471 1244 791 384.54 65 -6.06 1294
(msnm)
Elevación mínima
del río principal 2134 1978 1521 1328.97 834 452.92 234.25 234.11 234.11 211.99 196.72 44.05 43.72 13.39 7.36 -6.06 5.85 -6.06 -6.48 -7.24
(msnm)
Pendiente de media
11.63 7.69 11.21 9.74 8.39 4.09 5.47 0 0 8.88 5.12 7.21 15.82 13.75 3.91 5.9 4.88 0.51 0 1.6
del río principal (%)
Tiempo de
concentración 0.66 1.51 1.07 1.27 2.52 4.63 3.7 28.59 28.59 2.22 6.26 2.12 0.87 1.3 5.22 2.12 1.29 6.49 93.21 15.75
Kirpich (h)
Tiempo de
concentración
0.67 1.51 1.08 1.28 2.53 4.66 3.72 28.75 28.75 2.23 6.3 2.14 0.87 1.3 5.25 2.13 1.3 6.53 93.71 15.84
California Highways
and public works (h)
104
Tabla 4-22 Lámina por exceso de precipitación sobre la cuenca por periodo de retorno
Caudal Base 124 M3/s
3
T (años)r Q (m /s) Q-Qbase Vol. Esc. (m3) Lámina (m) Lámina(mm)
5 923.9849 799.9849 69118697.9362 0.0158 15.8307
10 967.8091 843.8091 72905103.3783 0.0167 16.6979
20 1005.0160 881.0160 76119781.5127 0.0174 17.4342
25 1016.0288 892.0288 77071288.6668 0.0177 17.6521
50 1048.0026 924.0026 79833824.8049 0.0183 18.2848
100 1077.3192 953.3192 82366779.8069 0.0189 18.8649
Tabla 4-23 Distribución acumulada de las aportaciones de las estaciones ficticias al Qmax medio diario
esperado por periodo de retorno
Tr
No. Est. Fic. 5 10 20 25 50 100 Promedio
%Qmax acumulado
1 0.34% 0.52% 0.33% 0.41% 0.37% 0.46% 0.42%
2 0.67% 1.19% 0.65% 0.99% 0.99% 1.07% 0.98%
3 0.87% 1.71% 0.85% 1.42% 1.53% 1.48% 1.40%
4 2.65% 3.99% 2.70% 3.46% 3.62% 3.76% 3.51%
5 2.82% 4.48% 2.92% 3.86% 4.11% 3.99% 3.87%
6 4.58% 8.03% 4.95% 6.77% 8.01% 6.93% 6.94%
7 6.76% 10.43% 7.28% 8.80% 10.62% 9.19% 9.26%
8 9.09% 13.50% 9.81% 11.09% 13.70% 11.95% 12.01%
9 10.89% 16.04% 11.80% 12.87% 16.12% 14.11% 14.19%
10 11.49% 17.46% 12.59% 13.76% 17.36% 15.14% 15.26%
11 26.17% 37.47% 28.44% 28.84% 37.27% 37.89% 33.98%
12 33.99% 50.08% 37.01% 37.17% 48.39% 47.62% 44.05%
13 36.66% 53.91% 39.94% 39.63% 51.97% 49.99% 47.09%
14 43.96% 62.64% 47.76% 45.22% 60.24% 55.76% 54.32%
15 50.87% 71.24% 55.44% 52.12% 64.92% 59.44% 60.63%
16 52.95% 74.65% 57.89% 55.07% 65.29% 60.93% 62.77%
17 54.76% 77.69% 59.96% 57.30% 65.76% 62.39% 64.62%
18 55.93% 79.70% 61.32% 58.62% 66.30% 62.99% 65.79%
19 61.76% 89.34% 67.91% 66.01% 66.48% 67.51% 71.45%
20 100.00% 100.00% 100.00% 100.00% 100.00% 100.00% 100.00%
Para cada una de las estaciones ficticias se calcula el caudal que aporta al caudal total, ya que ese es
el caudal con el que serán construidos sus hidrogramas. Por cuestiones de simulación, el caudal base
105
no se tomó en cuenta (porque debe ser ingresado por separado en los datos). Por esto, los caudales
considerados fueron los que se muestran en la columna 3 de la Tabla 4-22 y su distribución se tabula
en la Tabla 4-24.
Tabla 4-24 Caudales medios diarios aportados y lámina escurrida por cada estación ficticia por periodo
de retorno (Restando el flujo base) (Pte. 1)
106
Tabla 4-25 Caudales medios diarios aportados y lámina escurrida por cada estación ficticia por periodo
de retorno (Restando el flujo base) (Pte. 2)
En algunas de las subcuencas se presentaron casos especiales en los que el río principal de la
subcuenca es el fragmento del río principal de la cuenca general (se puede ver el ejemplo de esto en
la Figura 4-30). Para estos casos la fórmula de Kirpich arroja valores desproporcionados (Pueden
verificarse en la Tabla 4-21) así que se optó por no utilizarlos y utilizar el tiempo de escurrimiento
del tramo del río principal.
Ahora, utilizando las fórmulas descritas en 2.9 Hidrograma sintético de la SCS, se estiman para cada
estación ficticia los valores de Tp y qp.
107
Tabla 4-26 Tiempo de concentración total para las estaciones ficticias
Calculado con la
No. Est. Fic. Tc local (Hrs.) T esc. (Hrs.) TcTotal
1 0.66 27.375 28.0350933 fórmula de Manning
2 1.51 27.080 28.5904618
3 1.07 26.775 27.8447025
4 1.27 26.591 27.8605968
5 2.52 25.655 28.1746864
6 4.63 23.528 28.1584386
7 3.7 22.671 26.3708182
8 0.65 22.020 22.6700066
9 2.35 20.566 22.9157204
10 2.22 19.862 22.0815264
11 6.26 18.034 24.2935835
12 2.12 16.558 18.6778335
13 0.87 16.176 17.0461618
14 1.3 15.249 16.5492359
15 5.22 13.243 18.4629309
16 2.12 11.758 13.8778086
17 1.29 11.300 12.5901766
18 6.49 7.796 14.2859654
19 6.116 1.680 7.7959238
20 15.75 0.000 15.75
108
Subcuenca
Río principal subcuenca
Río principal cuenca gen.
Figura 4-30 Subcuenca de la estación ficticia 8, con detalle del río principal de la subcuenca y de la cuenca
general
Tabla 4-27 Factores de forma para el hidrograma sintético de la SCS por estación ficticia
No. Est. Fic. Tp qp
1 18.696 1.027
2 19.029 4.323
3 18.582 4.871
4 18.592 6.412
5 18.780 22.397
6 18.770 52.007
7 17.698 32.543
8 15.477 3.869
9 15.625 22.008
10 15.124 18.019
11 16.451 157.708
12 13.082 20.555
13 12.103 5.486
14 11.805 7.184
15 12.953 42.829
16 10.202 17.312
17 9.429 3.543
18 10.447 8.010
19 6.553 32.454
20 11.325 179.697
109
Con éstos datos es posible determinar el hidrograma de las subcuencas, considerando un
escurrimiento esta hasta la desembocadura (es decir, el hidrograma producido por la subcuenca, si se
midiera en la desembocadura), en su marco de tiempo local (Un ejemplo es el tabulado en Tabla 4-28
y mostrado en la Figura 4-31) . Como paso final, se modifica el marco de referencia de tiempo al de
la desembocadura, haciendo que el escurrimiento comience con un retraso igual al tiempo de
escurrimiento desde la estación ficticia hasta la desembocadura, que puede verse en la Tabla 4-26.
110
3.00
2.50
2.00
Q (m3/s)
1.50
1.00
0.50
0.00
0.00 20.00 40.00 60.00 80.00 100.00
T (Horas)
Una vez determinados para cada estación, para realizar la convolución (dado que tienen pasos de
tiempo que no empatan) se elaboró la herramienta EXTHID (Ver 6.1.21 EXTHID), que cambia el
paso de tiempo de los hidrogramas interpolando linealmente para que puedan unirse con facilidad
(Ver Figura 4-32).
Figura 4-32 Hidrogramas por estación ficticia y su convolución para Tr=5, correspondientes a un día
de precipitación.
700
600
500
400
300
200
100
0
-100 0 1000 2000 3000 4000 5000 6000 7000 8000
Est. Fic. 1 Est. Fic. 2 Est. Fic. 3 Est. Fic. 4 Est. Fic. 5
Est. Fic. 6 Est. Fic. 7 Est. Fic. 8 Est. Fic. 9 Est. Fic. 10
Est. Fic. 11 Est. Fic. 12 Est. Fic. 13 Est. Fic. 14 Est. Fic. 15
Est. Fic. 16 Est. Fic. 17 Est. Fic. 18 Est. Fic. 19 Est. Fic. 20
Combol.
Se repite este proceso para los periodos de retorno considerados, y se obtienen los hidrogramas por
convolución que se muestran en la Figura 4-33.
111
700
600
500
Tr= 5
400
Q (m3/s)
Tr= 10
300 Tr= 20
Tr= 25
200 Tr= 50
0
0 1000 2000 3000 4000 5000 6000 7000 8000
T (mins)
Figura 4-33 Convoluciones de los hidrogramas de las estaciones ficticias por periodo de retorno
Analizando los datos de las convoluciones, determinamos los caudales pico para cada uno de los
periodos de retorno, que se enlistan en la Tabla 4-29. En el análisis de la gráfica de caudales
instantáneos se observa que en tormentas aisladas, el comportamiento de la cuenca es similar al de
los hidrogramas mostrados en la Figura 4-33, une ejemplo de esto puede verse en la Figura 4-34.
112
Hidrograma de la tormenta del 26/06/2005 al
01/07/2005
510
490
470
450
Q (M3/s)
430
410 Series1
390
370
350
0 1000 2000 3000 4000 5000 6000 7000 8000
T (mins)
4.4 Simulaciones 1D
Para esto, necesitamos obtener secciones transversales acordes a la zona de estudio. Para esto, nos
servimos de las herramientas elaboradas MAPMANNING, SECCTRASLAPE y EXPGEO (Ver
6.1.10 MAPMANNING, 6.1.6 SECCTRASLAPE, 6.1.14 EXPORTGEO )
MAPMANNING, sirve para importar un ráster de uso de suelo y cambiar las claves de uso de suelo
por el valor del coeficiente n de Manning que les corresponde, tomando éstos datos de un archivo de
texto proporcionado al programa (en una columna el valor de la clave, en la segunda el valor del valor
de n, separados por un tabulador, Figura 4-35). Se escogió este formato porque es el estándar al copiar
y pegar desde Excel) y referenciando cada coordenada matricial con su respectivo rango de
coordenadas UTM (Ver los diagramas de flujo de la Figura 3-10 a la Figura 3-12).
113
Figura 4-35 Ejemplo de archivo de texto para la conversión de la matriz de claves de uso de suelo por su
respectivo valor n
SECTRASLAPE sirve para seccionar a lo largo del eje del río proporcionando un espaciado mínimo
entre secciones y la longitud de éstas, verificando que no existan colisiones entre secciones, y
agregando datos a la matriz de secciones como el coeficiente de Manning, la posición de los puntos
de bancos y diques, y las coordenadas en formato HEC-RAS.
EXPGEO sirve para exportar los datos de las secciones transversales a un formato *.GEO como lo
requiere HEC-RAS. Para realizar la conversión del ráster de uso de suelo a n de Manning, se utilizan
los valores medios propuestos por Manning (1994) y empatados con las claves de uso de suelo usados
(los datos fueron obtenidos de CONABIO, como se mencionó en 3.2.4 Datos de uso de suelo, optando
por la capa denominada “Cobertura del suelo de México a 30 metros, 2010”, que se puede ver en la
Figura 4-36), mostrados en la Tabla 4-30.
114
Tabla 4-30 Coeficientes de rugosidad de Manning para las claves de uso de suelo
Clave Descripción n Fuente*
1 Bosque de coníferas templado o subpolar 0.1 D-2.d.4
2 Bosque de coníferas (taiga) subpolar 0.1 D-2.d.4
3 Bosque de latifoliadas perennifolio tropical o subtropical 0.1 D-2.d.4
4 Bosque de latifoliadas caducifolio tropical o subtropical 0.1 D-2.d.4
5 Bosque de latifoliadas caducifolio templado o subpolar 0.1 D-2.d.4
6 Bosque mixto 0.1 D-2.d.4
7 Matorral tropical o subtropical 0.07 D-2.c.4
8 Matorral templado o subpolar 0.07 D-2.c.4
9 Pastizal tropical o subtropical 0.035 D-2.a.2
10 Pastizal templado o subpolar 0.035 D-2.a.2
11 Matorral con líquenes y musgos subpolar o polar 0.05 D-2.c.2
12 Pastizal con líquenes y musgos subpolar o polar 0.035 D-2.a.2
13 Suelo desnudo con líquenes y musgos subpolar o polar 0.04 D-3.a
14 Humedal 0.08 C.e.1
15 Suelo agrícola 0.04 D-2.b.3
16 Suelo desnudo 0.04 D-2.c.2
17 Asentamiento humano 0.016 B-2.i.2
18 Cuerpo de agua 0.035 D-3.b
19 Nieve y hielo 0.01 A-2.b
* Sección de la tabla 5.6 de Chow (1994) de la que se tomó el valor n.
115
Figura 4-36 Distribución de uso de suelo en la zona de estudio
La función SECCTRASLAPE usa el algoritmo descrito en la Figura 4-13 al mismo tiempo en las
matrices de elevación y de uso de suelo, de forma que, para cada coordenada matricial, tenemos un
valor de elevación y uno de n de Manning. Además, hace uso de la función BANYDIQS, para
determinar los puntos de banco y dique, de forma que de los puntos que tenemos para cada sección,
podamos diferenciar cuales pertenecen al canal principal y cuales a los bancos. Esto es con el fin de
determinar el coeficiente n de Manning compuesto para cada una de estas secciones, usando cada
punto de la sección como una subdivisión, dentro de la Ecuación 4-5.
. /
∑ 𝑃𝑛
𝑛 = Ecuación 4-5
𝑃
Donde:
nc: coeficiente de Manning compuesto
N: número de subdivisiones de cada sección del río
116
Pi: Perímetro mojado para la subdivisión i de la sección
ni: Coeficiente n de Manning para la subdivisión i de la sección
P: Perímetro mojado total de la sección del río
Para la función EXPGEO, exportamos los datos calculados y obtenidos para cada sección a
un archivo de texto que tiene la estructura mostrada en la Figura 4-37. La estructura detallada de los
archivos puntos geo puede verse en el Anexo 6.2 Estructura de archivos *.geo)
El archivo *.geo puede importarse directamente a HEC-RAS mediante la opción: importar formato
GIS.
En cuanto a los caudales a transitar son los obtenidos por el método de la SCS, mostrados en la Tabla
4-29, son divididos de acuerdo a la distribución señalada en la Tabla 4-23 y se agrega a cada uno el
caudal base. Pueden verse en la Tabla 4-31 la distribución empleada para la simulación 1D.
117
Tabla 4-31 Caudales pico (en m3/s) distribuidos por estación ficticia y periodo de retorno
No. Est. Fic. Tr=5 años Tr=10 años Tr=20 años Tr=25 años Tr=50 años Tr=100 años
1 127.22 127.37 127.56 127.56 127.75 127.80
2 131.51 131.87 132.13 132.26 132.65 132.85
3 134.67 135.24 135.73 135.86 136.31 136.69
4 150.76 152.29 153.55 153.92 154.97 155.95
5 153.59 155.15 156.56 156.88 158.09 159.21
6 175.63 178.43 180.82 181.55 183.62 185.45
7 192.65 196.38 199.59 200.58 203.31 205.85
8 213.80 218.62 222.86 224.10 227.65 230.90
9 229.82 235.59 240.47 241.99 246.18 250.03
10 237.65 243.84 249.13 250.75 255.22 259.41
11 383.09 397.30 409.37 412.88 423.29 432.77
12 460.86 479.36 495.06 499.67 513.13 525.43
13 484.06 503.77 520.51 525.44 539.88 553.09
14 535.91 558.52 577.68 583.31 599.76 614.92
15 585.80 611.14 632.62 638.93 657.46 674.37
16 596.81 622.70 644.67 651.18 670.19 687.42
17 604.59 630.96 653.30 659.89 679.06 696.74
18 608.29 634.84 657.36 663.99 683.34 701.15
19 609.40 636.01 658.59 665.24 684.61 702.48
20 649.01 677.81 702.21 709.41 730.33 749.67
Los resultados obtenidos demuestran que para el escurrimiento calculado para un periodo de retorno de hasta
cien años (Ver Tabla 4-31) no existe desbordamiento importante en el canal del rio (Ver
Figura 4-38).
118
Figura 4-38 Situación de la superficie de agua para un Tr=100 años para simulación en modelo 1D
Ya que en tal caso no resulta relevante comparar la extensión de la planicie de inundación entre los
métodos 1D y 2D de simulación, determinamos las elevaciones de la superficie de agua en la
desembocadura de la cuenca (Ver Figura 4-39) a través del tiempo para poder llevar a cabo la
comparación respecto al modelo 2D.
119
Elevaciones de la superficie de agua
7
6.5
6
5.5 5
5 10
H (m)
4.5 20
4 25
3.5 50
3 100
2.5
0 100 200 300 400 500 600 700
T (Min)
Figura 4-39 Elevaciones de la superficie del agua determinadas por modelación 1D en la desembocadura
de la cuenca.
4.5 Simulaciones 2D
Para las simulaciones 2D es necesario realizar un mallado sobre la superficie de estudio (Ver
Figura 4-40). El procedimiento para el trazado de esta mala puede realizarse manualmente, por
medio de coordenadas o por una combinación de ambas.
Es recomendable usar la simulación 1D inicial como punto de referencia, para evitar los problemas
que se ilustra en la Figura 3-18 y laFigura 3-19 al dimensionar inadecuadamente la malla. De esta
forma, en las secciones donde el río no es susceptible a inundaciones, la malla puede ser más estrecha,
apegándose al cauce del río. Mientras que, en zonas con zonas de inundación, el área de la malla debe
ampliarse.
120
Figura 4-40 Mallado sobre la zona de estudio
Es necesario de igual forma importar los datos de uso de suelo, procesados de acuerdo con el
algoritmo descrito en la Figura 3-11 con los valores de la Tabla 4-30, con los datos que han sido
determinados para este propósito en 4.4 Simulaciones 1D.
Los datos son ingresados dentro del sistema 2D mediante celdas en la malla 2D, como se describe en
3.4.3 Simulaciones hidrodinámicas de avenidas en 2D usando HEC-RAS (Ver Figura 4-41) ,
ingresando los hidrogramas para cada estación ficticia y periodo de retorno que fueron calculados de
acuerdo con 4.3 Hidrogramas sintéticos por el método de la SCS (Ver Figura 4-32).
121
Caudal a través de la sección 6+688 para un Tr=5
Q (m3/s)
T (Dias)
T (Dias)
122
Figura 4-42 Situación de la superficie de agua para un Tr=100 años para simulación en modelo 2D
Los resultados obtenidos demuestran que para el escurrimiento calculado para un periodo de retorno de hasta
cien años (Ver Tabla 4-31) no existe desbordamiento importante en el canal del rio (Ver Figura 4-42).
123
Elevaciones de la superficie de agua
7
6.5
5.5 5
5 10
H (m)
4.5 20
4 25
50
3.5
100
3
2.5
0 100 200 300 400 500 600 700
T (Min)
Figura 4-43 Elevaciones de la superficie del agua determinadas por modelación 2D en la desembocadura
de la cuenca por periodo de retorno.
Comparando los flujos y las elevaciones entre las simulaciones 1D y 2D obtenemos la Figura 4-44
400
Sim. 1D
300
Sim. 2D
200
100
0
0 100 200 300 400 500 600 700
T (mins)
Figura 4-44 Caudal a través de la desembocadura de la cuenca, calculados con modelos 1D y 2D para
Tr=5 años
124
Caudales en la desembocadura para Tr=10 años
800
700
600
500
Q (m3/s)
400
Series1
300
Series2
200
100
0
0 100 200 300 400 500 600 700
T (mins)
Figura 4-45 Caudal a través de la desembocadura de la cuenca, calculados con modelos 1D y 2D para
Tr=10 años
400
Series1
300
Series2
200
100
0
0 100 200 300 400 500 600 700
T (mins)
Figura 4-46 Caudal a través de la desembocadura de la cuenca, calculados con modelos 1D y 2D para
Tr=20 años
125
Caudales en la desembocadura para Tr=25 años
800
700
600
500
Q (m3/s)
400
Series1
300
Series2
200
100
0
0 100 200 300 400 500 600 700
T (mins)
Figura 4-47 Caudal a través de la desembocadura de la cuenca, calculados con modelos 1D y 2D para
Tr=25 años
500
400 Series1
300 Series2
200
100
0
0 100 200 300 400 500 600 700
T (mins)
Figura 4-48 Caudal a través de la desembocadura de la cuenca, calculados con modelos 1D y 2D para
Tr=50 años
126
Caudales en la desembocadura para Tr=100 años
900
800
700
600
Q (m3/s)
500
400 Series1
300 Series2
200
100
0
0 100 200 300 400 500 600 700
T (mins)
Figura 4-49 Caudal a través de la desembocadura de la cuenca, calculados con modelos 1D y 2D para
Tr=100 años
Realizando un análisis por mínimos cuadrados, se determina que la diferencia en flujo calculado
entre la simulación 1D y 2D para cada periodo de retorno y se tabula en la Tabla 4-32.
.
Tabla 4-32 Diferencia de flujo promedio en la desembocadura entre simulaciones 1D y 2D
Tr (años) Error (m3/s)
5 6.08
10 6.93
20 6.70
25 6.78
50 6.99
100 7.22
127
7
5
H (msnm)
3 Elev. 1D
Elev. 2D
2
0
0 100 200 300 400 500 600 700
T (min)
Figura 4-50 Comparación de las elevaciones del espejo de agua en el exutorio para modelos 1D y 2D
Por otro lado, la diferencia en los niveles calculados entre los modelos 1D y 2D en el exutorio a través
del tiempo es de 0 m, como puede apreciarse en la Figura 4-50.
Dados los resultados obtenidos y mostrados en la Tabla 4-32 y Figura 4-50, se considera que las
simulaciones en 1D y 2D proporcionan resultados equiparables, lo que señala que la elección de datos
para los sistemas fue apropiada.
Ya que, como puede apreciarse en las figura Figura 4-38 y Figura 4-42, con los caudales considerados
dentro de los periodos de retorno en estudio no produjo desbordamientos considerables, que permitan
analizar el comportamiento de las simulaciones usando esquemas combinados 1D/2, se optó por
provocar desbordamientos en secciones determinadas del rio aumentando el caudal específicamente
para observar el comportamiento de la combinación de modelos 1D/2D.
Dados los resultados obtenidos para simulaciones 1D y 2D sin que exista un desbordamiento del
cauce principal, no existe una diferencia sensible en los resultados de las elevaciones de superficie
del agua entre uno y otro, por lo que es más conveniente utilizar una simulación 1D.
128
Una vez que el agua comienza a extenderse más allá del canal principal, debe considerarse la
complementación del modelo mediante una malla 2D. En particular, los casos en los que la dirección
del agua de la planicie de inundación se aparta considerablemente de la dirección del eje del rio.
Un análisis en 1D no es factible en zonas del río donde la planicie de inundación tiene un nivel por
debajo del nivel de desborde del río (Ver Figura 4-51). Concretamente, se debe a qué en tal caso, la
dirección del agua no se ajusta a la dirección del eje del río. En la figura Figura 4-52 se aprecia un
ejemplo de esto (en este caso, se generó el desbordamiento usando un hidrograma correspondiente a
Tr=100 años, pero aumentando el caudal base de 124 m3/s a 500 m3/s).
Elevación (msnm)
DO (m)
Figura 4-51 Ejemplo de planicie de inundación debajo del nivel de desborde del río
En tales casos, es indispensable complementar el modelo 1D con una malla 2D unida de forma lateral,
de manera que el exceso de volumen sea distribuido de manera adecuada al terreno adyacente.
129
A)
B)
Figura 4-52 Error en la simulación para terrenos con grandes planicies de inundación por debajo del
nivel de desborde del río. A) Planicie de inundación obtenida por simulación 1D B) Planicie de inundación
obtenida por simulación combinada 1D/2D.
130
Los modelos 1D también presentan problemas de modelado cuando dentro de una misma sección
tiene más de un canal, sea por la existencia de canales adyacentes al canal principal o bien, por la
existencia de meandros que causen que una sección intercepte dos veces con el eje del río.
El problema de los canales adyacentes puede resolverse moviendo los puntos de dique dentro de la
sección a lugares apropiados, lo que permite que la simulación se corresponda con la realidad. Sin
embargo, esto implica que, para cada caudal simulado, sea necesario revisar individualmente cada
sección transversal para determinar si la posición de tales puntos es la correcta (Ver Figura 4-53 A.1)
y A.2), debe determinarse y señalar manualmente que canales deben estar inundados), lo que en ríos
donde exista una gran cantidad de secciones transversales, este sea un trabajo extremadamente
laborioso, propenso a error humano y en algunas situaciones, ni siquiera exista una posición de estos
puntos que permita modelar apropiadamente estos eventos.
Es por esto, que una simulación combinada, utilizando una unión frontal se considera la opción más
pertinente en estos casos (Ver Figura 4-53 B)).
En casos donde una sola sección transversal es propensa a presentar desbordamiento (sin que esto
afecte las demás secciones, ver Figura 4-54 A) ), los modelos 1D fallan en el modelado correcto de
las consecuencias de este desbordamiento puntual, principalmente porque no puede representar el
vector de velocidad en sentido diferente al que tiene el eje del rio en ese punto. En estos casos, es
pertinente realizar una combinación 1D/2D de unión lateral que permita conocer la extensión de la
inundación provocada (Ver Figura 4-54 B))
131
A.1)
A)
A.2)
B)
Figura 4-53 Meandros y canales adyacentes A) Simulación 1D A.1) Sección con dos canales inundados
A.2) Sección con tres canales inundados B) Simulación 2D (Se utilizó en la sección 62+675 un Qp=1000
m3/s)
132
A)
A)
A.2)
B)
A.4)
Se cambia la posición de los diques para
permitir que la inundación se
represente adecuadamente
Figura 4-55 Comparación de inundaciones pequeñas con reingreso de agua al sistema simulados con
métodos 1D y 2D. A) Simulación 1D B) Simulación 2D
134
En general, de las simulaciones realizadas en 1D, 2D y 1D/2D, los tiempos de simulación observados
para el proyecto se resumen en la Tabla 4-33.
Tabla 4-33
de esta Tiempos de procesamiento
investigación, de las simulaciones
se exponen a continuación
135
Conclusiones y recomendaciones
Los datos faltantes correspondientes a los cuerpos de agua profundos encontrados en el DEM
obtenido desde INEGI fueron sustituidos satisfactoriamente por un lecho sintético por medio de
secciones hidráulicamente equivalentes a las originales utilizando el método de Corum. Se determinó
que la mejor manera para implementar el método de Corum consiste en encontrar el espejo de agua
mediante un algoritmo que combine interpolaciones vectoriales con los datos matriciales del DEM.
Las simulaciones (1D, 2D y 1D/2D) fueron realizadas utilizando el DEM obtenido desde INEGI,
complementado con el lecho sintético determinado por el método de Corum, distribuyendo los
caudales obtenidos mediante un análisis estadístico a lo largo del río apoyándonos en las simulaciones
de estaciones ficticias obtenidas mediante el modelo y software CEQUEau, y obteniendo los
hidrogramas correspondientes a tales estaciones ficticias utilizando el método dela SCS, que a su vez
fue complementado con tiempos de escurrimiento hasta la desembocadura determinados con ayuda
del software HEC-RAS (Módulo de flujo 1D estable) para determinar la velocidad de escurrimiento
en cada tramo del río. Las elevaciones sobre el nivel del mar del espejo de agua obtenido en las
simulaciones en la desembocadura fueron comparadas con los datos observados para la estación
30016 (Curvas de aforo) para calibrar el coeficiente de Manning en el lecho del río y se determinó
que las diferencias entre las elevaciones de espejo de agua esperadas y observadas fueron menores a
0.006 m, que empleamos como umbral de tolerancia por lo que se concluye que los datos empleados
(Batimétricos, uso de suelo y distribución de caudales) fueron adecuados para las simulaciones.
De acuerdo a los resultados obtenidos, la respuesta básica a que modelo de simulación es más
conveniente utilizar es que los modelos 1D son recomendables cuando el todo el flujo dentro de un
tramo tenga la misma dirección que el eje del río, los modelos 2D resultan más convenientes cuando
en un tramo se presentan canales adyacentes, flujos importantes en sentidos diferentes al del eje del
río y sistemas complejos de salida y reingreso de flujo al canal principal.
Se evaluó el desempeño de las simulaciones usando modelos 1D, 2D y la combinación de estos, para
determinar la manera más eficiente de emplearlos.
136
En el caso particular del presente estudio el uso de un esquema combinado de simulación 1D/2D de
flujo no estable probó poder realizar la simulación hasta seis veces en comparación con una
simulación únicamente en 2D.
Una aproximación de los tiempos de proceso requeridos para cada tipo de simulación en el transcurso
de esta investigación, se exponen en la Tabla 4-33.
Es importante contar con una batimetría adecuada para la simulación de avenidas, y en los lugares
donde no se cuenta con esta, la creación de una batimetría sintética que sea hidráulicamente
equivalente a la original es una actividad crucial para la simulación. Este puede probar ser un proceso
complicado y laborioso, pero en 3.2.1 Datos batimétricos, se detalla un proceso que en el presente
trabajo proporcionó resultados confiables.
137
Los resultados obtenidos verifican que un esquema combinado de simulación puede incrementar la
velocidad de cálculo por más de veinte veces, siguiendo las recomendaciones siguientes
En cada caso en el que sea necesario emplear una malla 2D, debe tenerse en consideración que en
cada iteración, todas las celdas presentan carga computacional, incluso si no hay flujo a través de
ellas (celdas secas), por lo es deseable que las mallas no se extiendan innecesariamente. La mejor
manera de determinar las zonas en las que es necesario realizar una combinación de modelos 1D y
138
2D, es realizar previamente una simulación puramente 1D, transitando los caudales correspondientes
a un periodo de retorno considerado adecuado por el proyectista para verificar las posibles áreas de
inundación.
Es importante señalar que dentro de la simulación existen diversos factores que pueden conducir a
una inestabilidad general del sistema que, al ser un proceso global, vuelve difícil la tarea de encontrar
el punto exacto en que se produce la inestabilidad. Algunos de los factores comunes con que dentro
de la investigación que conducen hacia una inestabilidad de en el modelo son los siguientes.
- Presencia de saltos hidráulicos. En general, no es deseable que entre una sección transversal
y otra exista un desnivel de más de cinco metros. Cuando se presenta este caso, es conveniente
agregar secciones intermedias para evitar desniveles tan grandes. También es aconsejable
verificar que no se hayan presentado errores en la batimetría y que el canal principal esté
correctamente limitado por puntos de bancos y en su caso, puntos de dique.
- Cambios bruscos del área mojada de la sección transversal. Esto puede ocurrir por una
variedad de factores, desde tamaño insuficiente de la sección transversal, mal
posicionamiento de ellos puntos de diques.
- Mala elección de los valores del coeficiente n de Manning.
- Cambio brusco del ancho del canal entre dos secciones. Es conveniente agregar secciones
intermedias cuando se presenten estos casos.
Como nota adicional sobre la detección de errores que conducen a la inestabilidad del sistema, puede
ser más eficiente comenzar con un modelo en blanco y agregar poco a poco los elementos de cálculo
(unas cuantas secciones a la vez, agregar una por una las estructuras laterales, probar con mallas
pequeñas e ir incrementando el tramo que abarcan) y realizar simulaciones cada vez que se agreguen
elementos hasta determinar exactamente que donde se produce la falla.
139
Referencias bibliográficas
Barkau, R.L (1992): “UNET, One-Dimensional Unsteady Flow Through a full network of
open channels, Computer Program”. St. Louis, MO. US Army Corps of
Engineers, Manual de Usuario.
Barnard, T, Kuch, A., Thompson, G., Mudaliar, S. y Philips, Brett.(2007): “Evolution of an Integrated
1D/2D Modeling Package for Urban Drainage”. Contemporary Modeling of Urban Water
Systems, Monograph 15.
Baró-Suárez, J, et al (2012): “Metodología para la valoración económica de daños potenciales
tangibles directos por inundación”. Universidad Autónoma del Estado de México, Edo. De
México, México.
Bitrán, D., et al, (2003), “Impacto Socioeconómico de los principales desastres Ocurridos en la
República Mexicana en el año 2002”. CENAPRED, Coordinación de Investigación, Riesgos
Hidrometeorológicos - Ingeniería Estructural y Geotécnia y Estudios Sociales y Económicos,
México.
Bladé, E. et al: “Desarrollo de un modelo de simulación de flujo en ríos. Convenios de colaboración
CEDEX – UPC – UdC”. Universitat Politecnica de Catalunya, Universidade da Coruña,
CEDEX
Bladé, E., (2005): “Modelación del flujo en lámina libre sobre cauces naturales. Análisis integrado
con esquemas en volúmenes finitos en una y dos dimensiones”. Universitat Politècnica De
Catalunya Departament D’enginyeria Hidràulica, Marítima I Ambiental Escola Tècnica
Superior D’enginyers De Camins Canals I Ports, Barcelona.
Bladé, E., Cea, L., Corestein, G., Escolano, E., Puertas, J., Vázquez-Cendón, M.E., Dolz, J., Coll, A.
(2014). "Iber: herramienta de simulación numérica del flujo en ríos". Revista Internacional
de Métodos Numéricos para Cálculo y Diseño en Ingeniería, Vol.30(1) pp.1-10
Borrel i Nogueras, G (2008): “Introducción informal a Matlab y Octave”. Recuperado de
http://webserver.dmt.upm.es/media/files/cursomo.pdf
Brunner, G. (2016): “HEC-RAS. River Analysis System, 2D Modeling User’s Manual Versión 5.0”.
US Army Corps of Engineers ,Manual de Usuario.
Brunner, G. (2016): “HEC-RAS. River Analysis System, Reference Manual”. US Army Corps of
Engineers , Manual de Referencias.
Brunner, G.(2016). “HEC-RAS. River Analysis System, User’s Manual Versión 5.0”. US Army Corps
of Engineers ,Manual de Usuario.
140
Brunner, G., (2014): “Combined 1D and 2D Modeling with HEC-RAS”. Manual de Usuario
Brunner, G., Warner, J., Wolfe, B., Piper, S., y Marston, L (2016): “HEC-RAS. River Analysis System,
Applications Guide”. US Army Corps of Engineers Guía de Aplicaciones.
Cannon, T. y Schipper, L. (2014): “World Disasters Report, 2014. Focus on culture and risk”.
International Federation of Red Cross and Red Crescent Societies, Lyon, Francia.
Charbonneau, R., Fortin, Jp. y Morin, G. (1977): “THE CEQUEAU MODEL: DESCRIPTION AND
EXAMPLES OF ITS USE IN PROBLEMS RELATED TO WATER RESOURCE
MANAGEMENT/Le modéle CEQUEAU: description et exemples d'utilisation dans le cadre de
problémes reliés I'aménagement”, Hydrological Sciences Bulletin, 22:1, 193-202, DOI.
Coleman, J. (2006): “Holes in te ocean: Filling data voids in bathymetric lidar, a case study in dry
tortugas national park, Florida”. Tesis de Maestría en Ciencias, B.S., The Univeristy of
Georgia, 2006
Coll, A., Ribó, R., Pasenau, M., Escolano, E., Perez, J.Suit., Melendo, A. y Monros, A. (2014): “GiD
v. 12 Reference Manual”. CINME.
Coll, A., Ribó, R., Pasenau, M., Escolano, E., Perez, J.Suit., Melendo, A. y Monros, A. (2014): “GiD
v. 12. User Manual”. CINME
Cook, A. (2008): “Comparison of one-dimensional HEC-RAS with two-dimensional FESWMS model
in flood inundation mapping”. Purdue University, West Lafayette, Indiana
Corum, Zachary P. (2015): “synthetic bathymetry method development, validation and application to
five pacific northwest rivers”. Joint Federal Interagency Conference, Abril, 2015.
Diaz, V. (2010): “Diseño Geomático del modelo Hidrológico CEQUEau para cuencas no
controladas”. (Tesis de Maestría en Ciencias del Agua), Centro Interamericano de
Recursos del Agua, Toluca, México.
Díaz-Delgado, C., Gómez-Albores, M., Hernández-Perez, J.(2012): “Modelado Hidrológico-
Hidráulico de Inundaciones con Estimación de Daños Directos Tangibles”. XXII Congreso
Nacional de Hidráulica, Acapulco, Gro., México.
Díaz-Delgado,C., Gaytán, J. (2014): “Flood risk assessment in humanitarian logistics process
desing”. Centro Interamericano de Recursos del Agua (CIRA), F.I. Universidad Autónoma
del Estado de México, Edo. de México, México, Vol. 12, Octubre 2014, pp. 976-984.
Fenoglio, E. et al (2015): “Inundaciones Urbanas y Cambio Climático. Recomendaciones para la
Gestión”. Secretaría de Ambiente y Desarrollo Sustentable, Buenos Aires, Argentina.
García, N., Mendez, K. y Sánchez, S. (2015): “Impacto Socioeconómico de los Desastres en México
durante el año 2014 (Resumen ejecutivo)”. CENAPRED, D.F., México.
141
Gibbs, M., Clarke, K. y Taylor, B. (2015): “Linking spatial inundation indicators and hydrological
modelling to improve assessment of inundation extent”. Ecological Indicators.
Guerra-Cobian, V. (2007). “Análisis del efecto de discretización espacial en el modelado de cuencas
hdrológicas utilizando el modelo distribuido CECUEAU-ONU”. Tesis de doctorado. Centro
Interamericano de Recursos del Agua de la UAEMéx. 231 p.
Guerra-Cobian, V., Bâ, K., Quentin, E., (2013): “Efecto de la discretización Espacial sobre las
Simulaciones de Caudal con el Modelo distribuido CEQUEAU”. Tecnología y Ciencias del
Agua, Vol. IV, núm. 5, pp. 33-53.
Hartnack, J., Madsen, H. y Tornfeldt, J. (2005): “Data Assimilation in a combined 1D-2D Flood
Model”. Proceedings of the International Conference “Innovation, Advances and
Implementation of Flood Forecasting Technology”, Norway, 17-19 Octubre 2005.
Karim, F., Petheram, C., Marvanek, S., Ticehurst, C., Wallace, J. y Gouweleeuw, B. (2001): “The use
of hydrodynamic modelling and remote sensing to estimate floodplain inundation and flood
discharge in a large tropical catchment”. 19th International Congress on Modelling and
Simulation, Perth, Australia, 12–16 December 2011.
Kuiry, S., Sen, D y Bates, P. (2010): “Coupled 1D-Quasi-2D Flood Inundation Model with
Unstructured Grids”. Journal of Hydraulic Engineering, Vol. 136, No. 8, pp. 493-506
Lee, J (1991): “Comparison of existing methos for building triangular irregular network,
models of terrain from grid digital elevation models”. International Journal of
Geographical Information Systems, 5:3, 267-285.
Llamas, J.(1993): “Hidrología General”. Universidad del País Vasco
Lluén, W. (2015): “Aplicacion de la nueva herramienta HEC-RAS 5.0 para cálculos bidimensionales
del flujo de agua en ríos “. Escola de Camins, ,UPC BARCELONATECH, Tesis de Máster
en Ingeniería Civil
Morales-Hernández, M., Petaccia, G., García-Navarro, P. (2014): “Conservative 1D-2D coupled
numerical strategies applied to river flooding”. Informatics, Networking and Intelligent
Computing: Proceedings of the 2014 International Conference on Informatics, Networking
and Intelligent Computing, November 2014, Shenzhen, China, pp. 605-613
Morales-Hernández, M., Petaccia, G., García-Navarro, P. (2015): “Conservative 1D-2D coupled
numerical strategies applied to river flooding: The Tiber (Rome)”. Applied Mathematical
Modelling, No. 40, pp. 2087-2105.
Morin, G., Paquet, P. (2007): “Modèle hydrologique CEQUEAU”. INRS-ETE, rapport de recherche
no R000926, 458p.
142
Olcina, J. (2008): “Prevención de Riesgos: Cambio Climático,Sequías e Inundaciones”. (Panel
científico-técnico de seguimiento de la política del agua ). Zaragoza : Fundación Nueva
Cultura del Agua.
Overton, I. (2005): “Modelling floodplain inundation on a regulated river: integrating gis, remote
sensing and hydrological models”. River Research and Applications No. 21, pp. 991-1001.
Paredes, P., Pedrozo, A., y Mejía, P., (2014): “Evaluación de Modelos Numéricos 1D y 2D para
Predecir Inundaciones”. XXIII Congreso Nacional de Hidráulica, Jalisco, México.
Powell, S., Jakeman, A. y Croke., B (2014): “Can NDVI response indicate the effective flood extent
in macrophytedominated floodplain wetlands?”. Ecological Indicators.
Salas, M. A. y Jimenez, M. (2004): “Inundaciones”. CENAPRED, Serie Fascículos, 1ra
edición, Octubre 2004.
Schenk, A y Shan, J (2002): “Fusion of LiDAR data Aerial imagery, for mor complete
Surface descriptions. International Archives of Photogrametry, Remote Sensing
and Spatial information sciences
Scleiss et al (2013): “A conservative strategy to couple 1D and 2D models for shallow wáter flow
simulation”. Computers and fluids, Vol. 81, 26-44
Shaikh, M., Green, D. y Cross, H. (2001): “A remote sensing approach to determine environmental
flows for wetlands of the Lower Darling River, New South Wales, Australia”. International
Journal of Remote Sensing, Vol. 22, No. 9, pp. 1737-1751, Nov. 2010.
Shwanghart, W. y Kuhn, N. (2010): “TopoToolBox: A set of Matlab functions for
topographic analisis”. Enviromenmental Modelling & software, Vol. 25, pp. 770-
781
Thomas, R., Kingsford, R., Lu, Y., y Hunter, S. (2010): “Landsat mapping of annual inundation
(1979–2006) of the Macquarie Marshes in semi-arid Australia”. International Journal of
Remote Sensing, Vol. 32, No. 16, pp. 4545-4569, Agosto 2011.
Timbe, L. (2011): “Desempeño de modelos hidráulicos 1D y 2D para la simulación de
inundaciones”. Revista MASKANA, Vol.2, Núm. 1, pp. 91-98.
van Dijk, E., van der Meulen, J., Kluck, J. y Straatman, J (2014): “Comparing modelling techniques
for analysing urban pluvial flooding”. Water Science and Technology, Vol., 69, Num. 2, pp.
305-311
Ven Te, C. (1994): “Hidráulica de canales abiertos”. Santa Fe de Bogotá, Colombia:
McGraw-Hill
Ven Te, C. (1994): “Hidrología Aplicada”. Santa Fe de Bogotá, Colombia: McGraw-Hill.
143
Wang, R. y Ateljevich, E. (2012): “Methodology for Flow and Salinity Estimates in the Sacramento-
San Joaquin Delta and Suisun Marsh”. California Departmen of Water Resources, 33rd
Annual Progress Report, Junio, 2012.
Werner, M., Blazkova, S. y Petr, J. (2005): “Spatially distributed observations in constraining
inundation modelling uncertainties”. Hidrológical Processes, Vol. 19, pp. 3081-
3096, Mayo del 2005.
Cano, A., (2005): “Solución Informática para el análisis de frecuencia de eventos máximos
hidrológicos”. Tesis de licenciatura en ingeniería en computación, Estado de México,
Universidad Autónoma del Estado de México, 2005.
144
Anexos
6.1 Códigos fuente de las herramientas desarrolladas
6.1.1 IMPTIF
function dem = imptif()
%Esta función importa un archivo tiff desde la interfaz de usuario y
%devuelve la matriz de datos del archivo, los datos de referencia y la
%información de la imagen, que se utilizaran después para poder exportar
la
%matriz de vuelta a una imagen TIFF
[archivo directorio]=uigetfile('*.tif');
6.1.2 CALRIO
function [rio dem]=calrio(dem)
%Esta función genera una estructura con varios valores concernientes al
rio principal, utilizando la herramienta "topotoolbox"
for a=2:dem.ref.RasterSize(2)
Y(:,a)=Y(:,1);
end
DEM = GRIDobj(X,Y,dem.mat);
DEM.name='DEM';
DEM.zunit='m';
DEM.xyunit='m';
FD = FLOWobj(DEM,'preprocess','carve');
S = STREAMobj(FD,'minarea',100);
D=flowdistance(FD,'upstream');
MaxDisFlo=find(D.Z==max(max(D.Z)));
[IX,distance]=flowpathextract(FD,MaxDisFlo(1));
RP=zeros(dem.ref.RasterSize(1),dem.ref.RasterSize(2));
RP(IX)=distance;
%distancia de flujo
rio.dflujo=RP;
%RioPrincipal
rioprin=RP;
rioprin(rioprin>0)=1;
rio.rioprin=rioprin;
for cont=1:length(a)
tic
[Y X]=find(rio.dflujo==a(cont));
DIRIO(cont,1)=Y;
DIRIO(cont,2)=X;
DIRIO(cont,3)=cont;
toc*(length(a)-cont)/60
end
rio.riocal=DIRIO;
%Determinamos la dirección de cada pixel del rio principal
for cont=1:(size(DIRIO,1)-1)
146
DIRIO(cont,4)=direccion(DIRIO,cont);
end
%DIRIO=[DIRIO,direcc];
rio.riocal=DIRIO;
%para el ultimo pixel, dejamos la dirección del rio tal como venía del
%pixel anterior (porque no tenemos un pixel "siguiente" para determinar
la
%dirección
rio.riocal(length(rio.riocal),4)=rio.riocal(length(rio.riocal)-1,4);
%Colocamos las secciones del rio de acuerdo a la dirección de cada pixel
rio.seccrio=seccionar(dem,rio,20);
rmpath(genpath('topotoolbox-master'))
rmpath(genpath('Mapping toolbox2'))
6.1.3 EXP2TIF
function g=exp2tif(A,R,info)
%Esta función exporta una imagen TIFF previamente cargada con impttif
[archivo ruta] = uiputfile('*.tif');
A=flip(A);
geotiffwrite([ruta
archivo],A,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag,'Ti
ffTags',struct('Compression',Tiff.Compression.None,'PhotometricInterpreta
tion',Tiff.Photometric.MinIsBlack,'RowsPerStrip',1));
6.1.4 CLOSEGAP
demcerrado=flip(dem.mat);
[filas,cols]=size(demcerrado);
dem = cerrarsecciones(dem,rio,valnulo);
secerrado=dem.secerrado;
demcerrado(secerrado(cont2,6),secerrado(cont2,5))=secerrado(cont2,4);
end
147
end
[y x]=find(demcerrado==valnulo);
while length(y)>9
for cont=1:length(y)
if y(cont)-1>0 & x(cont)-1>0 & y(cont)+1<=filas & x(cont)+1<=cols
sub=demcerrado([y(cont)-1:y(cont)+1],[x(cont)-1:x(cont)+1]);
sub(sub==valnulo)=NaN;
med=nanmean(nanmean(sub));
if med~=NaN
demcerrado(y(cont),x(cont))=med;
end
end
end
[y x]=find(demcerrado==valnulo);
length(y)
end
dem.demcerrado=demcerrado;
dem.secerrado=secerrado;
dem2=dem;
dem2.mat=dem.demcerrado;
y=find(secerrado(:,5)==riocal(cont,2)&secerrado(:,6)==riocal(cont,1)&sece
rrado(:,1)==cont);
riocal(cont,5)=secerrado(y,2);
riocal(cont,6)=secerrado(y,3);
%para el resto de las secciones
for cont=(largo-1):-1:1
y=find(secerrado(:,5)==riocal(cont,2)&secerrado(:,6)==riocal(cont,1)&sece
rrado(:,1)==cont);
riocal(cont,5)=secerrado(y,2);
riocal(cont,6)=secerrado(y,3);
riocal(cont,7)=((riocal(cont,5)-riocal(cont+1,5))^2+(riocal(cont,6)-
riocal(cont+1,6))^2)^.5;
end
148
rio.riocal=riocal;
%Colocamos las elevaciones de las secciones cerradas dentro de la matriz
%del rio calibrado
for cont=1:length(rio.riocal)
yelev=find(dem.secerrado(:,5)==rio.riocal(cont,2) &
dem.secerrado(:,6)==rio.riocal(cont,1));
rio.riocal(cont,8)=dem.secerrado(min(yelev),4);
end
%Entonces el formato de la matriz quedará:
%Ymatlab,Xmatlab,No.Sección,Dirección del pixel,Xreal,Yreal,Z,
%Finalmente, obtenemos las secciones transversales del rio una vez que ha
%sido cerrado
toc
6.1.5 CERRARSECCIONES
cont=1;
seccs=rio.seccrio;
seccs2=rio.seccrio;
while cont<=length(seccs);
%Buscamos valores nulos
while seccs(cont,4)~=valnulo
if cont+1<=length(seccs)
cont=cont+1;
else
break
end
end
if cont+1>length(seccs)
break
end
%Determinamos la sección a la que pertenece el valor nulo
seccion=seccs(cont,1);
y=find(seccs(:,1)==seccion);
cont4=1;
cont5=1;
inicio=[];
fin=[];
149
y2=length(find(seccs(y,4)==valnulo));
for cont2=1:length(y)
if y(cont2)+1<=length(seccs)
if seccs(y(cont2),4)~=valnulo & seccs(y(cont2)+1,4)==valnulo
& seccs(y(cont2),1)==seccion & seccs(y(cont2)+1,1)==seccion
inicio(cont4)=seccs(y(cont2),4);
cont4=cont4+1;
end
try
if isempty(inicio)==0 || isempty(fin)==0
if isempty(inicio)==1
inicio=99999;
end
if isempty(fin)==1
fin=99999;
end
elevesp=min(min(inicio),min(fin));
for cont2=1:length(y)
if seccs(y(cont2),4)==valnulo
seccs2(y(cont2),4)=elevesp;
end
end
end
catch
%cont+cont2
y
seccion
inicio
fin
end
150
cont=cont+y2;
end
y=seccs2(:,6);
x=seccs2(:,5);
seccbool(seccbool(:,4)~=valnulo,4)=0;
seccbool(seccbool(:,4)==valnulo,4)=1;
%creamos una matriz que solo contenga valores anómalos en las elevaciones
seccanom=secerrado;
seccanom(:,4)=seccbool(:,4).*secerrado(:,4);
seccsvalnul(contbis,2)=max(unique(seccanom(find(seccanom(:,1)==seccsvalnu
l(contbis)),4)));
end
seccsvalnul(contbis+cont3bis,2)=seccsvalnul(contbis,2)+grad*cont3bis;
end
end
end
6.1.6 SECCTRASLAPE
function seccs=secctraslape(dem,rio,tam,esp,usosuelo)
%Esta función se encarga de crear una matriz con los datos de las
secciones
%transversales a intervalos cerrados, evitando que existan traslapes
entre
%ellas, espaciándolas de forma que esto no ocurra
%Dem: Matriz numérica del Dem. Obtenida al importar Rio: Matriz con los
%datos calibrados del rio principal
%
seccs.secs=dem.secerrado(find(dem.secerrado(:,1)==rio.riocalint(1,3)),:);
%En esta parte de aquí, irá el bucle que lo haga para todas las
secciones,
%para fines de prueba, lo haré para un solo par de puntos.
cont=1;
while (cont+1)<=size(rio.riocalint,1)
p2=[rio.riocalint(cont+1,5),rio.riocalint(cont+1,6),rio.riocalint(cont+1,
3)];
while p2(3)>rio.riocalint(secsig,3)
if secsig+1<=size(rio.riocalint,1)
secsig=secsig+1;
else
secsig=secsig+1;
break
end
end
else
p2=[rio.riocal(p2(3)+cont2,5),rio.riocal(p2(3)+cont2,6),rio.riocal(p2(3)+
cont2,3)];
cont2=cont2+1;
end
153
end
seccs.secs=pruebasec;
cont=secsig;
end
seccs=seccs.secs;
%-------------------------------------------------------------------
cmanning=mapmanning(usosuelo);
for cont=1:size(seccs,1)
seccs(cont,7)=cmanning(seccs(cont,6),seccs(cont,5));
end
for cont=1:size(listasec,1)
%Vamos trabajando sección por sección
y=find(seccs(:,1)==listasec(cont));
%Encontramos en que punto se encuentran los bancos
bi=find(seccs(:,2)==bancos.biz(cont,2) &
seccs(:,3)==bancos.biz(cont,3));
bd=find(seccs(:,2)==bancos.bder(cont,2) &
seccs(:,3)==bancos.bder(cont,3));
%Colocamos las claves a cada parte de la sección: 1 banco izq 2
%canal principal 3 banco derecho
if bi~=y(1)
for cont2=y(1):bi
seccs(cont2,8)=1;
end
else
seccs(bi,8)=2;
end
for cont2=bi+1:bd-1
seccs(cont2,8)=2;
end
if bd~=y(end)
for cont2=bd:y(end)
seccs(cont2,8)=3;
end
else
seccs(bd,8)=2;
end
154
end
%Ahora se determina para cada sección el coeficiente de rugosidad
%compuesta, con la fórmula 6.17 de Chow Hidráulica de Canales Abiertos
for cont=1:size(listasec,1)
%Vamos trabajando sección por sección
y=find(seccs(:,1)==listasec(cont));
%Encontramos los valores correspondientes a cada sección del rio
yi=find(seccs(y,8)==1);
ycp=find(seccs(y,8)==2);
yd=find(seccs(y,8)==3);
%Iniciamos/Reinicializamos el Perímetro Mojado Total
pmt=0;
sum=0;
%Determinamos los valores de los coeficientes de Manning compuestos
for cont2=1:size(yi,1)
%Perímetro Mojado Parcial
pmp=seccs(yi(cont2),9);
%Coeficiente de rgosidad de cada sección
np=seccs(yi(cont2),7);
%La suma de los coeficeinte ajustadas está dada por
sum=sum+pmp*(np^(3/2));
%Perímetro mojado total
pmt=pmt+pmp;
end
ncbiz=(sum/pmt)^(2/3);
155
ncbder=(sum/pmt)^(2/3);
seccs(y(1),10)=ncbiz;
seccs(y(1),11)=nccp;
seccs(y(1),12)=ncbder;
end
6.1.7 SECINT
dist=unique(rio.dflujo);
dist=(max(max(dist))-dist);
dist=sort(dist);
dist(1)=[];
nsec=floor((dist(end)-dist(1))/int);
%Ingresamos el primer valor del intervalo
secinte(1)=0;
for cont=1:nsec;
secinte(cont+1,1)=int*cont;
end
for cont=1:nsec
a=abs(dist-secinte(cont,1));
secinte(cont,2)= find(a==min(min(a)));
end
rio.dflujoint=secinte;
%Ahora, las secciones a los intervalos fijos las colocamos dentro de una
%matriz independiente dentro del objeto del dem
y=find(dem.secerrado(:,1)==secinte(1,2));
secerradoint=dem.secerrado(y,:);
for cont=2:nsec
y=find(dem.secerrado(:,1)==secinte(cont,2));
secerradoint=[secerradoint;dem.secerrado(y,:)];
end
dem.secerradoint=secerradoint;
156
%principal, donde quitamos los pixeles innecesarios
y=find(rio.riocal(:,3)==secinte(1,2));
riocalint=rio.riocal(y,:);
for cont=2:nsec
y=find(rio.riocal(:,3)==secinte(cont,2));
riocalint=[riocalint;rio.riocal(y,:)];
end
%Modificamos las longitudes de los tramos en cada pixel del rio
for cont=1:length(riocalint)-1
riocalint(cont,7)=int;
end
rio.riocalint=riocalint;
6.1.8 SECXPUNTOS
function sec=secxpuntos(dem,p1,p2,tam)
%V1.2. Se arregla la matriz para que sea compatible con las secciones que
%genera la función "seccionar" y por lo tanto, compatible con el resto
%de las funciones creadas.
%Esta función obtiene en el punto 2 (p2) una sección transversal
%perpendicular a la dirección p1,p2.
%d=[(p1(1)-p2(1))/abs(p1(1)-p2(1)),(p1(2)-p2(2))/abs(p1(2)-p2(2))]
%El vector unitario de dirección entre las secciones está definido por:
d=[(p1(1)-p2(1)),(p1(2)-p2(2))]/(((p1(1)-p2(1))^2+(p1(2)-p2(2))^2)^.5);
d(isnan(d))=0;
dem.mat=flip(dem.mat);
%por lo que el vector unitario de dirección
di=[d(2),-d(1)];
%determinamos donde se necesitan obtener puntos sobre la sección
%esquina inferior izquierda
Xref=dem.ref.XWorldLimits(1);
Yref=dem.ref.YWorldLimits(1);
%Resoluciones en m
ResX=dem.ref.CellExtentInWorldX;
ResY=dem.ref.CellExtentInWorldY;
157
ncolum=size(dem.mat,2);
%esquina superior izquierda (Para las referencias en formato Matlab)
Xml=Xref;
Yml=Yref+ResY*nfilas;
158
end
cont=cont+1;
% % sec.coor(1,1)=coop1.utm(1);
% % sec.coor(1,2)=coop1.utm(2);
% % sec.coor(1,3)=coop1.ml(1);
% % sec.coor(1,4)=coop2.ml(2);
% % sec.coor(1,5)=dem.mat(sec.coor(1,4),sec.coor(1,3))
% sec.coor(cont,1)=coop2.utm(1);
% sec.coor(cont,2)=coop2.utm(2);
% sec.coor(cont,3)=coop2.ml(1);
% sec.coor(cont,4)=coop2.ml(2);
% sec.coor(cont,5)=dem.mat(sec.coor(cont,4),sec.coor(cont,3))
a=[0,0,0,0,0];
sec.coor=[a;sec.coor];
%Agregamos los valores de las "puntas"
sec.coor(cont,1)=coop1.utm(1);
sec.coor(cont,2)=coop1.utm(2);
sec.coor(cont,3)=coop1.ml(1);
sec.coor(cont,4)=coop1.ml(2);
sec.coor(cont,5)=dem.mat(sec.coor(cont,4),sec.coor(cont,3));
sec.coor(1,1)=coop2.utm(1);
sec.coor(1,2)=coop2.utm(2);
sec.coor(1,3)=coop2.ml(1);
sec.coor(1,4)=coop2.ml(2);
sec.coor(1,5)=dem.mat(sec.coor(1,4),sec.coor(1,3));
sec2=sec;
sec.coor(:,1)=p2(3);
sec.coor(:,2)=sec2.coor(:,1);
sec.coor(:,3)=sec2.coor(:,2);
sec.coor(:,4)=sec2.coor(:,5);
sec.coor(:,5)=sec2.coor(:,3);
sec.coor(:,6)=sec2.coor(:,4);
6.1.9 VECUNI
function vec=vecuni(p1,p2)
d=[(p1(1)-p2(1)),(p1(2)-p2(2))]/(((p1(1)-p2(1))^2+(p1(2)-p2(2))^2)^.5);
%Consideramos la situación en que hay vecotres horizontales o verticales,
%por lo que las pendientes son infinias
d(isnan(d))=0;
vec=[d(2),-d(1)];
6.1.10 MAPMANNING
function manning=mapmanning(usosuelo)
159
%de Manning. Se necesita dar la ubicación de un archivo txt que tenga en
la
%primera columna la clave de uso de suelo y en la segunda el coeficiente
de
%Manning que le corresponda, separados por tabulador (de esta manera
pueden
%copiarse esaas columnas desde excel
for cont=1:size(manning,1)
tic;
manning(cont,2)
usosuelo.mat(find(usosuelo.mat==manning(cont,1)))=manning(cont,2);
toc;
tminus='Tiempo restante mapa de C manning: '
tminus= (size(manning,1)-cont)*toc/60
end
manning=flip(usosuelo.mat);
6.1.11 TIRANTES
function [A]=erosionarsteady()
%Esta función ejecuta una simulación en el programa HEC RAS y devuelve
las elevaciones del fondo, de la superficie del agua y el tirante.
h=actxserver('RAS500.HECRASCONTROLLER');
[archivo directorio]=uigetfile('*.prj');
ras_file=[directorio '\' archivo];
h.Project_Open(ras_file);
h.ShowRas;
h.Compute_ShowComputationWindow;
h.CurrentPlanFile;
h.Compute_CurrentPlan(0,0);
h.Compute_Complete();
z9 = 'error message';
ns1=h.Schematic_XSCount
V=zeros(ns1,2);
for cont=1:ns1
ws=h.Output_NodeOutput(1, 1, cont, 0, 1, 2);
fondocanal=h.Output_NodeOutput(1, 1, cont, 0, 1, 5);
tirante=ws-fondocanal;
V(cont,1)=cont;
V(cont,2)=fondocanal;
V(cont,3)=ws;
V(cont,4)=tirante;
end
A=V;
h.Project_Save; %Saves the project
h.Compute_HideComputationWindow;
h.QuitRas;
delete(h);
160
6.1.12 DIRECCIÓN
Ydespues=dirrio(find(dirrio(:,3)==secc+1),1);
Xdespues=dirrio(find(dirrio(:,3)==secc+1),2);
6.1.13 ESPBOOL
for cont=1:size(seccsespejo.seccsnulas,1)
tic
y=find(seccsespejo.seccsnulas(:,1)==cont);
162
%Determinamos cual es el espejo principal
%Para esto, suponemos como principal el espejo que está unido
al rio
%principal
yrp=rio.riocal(cont,1);
xrp=rio.riocal(cont,2);
%Encontramos el punto dentro de las secciones que corresponde
a las
%coordenadas del rio principal para la sección
yrpse=find(seccsespejo.seccsnulas(y,5)==xrp &
seccsespejo.seccsnulas(y,6)==yrp)+y(1)-1;
if seccsespejo.seccsnulas(yrpse+cont2,4)==0
yfin=yrpse+cont2-1;
break
end
cont2=cont2+1;
end
%Buscamos hacia arriba
cont2=0;
while yrpse+cont2>=y(1)
if seccsespejo.seccsnulas(yrpse+cont2,4)==0
yini=yrpse+cont2+1;
break
end
cont2=cont2-1;
end
else
%Estamos en el caso en que el eje del rio no se encuentre
%en un espejo de agua
163
end
if yrpse+cont2==y(1) &
seccsespejo.seccsnulas(yrpse+cont3,4)==0
cont3=9999999;
end
end
if abs(cont3)<cont2
yfin=yrpse+cont3+1;
cont2=cont3;
while yrpse+cont2>=y(1)
if seccsespejo.seccsnulas(yrpse+cont2,4)==0
yini=yrpse+cont2+1;
break
end
cont2=cont2-1;
end
else
yini=yrpse+cont2-1;
cont2=cont2;
while yrpse+cont2<=y(end)
if seccsespejo.seccsnulas(yrpse+cont2,4)==0
yfin=yrpse+cont2-1;
break
end
cont2=cont2+1;
end
end
end
%Columna 7, se usará para determinar el error
yep=[yini:yfin];
seccsespejo.seccsnulas(yep,7)=seccsespejo.seccsnulas(yep,4);
%Buscamos el número de intersecciones que tiene la sección con
%el eje principal
%Reseteamos el contador de colisiones
col=0;
for cont2=1:size(y,1)
ybus=find(seccsespejo.seccsnulas(y(cont2),5)==rio.riocal(:,2) &
seccsespejo.seccsnulas(y(cont2),6)==rio.riocal(:,1));
if isempty(ybus)==0
col=col+1;
end
end
if col>1
seccsespejo.seccsnulas(y,8)=seccsespejo.seccsnulas(y,7);
else
seccsespejo.seccsnulas(y,8)=seccsespejo.seccsnulas(y,4);
end
else
%Si solo existe un espejo, copiamos la información de la
sección
%nula a las columnas 7 y 8
seccsespejo.seccsnulas(y,7)=seccsespejo.seccsnulas(y,4);
seccsespejo.seccsnulas(y,8)=seccsespejo.seccsnulas(y,4);
164
end
%Si no hay espejos, saltamos a la siguiente sección sin hacer nada
end
toc*(size(seccsespejo.seccsnulas,1)-cont)/60
end
seccsnulas=seccsespejo;
6.1.14 EXPORTGEO
function exportgeo5(dem,rio,seccs,usosuelo)
%V2.En esta versión agrego soporte para los tramos en intervalos
regulares
%V3.En esta versión agrego las posiciones de los diques y la detección de
%V4. Se cambia el nombre a exportgeo2 para mantener intacta la función
%original y se hacen los cambios requeridos para que sea compatible con
las
%secciones espaciadas para evitar traslape que se utilizaron en la
función
%secctraslape
%V.5. Se agrega la compatibilidad para el coeficiente de Manning
rio.dflujoint(cont,2)=rio.dflujo(rio.riocalint(cont,1),rio.riocalint(cont
,2));
end
tem=rio.dflujoint;
rio.dflujoint(:,1)=tem(:,2);
rio.dflujoint(:,2)=tem(:,1);
%assignin('base', 'dflujo', rio.dflujoint)
%------------------------
dist=rio.dflujoint(:,1);
165
%---------------------
nsecc=rio.riocalint(:,3)';
A=[x;y;z;nsecc];
fprintf(fileID,formato,A);
fprintf(fileID,'%0s\r\n','');
%Datos del tramo
fprintf(fileID,'%1s\r\n','REACH:');
fprintf(fileID,'%0s\r\n','');
fprintf(fileID,'%1s %1s\r\n','STREAM ID:',nr);
fprintf(fileID,'%1s %1s\r\n','REACH ID:',nomt);
fprintf(fileID,'%1s %1d\r\n','FROM POINT',1);
fprintf(fileID,'%1s %1d\r\n','TO POINT:',length(rio.riocalint));
fprintf(fileID,'%0s\r\n','');
%Colocamos los datos de la línea central del tramo
fprintf(fileID,'%1s\r\n','CENTERLINE:');
%colocamos las distancias al origen de las secciones
dist=rio.dflujoint(:,1);
166
%ahora, como la distancia debe estar medida a partir de aguas abajo, lo
%invertimos
dist=max(max(dist))-dist;
dist=dist';
%dist(1)=[];
x2=rio.riocal(:,5)';
y2=rio.riocal(:,6)';
z2=rio.riocal(:,8)';
x2(end)=[];
y2(end)=[];
z2(end)=[];
dist2=unique(rio.dflujo);
dist2(1)=[];
dist2=dist2-dist2(1);
dist2=max(max(dist2))-dist2;
dist2=dist2';
dist2(end)=[];
assignin('base', 'dflujo',dist2);
assignin('base', 'x',x2);
assignin('base', 'y',y2);
assignin('base', 'z',z2);
A=[x2;y2;z2;dist2];
assignin('base', 'A',A);
formato = '%1.2f, %1.2f, %1.6f, %1.2f\r\n';
fprintf(fileID,formato,A);
fprintf(fileID,'%1s\r\n','END:');
fprintf(fileID,'%0s\r\n','');
fprintf(fileID,'%1s\r\n','END STREAM NETWORK:');
fprintf(fileID,'%0s\r\n','');
%COLOCAMOS LOS DATOS DE LAS SECCIONES TRANSVERSALES
fprintf(fileID,'%1s\r\n','BEGIN CROSS-SECTIONS:');
fprintf(fileID,'%0s\r\n','');
formato = '%1.9f, %1.8f, %1.14f\r\n';
for cont=1:length(rio.dflujoint)
y=find(dem.secerradoint(:,1)==rio.dflujoint(cont,2));
fprintf(fileID,'%0s\r\n','');
fprintf(fileID,'%1s\r\n','CROSS-SECTION:');
fprintf(fileID,'%0s\r\n','');
fprintf(fileID,'%1s %1s\r\n','STREAM ID:',nr);
fprintf(fileID,'%1s %1s\r\n','REACH ID:',nomt);
fprintf(fileID,'%1s %1.2f\r\n','STATION:',dist(cont));
fprintf(fileID,'%1s\r\n','NODE NAME:');
%Colocamos los bancos
distsec=[bancos.biz(cont,5),bancos.bder(cont,5)];
%rio.riocalint(cont,7),rio.riocalint(cont,7),rio.riocalint(cont,7));
end
%Inicia coeficiente de rugosidad de Manning-----------
167
%Un espacio para separar de la sección anterior
fprintf(fileID,'%0s\r\n','');
%Se imprime el encabezado de la sección
fprintf(fileID,'%1s\r\n','NVALUES:');
%Continuamos en la siguiente línea
fprintf(fileID,'%0s\r\n','');
%Listamos los valores de los coeficientes compuestos por sección
%Las columnas 10, 11 y 12, en la primer fila de la sección, de la
tabla
%'secciones' contienen los valores del coeficiente de Manning
compuesto
%para la sección
ycoefm=find(seccs(:,11)~=0);
%Creamos una tabla que contenga: No.Secc,Coef.Man.Banco.Iz,
%Coef.Canal Prin y Coef. Manning Banco Der
coefman=seccs(ycoefm,[1;10;11;12]);
%Determinamos si existe el canal del banco izquierdo. Si existe, lo
%imprimimos junto con el canal principal. En caso de no existir, solo
%se imprime el canal principal (en ubicación "0", ya que sería desde
el
%principio
if isnan(coefman(cont,2))==0
fprintf(fileID,'%1s, %1.3f\r\n','0.00',coefman(cont,2));
fprintf(fileID,'%1.7f,
%1.3f\r\n',bancos.biz(cont,5),coefman(cont,3));
else
fprintf(fileID,'%1s, %1.3f\r\n','0.00',coefman(cont,3));
end
%Verificamos que exista el canal derecho. De ser el caso, lo
imprimimos
if isnan(coefman(cont,4))==0
fprintf(fileID,'%1.7f,
%1.3f\r\n',bancos.bder(cont,5),coefman(cont,4));
end
fprintf(fileID,'%0s\r\n','');
%Agregamos las posiciones de los diques
diqizq=[1,bancos.biz(cont,5),bancos.biz(cont,4)];
diqder=[2,bancos.bder(cont,5),bancos.bder(cont,4)];
fprintf(fileID,'%1s\r\n','LEVEE POSITIONS:');
fprintf(fileID,'%1.0f, %1.7f, %1.5f\r\n',diqizq);
fprintf(fileID,'%1.0f, %1.7f, %1.5f\r\n',diqder);
fprintf(fileID,'%0s\r\n','');
fprintf(fileID,'%1s\r\n','CUT LINE:');
coord=[dem.secerradoint(y(1),2);dem.secerradoint(y(1),3)];
fprintf(fileID,'%1.9f, %1.8f\r\n',coord);
coord=[rio.riocalint(cont,5);rio.riocalint(cont,6)];
fprintf(fileID,'%1.9f, %1.8f\r\n',coord);
coord=[dem.secerradoint(y(end),2);dem.secerradoint(y(end),3)];
fprintf(fileID,'%1.9f, %1.8f\r\n',coord);
fprintf(fileID,'%0s\r\n','');
168
fprintf(fileID,'%1s\r\n','SURFACE LINE:');
A=[dem.secerradoint(y,2)';dem.secerradoint(y,3)';dem.secerradoint(y,4)'];
fprintf(fileID,formato,A);
fprintf(fileID,'%0s\r\n','');
fprintf(fileID,'%1s\r\n','END:');
fprintf(fileID,'%0s\r\n','');
end
fprintf(fileID,'%1s\r\n','END CROSS-SECTIONS:');
fprintf(fileID,'%0s\r\n','');
%Comenzamos las líneas de los diques----------------------------
fprintf(fileID,'%1s\r\n','BEGIN LEVEES:');
%Dique izquierdo
fprintf(fileID,'%1s %1s\r\n','LEVEE ID:','1');
fprintf(fileID,'%1s\r\n','SURFACE LINE:');
diqs=bancos.biz(:,[2:4]);
fprintf(fileID,'%1.8f %1.8f
%1.8f\r\n',[diqs(:,1)';diqs(:,2)';diqs(:,3)']);
fprintf(fileID,'%1s\r\n','END:');
%Dique derecho
fprintf(fileID,'%1s %1s\r\n','LEVEE ID:','2');
fprintf(fileID,'%1s\r\n','SURFACE LINE:');
diqs=bancos.bder(:,[2:4]);
fprintf(fileID,'%1.8f %1.8f
%1.8f\r\n',[diqs(:,1)';diqs(:,2)';diqs(:,3)']);
fprintf(fileID,'%1s\r\n','END:');
fprintf(fileID,'%1s\r\n','END LEVEES:');
%termina línea de los diques-----------------------------------
fclose(fileID);
6.1.15 BANYDIQ
169
rio.dflujoint=unique(seccs(:,1));
rio.riocalint=rio.riocal(rio.dflujoint,:);
clutch=1;
end
%
for cont=1:length(rio.dflujoint)
bancos.biz(cont,2)=dem.secerradoint(max,2);
bancos.biz(cont,3)=dem.secerradoint(max,3);
bancos.biz(cont,4)=dem.secerradoint(max,4);
170
%colocamos el porcentaje que le corresponde en cuanto a la longitud
de
%la sección
Xiz=bancos.biz(cont,2);
Yiz=bancos.biz(cont,3);
diz=((Xeje-Xiz)^2+(Yeje-Yiz)^2)^0.5;
if diz>301
bancos.biz(cont,1)
Xiz
Yiz
diz
y2-cont2-2
max
end
bancos.biz(cont,5)=(lonsec/2-diz)/lonsec;
%-------------------------------------------
%Encontramos el primer máximo a la derecha
max=0;
cont2=1;
while max==0
if y2+cont2+2<=y(end)
p1=dem.secerradoint(y2+cont2,4);
p2=dem.secerradoint(y2+cont2+1,4);
p3=dem.secerradoint(y2+cont2+2,4);
else
max=y2+cont2+1;
end
if p2>=p3 & p2>p1
max=y2+cont2+1;
end
cont2=cont2+1;
end
%Colocamos los datos del banco izquierdo encontrado en una matriz
%Colocamos primero la sección a la que pertenece el banco
bancos.bder(cont,1)=dem.secerradoint(y(1));
%Colocamos las tres coordenadas del banco, x,y,z
bancos.bder(cont,2)=dem.secerradoint(max,2);
bancos.bder(cont,3)=dem.secerradoint(max,3);
bancos.bder(cont,4)=dem.secerradoint(max,4);
%colocamos el porcentaje que le corresponde en cuanto a la longitud
de
%la sección
Xder=bancos.bder(cont,2);
Yder=bancos.bder(cont,3);
dder=((Xeje-Xder)^2+(Yeje-Yder)^2)^0.5;
bancos.bder(cont,5)=(lonsec/2+dder)/lonsec;
end
%Encontramos las coordenadas de los tercios centrales
try
for cont=1:size(bancos.biz,1)
171
%bancos.biz(cont,6)=((bancos.biz(cont,1)-
bancos.biz(cont+1,1))^2+(bancos.biz(cont,2)-bancos.biz(cont+1,2))^2)^.5;
%bancos.bder(cont,6)=((bancos.bder(cont,1)-
bancos.bder(cont+1,1))^2+(bancos.bder(cont,2)-
bancos.bder(cont+1,2))^2)^.5;
ys=find(seccs(:,1)==bancos.biz(cont,1));
if bancos.biz(cont,5)~=0
xini=seccs(ys(1),2);
yini=seccs(ys(1),3);
%la distancia del centroide
dci=(2/3)*((xini-bancos.biz(cont,2))^2+(yini-
bancos.biz(cont,3))^2)^.5;
%El vector de dirección
iiz=-1*[(xini-bancos.biz(cont,2));(yini-
bancos.biz(cont,3))]/abs(((xini-bancos.biz(cont,2))^2+(yini-
bancos.biz(cont,3))^2)^.5);
xci=xini+iiz(1)*dci;
yci=yini+iiz(2)*dci;
else
%Las coordenadas del centroide
xci=bancos.biz(cont,2);
yci=bancos.biz(cont,3);
end
if bancos.bder(cont,5)~=1
xfin=seccs(ys(end),2);
yfin=seccs(ys(end),3);
dcd=(1/3)*((xfin-bancos.bder(cont,2))^2+(yfin-
bancos.bder(cont,3))^2)^.5;
ider=[(xfin-bancos.bder(cont,2));(yfin-
bancos.bder(cont,3))]/abs(((xfin-bancos.bder(cont,2))^2+(yfin-
bancos.bder(cont,3))^2)^.5);
xcd=bancos.bder(cont,2)+ider(1)*dcd;
ycd=bancos.bder(cont,3)+ider(2)*dcd;
else
xcd=bancos.bder(cont,2);
ycd=bancos.bder(cont,3);
end
bancos.biz(cont,6)=xci;
bancos.biz(cont,7)=yci;
bancos.bder(cont,6)=xcd;
bancos.bder(cont,7)=ycd;
end
catch
cont
end
%Determinamos las distancias entre los bancos
for cont=1:size(bancos.biz,1)-1
bancos.biz(cont,8)=((bancos.biz(cont,6)-
bancos.biz(cont+1,6))^2+(bancos.biz(cont,7)-bancos.biz(cont+1,7))^2)^.5;
bancos.bder(cont,8)=((bancos.bder(cont,6)-
bancos.bder(cont+1,6))^2+(bancos.bder(cont,7)-
bancos.bder(cont+1,7))^2)^.5;
172
end
6.1.16 SECCIONAR
function V = seccionar(dem,rio,d)
%Obtengo las secciones transversales a partir de un punto
%Los datos de entrada serán: Coordenadas X e Y del punto
%Dirección de la sección
%Coordenadas de Referencia
%Resolución en m de los pixeles
%Tamaño de la sección (en pixeles, por ahora)
%Matriz del dem
%Variables de prueba
%Coordenadas de Referencia
%esquina inferior izquierda
Xref=dem.ref.XWorldLimits(1);
Yref=dem.ref.YWorldLimits(1);
%Resoluciones en m
ResX=dem.ref.CellExtentInWorldX;
ResY=dem.ref.CellExtentInWorldY;
%cont=2
for cont2=1:(length(rio.riocal))
direc=rio.riocal(cont2,4);
%direccion=2;
cont3=1;
V=zeros(tam*2+1,6);
switch direc
case {1,5}
for cont=1:(tam*2+1)
if X-tam-1+cont>0 & X-tam-1+cont<=cols
%El número de la sección
V(cont3,1)=cont2;
%La coordenada en X
V(cont3,2)=Xr-ResX*(tam+1-cont);
%La Coordenada en Y
V(cont3,3)=Yr;
%El valor de la sección
V(cont3,4)=matdem(Y,X-tam-1+cont);
%Las coordenadas dentro de la matriz de la sección
V(cont3,5)=(V(cont3,2)-Xml)/ResX;
V(cont3,6)=(Yml-V(cont3,3))/ResY;
cont3=cont3+1;
end
end
case {2,6}
for cont=1:(tam*2+1)
if Y+tam+1-cont>0 & Y+tam+1-cont<=filas
if X+tam+1-cont>0 & X+tam+1-cont<=cols
%El número de la sección
V(cont3,1)=cont2;
%La coordenada en X
V(cont3,2)=Xr+ResX*(tam+1-cont);
%La Coordenada en Y
V(cont3,3)=Yr-ResY*(tam+1-cont);
%El valor de la sección
V(cont3,4)=matdem(Y+tam+1-cont,X+tam+1-cont);
%V(cont,4)=M(Y+tam+1-cont,X-tam-1+cont);
%Las coordenadas dentro de la matriz de la sección
V(cont3,5)=(V(cont3,2)-Xml)/ResX;
V(cont3,6)=(Yml-V(cont3,3))/ResY;
cont3=cont3+1;
end
end
end
case {3,7}
174
for cont=1:(tam*2+1)
if Y-tam-1+cont>0 & Y-tam-1+cont<=filas
%El número de la sección
V(cont3,1)=cont2;
%La coordenada en X
V(cont3,2)=Xr;
%La Coordenada en Y
V(cont3,3)=Yr+ResY*(tam+1-cont);
%El valor de la sección
V(cont3,4)=matdem(Y-tam-1+cont,X);
V(cont3,5)=(V(cont3,2)-Xml)/ResX;
V(cont3,6)=(Yml-V(cont3,3))/ResY;
cont3=cont3+1;
end
%Las coordenadas dentro de la matriz de la sección
end
case {4,8}
for cont=1:(tam*2+1)
if Y-tam-1+cont>0 & Y-tam-1+cont<=filas
if X+tam+1-cont>0 & X+tam+1-cont<=cols
%El número de la sección
V(cont3,1)=cont2;
%La coordenada en X
V(cont3,2)=Xr+ResX*(tam+1-cont);
%La Coordenada en Y
V(cont3,3)=Yr+ResY*(tam+1-cont);
%El valor de la sección
V(cont3,5)=(V(cont3,2)-Xml)/ResX;
V(cont3,6)=(Yml-V(cont3,3))/ResY;
V(cont3,4)=matdem(Y-tam-1+cont,X+tam+1-cont);
cont3=cont3+1;
end
end
%Las coordenadas dentro de la matriz de la sección
end
end
V2=[V2;V];
end
% V(1,1)='Sección';
% V(1,2)='X';
% V(1,3)='Y';
% V(1,4)='Z';
V=V2;
for cont=1:10
beep
pause(1)
end
%Guaradmos en un archivo la matriz con:
%dlmwrite('seccion.txt',a,'newline','pc','precision',12)
175
6.1.17 EROSIONADOR
function [seccs p
dem4]=erosionador(dem,dem2,rio2,tam,esp,vn,usosuelo,tirantes)
tp=dem2.info.PixelScale(1);
%Determinamos cuantos pixeles corresponden al tamaño seleccionado
%(Observando que la función seccionar pide como dato el número de pixeles
%a cada lado del eje)
d=ceil(tam/(2*tp));
seccnulas=seccionar(dem,rio2,d);
%Dejamos los valores de elevaciones como booleano según pertenece o no al
%espejo de agua
nulo=seccnulas(:,4);
nulo(nulo~=vn)=0;
nulo(nulo==vn)=1;
seccnulas(:,4)=nulo;
176
%Quitamos los renglones en blanco dentro de las matrices de secciones
y=find(seccsespejo(:,1)==0);
seccsespejo(y,:)=[];
y=find(seccnulas(:,1)==0);
seccnulas(y,:)=[];
% assignin('base', 'secespejo',seccsespejo);
% assignin('base', 'seccnulas',seccnulas);
% assignin('base', 'seccstraslape',seccstraslape);
%Creamos suna matriz de superficies de agua para todas las secciones (por
%ahora, solo tenemos en las secciones que no se traslapan
p=[];
ws=[];
for cont=1:size(tirantes,1)
try
if cont+1>size(tirantes,1)
ysnu=find(seccnulas(:,1)>=nseccstras(cont));
else
ysnu=find(seccnulas(:,1)>=nseccstras(cont) &
seccnulas(:,1)<nseccstras(cont+1));
end
for cont=1:nseccrio
try
ysnu=find(seccsespejo(:,1)==cont);
p(ysnu,1)=ysnu;
esp=seccsespejo(ysnu,4).*seccnulas(ysnu,4);
p(ysnu,2)=esp;
tir=ws(ysnu)-esp;
p(ysnu,3)=tir;
tirm=max(tir);
tir=tirm*seccnulas(ysnu,4);
p(ysnu,4)=tir;
seccsespejo([ysnu(1):ysnu(end)],4)=(seccsespejo([ysnu(1):ysnu(end)],4)-
tir);
catch
cont
end
end
p(:,5)=ws;
%Las seccones erosionadas
seccs.seccer=seccsespejo;
dem3=dem;
%rio3=rio2;
[rio3 dem3]=closegap2b(dem3,rio2,-10,seccs.seccer);
dem4=dem3;
177
dem4.mat=flip(dem4.demcerrado);
%Las secciones erosionadas sin traslape
seccs.secctraslape=secctraslape(dem4,rio2,600,100,usosuelo);
%El espejo de agua conocido original
Sub coordenadassecciones()
Dim fd As Office.FileDialog
Dim myFile As String, text As String, textline As String, posLat As Integer, posLong As Integer
Set fd = Application.FileDialog(msoFileDialogFilePicker)
With fd
.AllowMultiSelect = False
.Title = "Importar Coordenadas Secciones HEC-RAS"
.Filters.Clear
.Filters.Add "Archivos .geo", "*.g??"
If .Show = True Then
archivo = .SelectedItems(1)
End If
End With
'Abrimos el archivo
Open archivo For Input As #1
cont = 1
cont3 = 0
Do Until EOF(1)
Line Input #1, textline
If Mid(textline, 1, 15) = "XS GIS Cut Line" Then
Line Input #1, textline
Range("a1").Item(cont, 1) = cont
cont3 = 0
cont = cont + 1
End If
text = text & textline
Loop
Close #1
'posLat = InStr(text, "XS GIS Cut")
End Sub
Sub puntossecciones()
Dim fd As Office.FileDialog
Dim myFile As String, text As String, textline As String, posLat As Integer, posLong As Integer
Set fd = Application.FileDialog(msoFileDialogFilePicker)
With fd
.AllowMultiSelect = False
.Title = "Importar Coordenadas Secciones HEC-RAS"
.Filters.Clear
.Filters.Add "Archivos .geo", "*.g??"
If .Show = True Then
archivo = .SelectedItems(1)
End If
End With
'Abrimos el archivo
Open archivo For Input As #1
cont = 1
cont3 = 0
Do Until EOF(1)
Line Input #1, textline
If Mid(textline, 1, 9) = "#Sta/Elev" Then
179
Line Input #1, textline
'Range("a1").Item(cont, 1) = textline
'cont = cont + 1
Line Input #1, textline
cont3 = cont3 + cont2 - 1
Loop
Range("a1").Item(cont, 1) = cont
cont3 = 0
cont = cont + 1
End If
text = text & textline
Loop
Close #1
'Ahora, no todas las secciones comienzan desde cero, es por esto vamos a arreglar las secciones de
forma que las distancias sean relativas al inicio
cont = 1
Do While Range("b1").Item(cont, 1) <> ""
ini = Range("b1").Item(cont, 1)
cont2 = 1
If ini <> 0 Then
Do While Range("b1").Item(cont, cont2) <> ""
Range("b1").Item(cont, cont2) = Range("b1").Item(cont, cont2) - ini
cont2 = cont2 + 2
Loop
End If
cont = cont + 1
Loop
End Sub
Sub distancias()
'Colocamos las distancias, emparejadas con los vectores unitarios
cont = 1
Do While Worksheets(1).Range("a1").Item(cont, 1) <> ""
'Colocamos el número de la sección
180
Range("a1").Item(cont, 1) = Worksheets(1).Range("a1").Item(cont, 1)
'Colocamos las parejas "distancia-vector" (vector tiene dos coordenadas
cont2 = 2
cont3 = 2
With Worksheets(1).Range("a1")
Do While .Item(cont, cont2 + 2) <> ""
dist = ((.Item(cont, cont2) - .Item(cont, cont2 + 2)) ^ 2 + (.Item(cont, cont2 + 1) - .Item(cont,
cont2 + 3)) ^ 2) ^ 0.5
Range("a1").Item(cont, cont3) = dist
'ALgunos puntos están repetidos, así que tomamos en cuenta el caso en que la distancia sea
cero
If dist = 0 Then
Range("a1").Item(cont, cont3 + 1) = 0
Range("a1").Item(cont, cont3 + 2) = 0
Else
With Worksheets(2).Range("a1")
Do While .Item(cont, 1) <> ""
cont2 = 2
Loop
Range("a1").Item(cont4, 1) = cont
Range("a1").Item(cont4, 2) = x
Range("a1").Item(cont4, 3) = y
Range("a1").Item(cont4, 4) = h
Range("a1").Item(cont4, 5) = d2
cont4 = cont4 + 1
cont2 = cont2 + 2
Loop
'cont4 = cont4 + 1
cont = cont + 1
Loop
End With
End Sub
Sub suavizar()
'Ya que hay una gran cantidad de picos falsos, trataremos de suavizar un poco las secciones de
forma que se reduzca la detección de margenes falsas
'Elegimos el nivel de suavizado (cuantas secciones se van a utilizar en el promedio
Worksheets(5).Activate
n=5
cont = 1
cont4 = 1
'Determinamos las dimensiones
Dim ndatossec As Long
Dim nsec As Long
With Worksheets(2).Range("a1")
nsec = .Item(1, 1).Cells(Rows.Count, 1).End(xlUp).Row
'Do While .Item(cont, 1) <> ""
For cont = 1 To nsec
ndatossec = .Item(cont, 1).Cells(1, Columns.Count).End(xlToLeft).Column
cont2 = 3
'Sacamos el promedio de los datos en los que es posible (antes del final de la media
For cont2 = 3 To (ndatossec - n * 2) Step 2
'Do While .Item(cont, cont2) <> ""
182
valor = 0
For cont3 = 0 To (n - 1) * 2 Step 2
valor = valor + .Item(cont, cont2 + cont3)
Next
valor = valor / n
Range("a1").Item(cont, 1) = .Item(cont, 1)
Range("a1").Item(cont, cont2 - 1) = .Item(cont, cont2 - 1)
'Range("a1").Item(cont, (cont2 - 1) / 2 + 1) = valor
Range("a1").Item(cont, cont2) = valor
'cont2 = cont2 + 2
Next
'Obtenemos los promedios restantes
'Loop
'cont = cont + 1
Next
'Loop
End With
End Sub
Sub centro()
Worksheets(6).Activate
'Encontramos a que distancia (dentro de la sección) en la que se encuentra el centro
cont = 1
Range("e:e").Clear
With Worksheets(1).Range("a1")
'Atrás
npuntos = npuntos - 1
dist = 0
For cont2 = 1 To npuntos Step 2
dist = dist + ((.Item(cont, 1 + cont2) - .Item(cont, 3 + cont2)) ^ 2 + (.Item(cont, 2 +
cont2) - .Item(cont, 4 + cont2)) ^ 2) ^ 0.5
Next
Range("a1").Item(cont, 1) = .Item(cont, 1)
Range("a1").Item(cont, 2) = dist
Range("a1").Item(cont, 3) = .Item(cont, npuntos + 2)
Range("a1").Item(cont, 4) = .Item(cont, npuntos + 3)
If (WorksheetFunction.RoundDown(Range("i1").Item(cont, 1), 0) = 0 And
WorksheetFunction.RoundDown(Range("i1").Item(cont, 2), 0) = 0) Then Range("a1").Item(cont, 5)
=1
'Adelante. Solo sustituimos si no dio resultados adecuados en a corrida pasada
If Not Range("a1").Item(cont, 5) = 1 Then
npuntos = npuntos + 2
dist = 0
For cont2 = 1 To npuntos Step 2
dist = dist + ((.Item(cont, 1 + cont2) - .Item(cont, 3 + cont2)) ^ 2 + (.Item(cont, 2 +
cont2) - .Item(cont, 4 + cont2)) ^ 2) ^ 0.5
Next
Range("a1").Item(cont, 1) = .Item(cont, 1)
Range("a1").Item(cont, 2) = dist
Range("a1").Item(cont, 3) = .Item(cont, npuntos + 2)
Range("a1").Item(cont, 4) = .Item(cont, npuntos + 3)
End If
End If
cont = cont + 1
Loop
End With
End Sub
Sub centro2()
'Ya que dentro de el archivo de secciones no hay pistas de cual corresponde al centro
184
'Decido importar directamente otro archivo con las coordenadas de los centros
'Que fueron obtenidas a partir de hec ras->shp->idrisi->xyz (vaya proceso...)
'Abrir archivo-------------------------------------
Dim fd As Office.FileDialog
Dim myFile As String, text As String, textline As String, posLat As Integer, posLong As Integer
Set fd = Application.FileDialog(msoFileDialogFilePicker)
With fd
.AllowMultiSelect = False
.Title = "Importar Coordenadas Secciones HEC-RAS"
.Filters.Clear
.Filters.Add "Archivos .geo", "*.*"
If .Show = True Then
archivo = .SelectedItems(1)
End If
End With
'Abrimos el archivo
Open archivo For Input As #1
'Abrir archivo fin------------------------------------
Range("a1").Item(cont, 1) = cont
valor = InStr(1, textline, " ")
Range("a1").Item(cont, 2) = Mid(textline, 1, valor)
valor2 = InStr(valor + 1, textline, " ")
Range("a1").Item(cont, 3) = Mid(textline, valor, (valor2 - valor))
cont = cont + 1
Loop
185
'Importar valores-------------------------------------------------
'Cerramos el archivo
Close #1
End Sub
Sub diques()
'Encontramos la posición inicial de los diques
cont = 1
With Worksheets(5).Range("a1")
Do While .Item(cont, 1) <> ""
cont2 = cont2 + 2
Loop
cont2 = cont2 - 2
Loop
biz = cont3
'Derecha
cont3 = cont2 + 1
Do While cont3 <= (ndatossec - 2)
If .Item(cont, cont3 - 2) <= .Item(cont, cont3) And .Item(cont, cont3 + 2) <= .Item(cont,
cont3) Then
Exit Do
186
End If
cont3 = cont3 + 2
Loop
bder = cont3
'Ahora, colocamos los valores de las secciones originales
Range("a1").Item(cont, 1) = Worksheets(2).Range("a1").Item(cont, 1)
Range("a1").Item(cont, 2) = Worksheets(2).Range("a1").Item(cont, biz - 1)
Range("a1").Item(cont, 3) = Worksheets(2).Range("a1").Item(cont, biz)
Range("a1").Item(cont, 4) = Worksheets(2).Range("a1").Item(cont, bder - 1)
Range("a1").Item(cont, 5) = Worksheets(2).Range("a1").Item(cont, bder)
cont = cont + 1
Loop
End With
End Sub
Sub diques2()
'Encontramos la posición inicial de los diques. Mejora en esta versión que el centro es encontrado
dentro de la sección y no deducido a la mitad de esta
Worksheets(7).Activate
cont = 1
With Worksheets(5).Range("a1")
Do While .Item(cont, 1) <> ""
If cont = 10 Then
a=1
End If
cont2 = cont2 + 2
Loop
cont2 = cont2 - 2
Loop
If cont3 < 5 Then biz = 3 Else biz = cont3
'biz = cont3
'Derecha
cont3 = cont2 + 1
Do While cont3 <= (ndatossec - 2)
If .Item(cont, cont3 - 2) <= .Item(cont, cont3) And .Item(cont, cont3 + 2) <= .Item(cont,
cont3) Then
Exit Do
End If
cont3 = cont3 + 2
Loop
Range("a1").Item(cont, 1) = Worksheets(2).Range("a1").Item(cont, 1)
Range("a1").Item(cont, 2) = Worksheets(2).Range("a1").Item(cont, biz - 1)
Range("a1").Item(cont, 3) = Worksheets(2).Range("a1").Item(cont, biz)
Range("a1").Item(cont, 4) = Worksheets(2).Range("a1").Item(cont, bder - 1)
Range("a1").Item(cont, 5) = Worksheets(2).Range("a1").Item(cont, bder)
cont = cont + 1
Loop
End With
End Sub
Sub perfiltra()
'Esta función quita los picos del eje
cont = 1
Do While Range("b2").Item(cont + 1, 1) <> ""
p1 = Range("b2").Item(cont, 1)
p2 = Range("b2").Item(cont + 1, 1)
If p1 < p2 Then
cont2 = 1
'Encontramos el siguiente valor menor
188
Do While p1 < p2
cont2 = cont2 + 1
p2 = Range("b2").Item(cont + cont2, 1)
Loop
'Corregimos todos los valores entre las secciones correctas
Range("d2").Item(cont, 1) = p1
dif = p1 - p2
inc = dif / cont2
cont4 = 1
For cont3 = cont + 1 To cont + cont2
cont = cont + 1
Loop
End Sub
Sub nombres()
'Obtenemos los nombres de las secciones
Dim fd As Office.FileDialog
Dim myFile As String, text As String, textline As String, posLat As Integer, posLong As Integer
Set fd = Application.FileDialog(msoFileDialogFilePicker)
With fd
.AllowMultiSelect = False
.Title = "Importar Coordenadas Secciones HEC-RAS"
.Filters.Clear
.Filters.Add "Archivos .geo", "*.g??"
If .Show = True Then
archivo = .SelectedItems(1)
End If
End With
'Abrimos el archivo
Open archivo For Input As #1
cont = 1
cont2 = 0
textline = ""
Do Until EOF(1)
Line Input #1, textline
If Mid(textline, 1, 14) = "Type RM Length" Then
Range("a1").Item(cont, 1) = nombre
cont = cont + 1
End If
cont2 = cont2 + 1
Range("b1") = cont2
'text = text & textline
Loop
Close #1
End Sub
Sub coordiques()
'Agregamos las coordenadas a los diques
Worksheets(7).Activate
cont = 1
cont2 = 1
Do While Range("a1").Item(cont, 1) <> ""
iz = Range("a1").Item(cont, 2)
der = Range("a1").Item(cont, 4)
'cont2 = 1
Do While Worksheets(4).Range("a1").Item(cont2, 1) <> cont &
Worksheets(4).Range("a1").Item(cont2, 1) <> ""
cont2 = cont2 + 1
Loop
Do While Worksheets(4).Range("a1").Item(cont2, 5) <> "" And
Worksheets(4).Range("a1").Item(cont2, 1) = Range("a1").Item(cont, 1) And
Worksheets(4).Range("a1").Item(cont2, 5) < iz
cont2 = cont2 + 1
Loop
190
Range("g1").Item(cont, 1) = Worksheets(4).Range("a1").Item(cont2, 2)
Range("g1").Item(cont, 2) = Worksheets(4).Range("a1").Item(cont2, 3)
cont = cont + 1
Loop
End Sub
Conjunto de funciones programadas en visual basic que ejecutan las pruebas estadísticas requeridas
en el proyecto
Sub hatanaca()
'Realiza la prueba de hatanaca para los datos seleccionados
'quitamos los valores de la columna de resultados
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Prueba de Hatanaca"
With Range("i3")
.Item(1, 1) = "Ns"
.Item(2, 1) = "Ni"
.Item(3, 1) = "w"
.Item(4, 1) = "e"
.Item(5, 1) = "Media"
.Item(6, 1) = "Sigma"
.Item(7, 1) = "Za/2"
.Item(8, 1) = "Resultado"
End With
cont1 = 1
Ns = 0
Ni = 0
Set a = Range("d2")
Set b = a
xj = Range("d2").Item(cont1 + 1, 1)
'Set a = Range("d2")
'Set b = a
'For c = 1 To cont1
' Set b = Union(b, a.Offset(c, 0))
191
'Next
rango = "d2:d" & cont1 + 2
Set b = Range(rango)
If xj = WorksheetFunction.Max(b) Then
Range("d2").Item(cont1 + 1, 3) = "S"
Ns = Ns + 1
ElseIf xj = WorksheetFunction.Min(b) Then
Range("d2").Item(cont1 + 1, 3) = "i"
Ni = Ni + 1
End If
cont1 = cont1 + 1
Loop
w = Ns - Ni
e = w / (2 * WorksheetFunction.Ln(cont1) - 0.845)
med = WorksheetFunction.Average(b)
desv = WorksheetFunction.StDev_S(b)
Z = Abs(WorksheetFunction.Norm_S_Inv(1 - Range("j2") / 200))
Range("j3") = Ns
Range("j4") = Ni
Range("j5") = w
Range("j6") = e
Range("j7") = med
Range("j8") = desv
Range("j9") = Z
'damos un veredicto
If (-Z <= e And e <= Z) Then
Range("j10") = "La muestra se considera aleatoria"
Else
Range("j10") = "La muestra no se considera aleatoria"
End If
End Sub
Sub waldywolf()
'Prueba de wald y wolfowitz
'limmpiamos la prueba anterior
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Prueba de Wald y Wolfowitz"
With Range("i3")
192
.Item(1, 1) = "R"
.Item(2, 1) = "Rachas+"
.Item(3, 1) = "Rachas-"
.Item(4, 1) = "medk"
.Item(5, 1) = "Vark"
.Item(6, 1) = "Z"
.Item(7, 1) = "Z crit"
.Item(8, 1) = "Resultado"
End With
cont = 1
Do While Range("d2").Item(cont, 1) <> ""
cont = cont + 1
Loop
mediana = WorksheetFunction.Median(b)
cont1 = 1
Else
Range("d2").Item(cont1, 3) = "-"
End If
cont1 = cont1 + 1
Loop
n = cont1
nmas = 0
nmenos = 0
cont1 = 1
193
End If
cont1 = cont1 + 1
Loop
cont1 = 1
'rachas -
Do While Range("f2").Item(cont1, 1) <> ""
End If
cont1 = cont1 + 1
Loop
alfa = Range("j2")
r = nmas + nmenos
medk = (2 * nmas * nmenos) / r + 1
vark = (2 * nmas * nmenos * (2 * nmas * nmenos - nmas - nmenos)) / (((nmas + nemenos) ^ 2) *
(nmas + nmenos - 1))
ze = (r - medk + 0.5) / (vark ^ 0.5)
zcrit = WorksheetFunction.Norm_S_Inv(1 - alfa / 200)
Range("j3") = r
Range("j4") = nmas
Range("j5") = nmenos
Range("j6") = medk
Range("j7") = vark
Range("j8") = ze
Range("j9") = zcrit
End Sub
194
Sub testiter()
'realiza el test de iteración para homogeneidad
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Prueba de Iteración"
With Range("i3")
.Item(1, 1) = "Ni0"
.Item(2, 1) = "Ni1"
.Item(3, 1) = "N0"
.Item(4, 1) = "N1"
.Item(5, 1) = "Xw"
.Item(6, 1) = "S2w"
.Item(7, 1) = "w"
.Item(8, 1) = "e"
.Item(9, 1) = "Za/2"
.Item(10, 1) = "Resultado"
End With
cont = 1
Do While Range("d2").Item(cont, 1) <> ""
cont = cont + 1
Loop
rango = "d2:d" & cont
Set b = Range(rango)
med = WorksheetFunction.Average(b)
n = cont - 1
'recomenzamos el conteo, porque el primer bucle era para determinar la media
cont1 = 1
ceros = 0
unos = 0
ceros = ceros + 1
Range("d2").Item(cont1, 3) = 0
Else
unos = unos + 1
Range("d2").Item(cont1, 3) = 1
End If
cont1 = cont1 + 1
195
Loop
Range("j5") = ceros
Range("j6") = unos
'Iteraciones de 0
cont1 = 1
ni0 = 0
ni1 = 0
If Range("f2").Item(cont1, 1) = 0 Then
ni0 = ni0 + 1
End If
cont1 = cont1 + 1
Loop
Range("j3") = ni0
'iteraciones de 1
cont1 = 1
If Range("f2").Item(cont1, 1) = 1 Then
ni1 = ni1 + 1
End If
cont1 = cont1 + 1
Loop
Range("j4") = ni1
196
alfa = Range("j2")
Xw = (2 * ceros * unos) / n + 1
s2w = ((2 * ceros * unos) * (2 * ceros * unos - n)) / (n ^ 2 * (n - 1))
w = ni0 + ni1
e = (w - Xw) / (s2w ^ 0.5)
Z = Abs(WorksheetFunction.Norm_S_Inv(1 - alfa / 200))
Range("j7") = Xw
Range("j8") = s2w
Range("j9") = w
Range("j10") = e
Range("j11") = Z
End If
End Sub
Sub mannwitney()
'Ejecuta una prueba de mann Witney
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Prueba de Mann Witney"
With Range("i3")
.Item(1, 1) = "n1"
.Item(2, 1) = "n2"
.Item(3, 1) = "R1"
.Item(4, 1) = "R2"
.Item(5, 1) = "U1"
.Item(6, 1) = "U2"
.Item(7, 1) = "U"
.Item(8, 1) = "med"
.Item(9, 1) = "Var"
.Item(10, 1) = "Desv"
.Item(11, 1) = "Z"
.Item(12, 1) = "Ucrit"
.Item(13, 1) = "p"
.Item(14, 1) = "Resultado"
End With
cont = 1
197
Do While Range("d2").Item(cont, 1) <> ""
cont = cont + 1
Loop
cont1 = 1
n = cont1
'determinamos si son pares o impares el número de datos
'0 par, 1 impar
If WorksheetFunction.Even(n) Then
par = 0
n1 = n / 2
Else
par = 1
n1 = (n - 1) / 2
End If
n2 = n - n1
For cont1 = 1 To n1
r1 = Range("e2").Item(cont1, 1) + r1
Next
For cont1 = n1 + 1 To n
r2 = Range("e2").Item(cont1, 1) + r2
Next
u1 = n1 * n2 + (n1 * (n1 + 1) / 2) - r1
u2 = n1 * n2 + (n2 * (n2 + 1) / 2) - r2
alfa = Range("j2")
u = WorksheetFunction.Min(u1, u2)
med = n1 * n2 / 2
vari = med * (n1 + n2 + 1) / 6
desv = vari ^ 0.5
ze = (Abs(u - med) - 0.5) / desv
Range("j3") = n1
Range("j4") = n2
Range("j5") = r1
Range("j6") = r2
Range("j7") = u1
Range("j8") = u2
Range("j9") = u
Range("j10") = med
Range("j11") = vari
Range("j12") = desv
Range("j13") = ze
Range("j14") = ucrit
Range("j15") = p
End Sub
Sub waldywolf2()
'Esta prueba viene descrita en la tesis AFA, no concuerda con la que encontré en la bibliografía,
pero la utilizaré
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Prueba de Wald y Wolfowitz"
With Range("i3")
.Item(1, 1) = "W"
.Item(2, 1) = "Xw"
.Item(3, 1) = "Sw^2"
.Item(4, 1) = "e"
.Item(5, 1) = "Za/2"
.Item(6, 1) = "Resultado"
End With
cont = 1
With Range("d2")
Do While .Item(cont, 1) <> ""
sumw = sumw + .Item(cont, 1) * .Item(cont + 1, 1)
t1 = t1 + .Item(cont, 1)
199
t2 = t2 + .Item(cont, 1) ^ 2
t3 = t3 + .Item(cont, 1) ^ 3
t4 = t4 + .Item(cont, 1) ^ 4
cont = cont + 1
Loop
cont = cont - 1
w = sumw + .Item(cont, 1) * .Item(1, 1)
Xw = (t1 ^ 2 - t2) / (cont - 1)
sw2 = (t2 ^ 2 - t4) / (cont - 1) + (t1 ^ 4 - 4 * t1 ^ 2 * t2 + 4 * t1 * t3 + t2 ^ 2 - 2 * t4) / ((cont - 1) *
(cont - 2)) - Xw ^ 2
e = (w - Xw) / (sw2 ^ 0.5)
alfa = Range("j2")
Z = Abs(WorksheetFunction.Norm_S_Inv(1 - alfa / 200))
Range("j3") = w
Range("j4") = Xw
Range("j5") = sw2
Range("j6") = e
Range("j7") = Z
If e > -Z And e < Z Then
Range("j8") = "La prueba se condiera independiente"
Else
Range("j8") = "La prueba no se considera independiente"
End If
End With
End Sub
Sub corror()
'ejecuta una prueba de correlación ordenada
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Correlación Ordenada"
With Range("i3")
.Item(1, 1) = "n"
.Item(2, 1) = "rs"
.Item(3, 1) = "t"
.Item(4, 1) = "tcrit"
.Item(5, 1) = "Resultado"
End With
cont = 1
Do While Range("d2").Item(cont, 1) <> ""
cont = cont + 1
Loop
cont1 = 1
alfa = Range("j2")
n = cont1
rs = 1 - (6 / (n * (n ^ 2 - 1))) * sumdicua
t = (((rs ^ 2) * (n - 2)) / (1 - rs ^ 2)) ^ 0.5
tcrit = WorksheetFunction.T_Inv(1 - alfa / 100, n - 2)
Range("j3") = n
Range("j4") = rs
Range("j5") = t
Range("j6") = tcrit
End Sub
Sub testinv()
'realiza una prueba de inversión a los datos
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Test de Inversión"
With Range("i3")
.Item(1, 1) = "SumInv"
.Item(2, 1) = "InvCrit"
.Item(3, 1) = "Resultados"
End With
201
cont = 1
Do While Range("d2").Item(cont, 1) <> ""
cont = cont + 1
Loop
cont1 = 1
cont2 = 1
inv = 0
suminv = 0
Do While Range("d2").Item(cont2 + 1, 1) <> ""
Do While Range("d2").Item(cont1 + 1, 1) <> ""
End If
cont1 = cont1 + 1
Loop
Range("d2").Item(cont2, 2) = inv
suminv = suminv + inv
inv = 0
cont2 = cont2 + 1
cont1 = cont2
Loop
n = cont1
alfa = Range("j2")
invcrit = ((n - 1) * n / 2) / 1.7
Range("j3") = suminv
Range("j4") = invcrit
'establecer resultado
If invcrit >= suminv Then
Range("j5") = "Datos independientes"
Else
Range("j5") = "Datos dependientes"
End If
End Sub
202
Sub geary()
'hace un test de normalidad de geary
Range("f1:f100000").Clear
Range("i3:j16").Clear
Range("f1") = "Prueba Geary"
With Range("i3")
.Item(1, 1) = "aj"
.Item(2, 1) = "Z"
.Item(3, 1) = "p"
.Item(4, 1) = "Resultado"
End With
cont = 1
Do While Range("d2").Item(cont, 1) <> ""
cont = cont + 1
Loop
n = cont - 1
med = WorksheetFunction.Average(b)
cont1 = 1
Loop
s1 = (suma / n)
s2 = (suma2 / n) ^ 0.5
Range("j3") = mu
Range("j4") = Z
Range("j5") = p
If p > mu Then
Range("j6") = "Se considera normal la muestra"
Else
Range("j6") = "La muestra no se consdiera normal"
End If
End Sub
203
6.1.20 EXPORTAR A CEQUEAU
Un conjunto de herramientas que prepara los datos hidrometeorológicos y los exporta a formato
legible para CEQUEau
Module completardiasfaltantes
Function compserie2(dir, vnulo) As ArrayList
'La función de esta aplicación es la de llenar los valores faltantes de algún
año
'colocamos una salida en caso de error
On Error GoTo er
'Aquí daré una dirección fija, en el programa real, será una variable
'dir = "C:\Users\Hidding\Box Sync\Box
Sync\Maestría\Tesis\Programas\Matlab\Pruebas.xlsx"
204
wb = exapp.Workbooks.Open(dir)
ws = wb.Worksheets(1)
With ws.Range("a9")
precipitación.Add(a)
Else
precipitación.Add(vnulo)
End If
fechag.AddDays(1)
cont2 = cont2 + 1
205
cont = cont + 1
End While
'Cada año debe tener 366 días, sea o no bisiesto. Aquí completamos
'los años no bisiestos con un valorextra
Dim numaños As Integer
Dim cont3 As Integer
numaños = añofin - añoini
Next
precipitacion.Insert(0, añoini)
precipitacion.Insert(1, añofin)
'colocamos el nombre de la estación a 8 caracteres
Dim estacion As String
estacion = ws.Range("B2").Text
estacion = estacion.Replace(" ", "")
precipitacion.Insert(2, estacion)
'intervalo de horas
precipitacion.Insert(3, 24)
'precensia de nieve (1 es ausencia)
precipitacion.Insert(4, 1)
'verificar este valor
206
precipitacion.Insert(5, 0)
End Function
207
End If
'agregamos el encabezado
' If cont Mod 366 = 0 Then
cadena = " " & tmin.item(2) & " " & (tmin.item(0) + (cont - 6) / 366) &
" " & tmin.item(3) & " " & tmin.item(4) & (" ") & tmin.item(5)
ordenar.Add(cadena)
cadena = ""
'End If
End While
'regresamos el valor que obtuvimos.
agrupar = ordenar
End Function
End Module
6.1.21 EXTHID
Sub exthid()
'Cambiamos el paso del hidrograma para poder compararlo
'ps es el paso de tiempo en minutos
Range("AT:CJ").Clear
ps = 10
cont = 1
cont2 = 1
cont4 = 1
cont5 = 0
209
'Range("d2").Item(cont2, 43 + cont5) = Range("a2").Item(cont, 1 + cont5) * 60
'Range("e2").Item(cont2, 43 + cont5) = qfin
'cont2 = cont2 + 1
cont = cont + 1
Loop
cont5 = cont5 + 2
Loop
'Call desconvol
Call Quitar2
End Sub
Los archivos *.geo son archivos de texto en un formato en el que HEC-RAS pueda importar toda la
información necesaria para un proyecto determinada. Ya que se están colocando los datos calculados
en este formato, se detalla cual es la estructura de tal archivo
%---------------------------Comienzan encabezados----------------------------------%
BEGIN HEADER:
END HEADER:
%-------------------------------Terminan encabezados-------------------------------%
REACH:
%Una Matriz de los datos del tramo (el eje del tramo)
X, Y, Z, Distancia al inicio
[----------]
END:
CROSS-SECTION:
%%%%%%%%%%%%%
Se colocan las posiciones de los diques, para evitar que el agua salga del canal principal sin rebasar
los bancos primero
Para esto, se necesitan tres datos:
El grupo de dique se define más adelante en el archivo, por medio de coordenadas X,Y,Z
ordenadas, que trazan una línea tridimensional que delimita el dique. En nuestro caso, hemos
usado dos de estos grupos, uno para la margen izquierda, y el otro para la derecha, en líneas
continuas por todo lo largo del rio
Posición del dique se define en una fracción decimal, donde 0 es el inicio de la sección y 1 es el final
de esta.
La elevación del dique se hace coincidir con la del terreno natural en este caso.
%%%%%%%%%%%%%
LEVEE POSITIONS:
%datos de ejemplo
211
1, 0.1750000, 361.00000
2, 0.9750000, 461.00000
CUT LINE:
X, Y
X, Y
X, Y
END CROSS-SECTIONS:
%-------------------------------Terminan Secciones transversales-------------------------------%
%-------------------------------Comienzan Diques-------------------------------%
BEGIN LEVEES:
LEVEE ID: Número de identificación del grupo de puntos de la línea que traza el dique.
SURFACE LINE:
%Se coloca cada punto de la línea del dique
X,Y,Z
[------se repite para cada punto del dique ----]
END:
[-------Se repite para cada dique-------]
END LEVEES:
%-------------------------------Terminan Diques-------------------------------%
212