Este cuaderno tiene como objetivo mostrar un pequeño ejemplo de la obtención de curvas Intensidad Duración Frecuencia sintéticas mediante el método de Diaz-Granados & Vargas. Y la estimación de caudales de diseño para diferentes periodos de retorno utilizando el 1) método racional y 2) el método del número de curva y convolución de caudales o hidrograma unitario.
Para el ejercicio se utilizarán los siguientes paquetes.
rm(list = ls())
library(dplyr) # Funciones auxiliares
library(terra) # Datos geográficos
library(raster) # Datos geográficos
library(sf) # Datos geográficos
library(leaflet) # Visor
library(phylin) # IDW
library(ggplot2) # Gráficas
library(xts) # Series temporales
# setwd("~/Diego/02.Profesionales/UPJ/")
work.dir <- "~/Diego/02.Profesionales/UPJ/"Nota: En cuadernos de R se suelen presentar problemas por el directorio de trabajo por lo que se guardó en la variable work.dir. Se recomienda que quien vaya a hacer uso de este cuaderno revise muy bien las rutas que se usan aca, la información necesaria se encuentra en este enlace
1 Zona de estudio
Para el ejemplo se utilizará una cuenca alta en los Andes colombianos, en el municipio de Guatavita, Cundinamarca.
cnca <- vect(paste0(work.dir, "GIS/CUENCA.shp")) # Lee la cuenca
drnjs <- vect(paste0(work.dir, "GIS/DRENAJES.shp")) # Drenajes
# Estensión de la cuenca en coordenadas geográficas
cnca.ext <- ext(project(cnca, "epsg:4326"))
# Cambia el tipo de variable para la visualización
cncaSpatial <- as(project(cnca, "epsg:4326"), "Spatial")
drnjSpatial <- as(project(drnjs, "epsg:4326"), "Spatial")
# Visualización interactiva
leaflet() %>% addTiles() %>%
setView(lat = mean(cnca.ext[3:4]), lng = mean(cnca.ext[1:2]), zoom = 15) %>%
addPolygons(data = cncaSpatial, color = "#830000", fillOpacity = 0.1, popup = "Cuenca") %>%
addPolylines(data = drnjSpatial, color = "#008B8B", weight = 3, fillOpacity = 1)cnca.tab <- as.data.frame(cnca)
cnca.tab[,1:5] <- round(cnca.tab[,1:5])
knitr::kable(cnca.tab, digits = 3, caption = 'Atributos de la Cuenca')| fid | Area | Perimetro | L | CMax | CMin | S | Sw |
|---|---|---|---|---|---|---|---|
| 1 | 3057067 | 11224 | 2169 | 3166 | 2784 | 0.115 | 19.778 |
2 Información de Entrada
Para el ejercicio se necesitan información de precipitación proveniente de estaciones del IDEAM a través de la plataforma DHIME; así como de coberturas y textura de los suelos proveniente del IDEAM y el IGAC, respectivamente.
La información de coberturas se agrego al segundo nivel y la de suelos por la textura. Posteriormente se cruzaron y se les asignó un coeficiente de escorrentía y un número de curva.
Nota: el catalogo nacional de estaciones que se carga en esta sección no se encuentra en el enlace compartido previamente, debe ser descargado del DHIME, en la sección de ‘Recursos’ del respectivo enlace.
# Catalogo nacional de estaciones
stn <- vect("~/Diego/04.InformaciónGeneral/Cartografia Colombia/SHP/SHP_IDEAM_ESTACIONES/CNE_IDEAM/CNE_IDEAM.shp")
# Información de precipitación descargada
data.dhime <- read.csv(paste0(work.dir, "excel.csv.csv"), encoding = "UTF-8")
data.dhime$CodigoEstacion <- as.character(data.dhime$CodigoEstacion)
# Filtro de las estaciones
stn <- stn[stn$CODIGO %in% unique(data.dhime$CodigoEstacion), ]
# Coberturas y textura de suelos
suelos <- vect(paste0(work.dir, "GIS/Cobertura&Suelos.shp"))
stnSpatial <- as(project(stn, "epsg:4326"), "Spatial")
sueloSpatial <- as(project(suelos, "epsg:4326"), "Spatial")
leaflet() %>% addTiles() %>%
setView(lat = mean(cnca.ext[3:4]), lng = mean(cnca.ext[1:2]), zoom = 13) %>%
addMarkers(data = stnSpatial, popup = ~nombre) %>%
addPolygons(data = sueloSpatial, fillOpacity = 0.7, weight = 2, color = "#8B7D7B", fillColor = "#FFE4E1",
popup = ~sprintf("Cobertura: %s\nTextura: %s", nivel_2, Textura)) %>%
addPolygons(data = cncaSpatial, color = "#830000", fillOpacity = 0.1, popup = "Cuenca") suelos.tab <- as.data.frame(suelos)
knitr::kable(suelos.tab, digits = 3, caption = 'Cobertutra y Textura de los Suelos')| codigo | leyenda | nivel_1 | nivel_2 | PAISAJE | CARACTER_1 | Textura | CN | Ce |
|---|---|---|---|---|---|---|---|---|
| 242 | 2.4.2. Mosaico de pastos y cultivos | 2. Territorios agrícolas | 2.4. Áreas agrícolas heterogéneas | Montaña | Suelos profundos a superficiales, bien a excesivamente drenados, con texturas finas a moderadamente gruesas, reacción extremada a muy fuertemente ácida, mediana saturación de aluminio y fertilidad moderada a baja | Media | 79 | 0.42 |
| 242 | 2.4.2. Mosaico de pastos y cultivos | 2. Territorios agrícolas | 2.4. Áreas agrícolas heterogéneas | Montaña | Suelos profundos a moderadamente profundos, bien drenados, con texturas medias a moderadamente gruesas, reacción muy fuerte a medianamente ácida, baja a media saturación de aluminio y fertilidad baja a moderada | Gruesa | 69 | 0.22 |
| 231 | 2.3.1. Pastos limpios | 2. Territorios agrícolas | 2.3. Pastos | Montaña | Suelos profundos a superficiales, bien a excesivamente drenados, con texturas finas a moderadamente gruesas, reacción extremada a muy fuertemente ácida, mediana saturación de aluminio y fertilidad moderada a baja | Media | 71 | 0.42 |
| 231 | 2.3.1. Pastos limpios | 2. Territorios agrícolas | 2.3. Pastos | Montaña | Suelos profundos a moderadamente profundos, bien drenados, con texturas medias a moderadamente gruesas, reacción muy fuerte a medianamente ácida, baja a media saturación de aluminio y fertilidad baja a moderada | Gruesa | 58 | 0.42 |
| 3231 | 3.2.3.1. Vegetación secundaria alta | 3. Bosques y áreas seminaturales | 3.2. Áreas con vegetación herbácea y/o arbustiva | Montaña | Suelos profundos a superficiales, bien a excesivamente drenados, con texturas finas a moderadamente gruesas, reacción extremada a muy fuertemente ácida, mediana saturación de aluminio y fertilidad moderada a baja | Media | 70 | 0.42 |
| 3231 | 3.2.3.1. Vegetación secundaria alta | 3. Bosques y áreas seminaturales | 3.2. Áreas con vegetación herbácea y/o arbustiva | Montaña | Suelos profundos a moderadamente profundos, bien drenados, con texturas medias a moderadamente gruesas, reacción muy fuerte a medianamente ácida, baja a media saturación de aluminio y fertilidad baja a moderada | Gruesa | 56 | 0.22 |
| 31121 | 3.1.1.2.1. Bosque denso bajo de tierra firme | 3. Bosques y áreas seminaturales | 3.1. Bosques | Montaña | Suelos profundos a superficiales, bien a excesivamente drenados, con texturas finas a moderadamente gruesas, reacción extremada a muy fuertemente ácida, mediana saturación de aluminio y fertilidad moderada a baja | Media | 73 | 0.50 |
| 31121 | 3.1.1.2.1. Bosque denso bajo de tierra firme | 3. Bosques y áreas seminaturales | 3.1. Bosques | Montaña | Suelos profundos a moderadamente profundos, bien drenados, con texturas medias a moderadamente gruesas, reacción muy fuerte a medianamente ácida, baja a media saturación de aluminio y fertilidad baja a moderada | Gruesa | 60 | 0.30 |
knitr::kable(head(data.dhime), digits = 3, caption = "Datos de Precipitación")| CodigoEstacion | NombreEstacion | Latitud | Longitud | Altitud | Categoria | Entidad | AreaOperativa | Departamento | Municipio | FechaInstalacion | FechaSuspension | IdParametro | Etiqueta | DescripcionSerie | Frecuencia | Fecha | Valor | Grado | Calificador | NivelAprobacion |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 35060200 | AMOLADERO EL [35060200] | 4.858 | -73.745 | 2963 | Pluviométrica | INSTITUTO DE HIDROLOGIA METEOROLOGIA Y ESTUDIOS AMBIENTALES | Area Operativa 11 - Cundinamarca-Amazonas-San Andrés | Cundinamarca | Guatavita | 15/04/1972 00:00 | NA | PRECIPITACION | PTPM_CON | Día pluviométrico (convencional) | Diaria | 1972-05-01 00:00 | 2 | 50 | 900 | |
| 35060200 | AMOLADERO EL [35060200] | 4.858 | -73.745 | 2963 | Pluviométrica | INSTITUTO DE HIDROLOGIA METEOROLOGIA Y ESTUDIOS AMBIENTALES | Area Operativa 11 - Cundinamarca-Amazonas-San Andrés | Cundinamarca | Guatavita | 15/04/1972 00:00 | NA | PRECIPITACION | PTPM_CON | Día pluviométrico (convencional) | Diaria | 1972-05-02 00:00 | 18 | 50 | 900 | |
| 35060200 | AMOLADERO EL [35060200] | 4.858 | -73.745 | 2963 | Pluviométrica | INSTITUTO DE HIDROLOGIA METEOROLOGIA Y ESTUDIOS AMBIENTALES | Area Operativa 11 - Cundinamarca-Amazonas-San Andrés | Cundinamarca | Guatavita | 15/04/1972 00:00 | NA | PRECIPITACION | PTPM_CON | Día pluviométrico (convencional) | Diaria | 1972-05-03 00:00 | 15 | 50 | 900 | |
| 35060200 | AMOLADERO EL [35060200] | 4.858 | -73.745 | 2963 | Pluviométrica | INSTITUTO DE HIDROLOGIA METEOROLOGIA Y ESTUDIOS AMBIENTALES | Area Operativa 11 - Cundinamarca-Amazonas-San Andrés | Cundinamarca | Guatavita | 15/04/1972 00:00 | NA | PRECIPITACION | PTPM_CON | Día pluviométrico (convencional) | Diaria | 1972-05-04 00:00 | 0 | 50 | 900 | |
| 35060200 | AMOLADERO EL [35060200] | 4.858 | -73.745 | 2963 | Pluviométrica | INSTITUTO DE HIDROLOGIA METEOROLOGIA Y ESTUDIOS AMBIENTALES | Area Operativa 11 - Cundinamarca-Amazonas-San Andrés | Cundinamarca | Guatavita | 15/04/1972 00:00 | NA | PRECIPITACION | PTPM_CON | Día pluviométrico (convencional) | Diaria | 1972-05-05 00:00 | 10 | 50 | 900 | |
| 35060200 | AMOLADERO EL [35060200] | 4.858 | -73.745 | 2963 | Pluviométrica | INSTITUTO DE HIDROLOGIA METEOROLOGIA Y ESTUDIOS AMBIENTALES | Area Operativa 11 - Cundinamarca-Amazonas-San Andrés | Cundinamarca | Guatavita | 15/04/1972 00:00 | NA | PRECIPITACION | PTPM_CON | Día pluviométrico (convencional) | Diaria | 1972-05-06 00:00 | 5 | 50 | 900 |
La estructura en la que vienen los datos es algo engorrosa, por lo cual es bueno reorganizarlos para que sea más como trabajar con estos. Para esto se usará una variable tipo xts.
data.precip <- NULL # Variable donde se guardará la información de precipitación
# Ciclo por cada una de las estaciones en el archivo del DHIME
for (i in unique(data.dhime$CodigoEstacion)){
di <- subset(data.dhime, CodigoEstacion == i) # Extrae la información correspondiente a la estación i
di <- as.xts(di$Valor, order.by = as.Date(di$Fecha)) # La convierte en xts
data.precip <- cbind(data.precip, di) # Crea progresivamente la matriz con una columna por estación
}
colnames(data.precip) <- unique(data.dhime$CodigoEstacion)
precip.tab <- as.data.frame(tail(data.precip))
knitr::kable(precip.tab, digits = 3, caption = 'Precipitación Ordenada')| 35060200 | 35060160 | 35060020 | |
|---|---|---|---|
| 2022-08-26 | 6.8 | 6.5 | 3.8 |
| 2022-08-27 | 0.8 | 3.0 | 2.2 |
| 2022-08-28 | 2.4 | 1.5 | 0.0 |
| 2022-08-29 | 5.4 | 2.5 | 0.0 |
| 2022-08-30 | 0.8 | 4.0 | 0.0 |
| 2022-08-31 | 0.2 | 0.0 | 0.0 |
3 Calculo de las IDF
Para el cálculo de las IDF se va a utilizar el método de Diaz-Granados & Vargas, el cual da la intensidad \(I\) en mm/hr de acuerdo con la siguiente ecuación.
\[I=\frac{a T^{b} M^{d}}{(t/60)^{c}}\]
Donde, \(T\) es el periodo de retorno en años, \(M\) en el promedio anual de la precipitación máxima en 24 horas, \(t\) el tiempo de la lluvia en minutos y \(a\), \(b\), \(c\) y \(d\) son coeficientes del método que dependen de la región en la que se trabaje. Para la región andina estos coeficientes son: 0.94, 0.18, 0.66 y 0.83, respectivamente.
Inicialmente, se calcula \(M\) para cada estación como el promedio de los máximos anuales. Interpolando espacialmente mediante IDW se obtuvo el \(M\) para la cuenca de estudio.
# Función para máximo anual, con más de 30% de datos faltantes devuelve NA.
maxAnual <- function(x){
if(sum(!is.na(x)) < 0.3*365){
return(NA)
}else{
return(max(x, na.rm = T))
}
}
# Máximo anual
data.precip.anual <- aggregate(data.precip, by = lubridate::year(index(data.precip)), maxAnual) %>%
tail(30)
knitr::kable(as.data.frame(data.precip.anual), digits = 1, caption = "Precipitación Máxima en 24 Horas Anual")| 35060200 | 35060160 | 35060020 | |
|---|---|---|---|
| 1993 | 34.3 | 31.0 | 26.8 |
| 1994 | 43.3 | 45.0 | 45.0 |
| 1995 | 43.3 | 41.0 | 35.8 |
| 1996 | 59.5 | 41.0 | 45.0 |
| 1997 | 66.6 | 50.0 | 35.0 |
| 1998 | 50.8 | 44.9 | 31.5 |
| 1999 | 36.8 | 38.2 | 40.1 |
| 2000 | 40.8 | 55.0 | 34.2 |
| 2001 | 40.9 | 50.9 | 22.4 |
| 2002 | 54.5 | 55.6 | NA |
| 2003 | 45.3 | 32.0 | 23.2 |
| 2004 | 43.6 | 40.5 | NA |
| 2005 | 36.2 | 32.8 | 22.0 |
| 2006 | 33.5 | 35.6 | 28.0 |
| 2007 | 37.3 | 35.0 | 28.3 |
| 2008 | 34.4 | 40.0 | 39.8 |
| 2009 | 33.3 | 28.8 | 43.3 |
| 2010 | 33.0 | 47.7 | 35.2 |
| 2011 | 62.5 | 51.0 | 47.5 |
| 2012 | 40.5 | 50.0 | 45.2 |
| 2013 | 51.5 | 33.5 | 58.5 |
| 2014 | 50.5 | 40.0 | 32.3 |
| 2015 | 46.0 | 29.3 | 40.0 |
| 2016 | 40.2 | 31.0 | 38.7 |
| 2017 | 54.8 | 33.0 | 42.2 |
| 2018 | 38.6 | 29.5 | 44.7 |
| 2019 | 64.8 | 38.0 | 37.0 |
| 2020 | 41.5 | 34.4 | 30.5 |
| 2021 | 45.5 | 61.1 | 36.2 |
| 2022 | 46.2 | 30.5 | 49.6 |
# Promedio de los máximos anuales
data.precip.anual <- colMeans(data.precip.anual, na.rm = T)
# Interpolación al centroide de la cuenca
stn.coords <- geom(project(stn, "epsg:9377"))[,c("x", "y")] %>% as.data.frame()
cnca.coords <- geom(centroids(cnca))[,c("x", "y")]
cnca.coords <- data.frame(x = cnca.coords[1], y = cnca.coords[2])
# Parámetro M de la cuenca de estudio
m <- idw(data.precip.anual, stn.coords, cnca.coords) %>% as.numeric()El cálculo de las IDF se va a realizar hasta las 3 horas de duración de la lluvia en intervalos de 15 minutos. Y se tendrán en cuenta los periodos de retorno de 2, 5, 10, 25, 50 y 100 años.
a <- 0.94; b <- 0.18; c <- 0.66; d <- 0.83 # Coeficientes
t <- seq(15, 180, 15) # Tiempos de lluvia
tr <- c(2, 5, 10, 25, 50, 100) # Periodos de retorno
idf <- NULL
for (i in tr){ # Calcula para cada periodo la intesidad y la guarda en idf
idfi <- a * i^b * m^d / (t/60) ^ c
idf <- rbind(idf, idfi)
}
idf <- cbind(tr, idf) %>% as.data.frame()
colnames(idf) <- c("Tr", t)
rownames(idf) <- as.character(tr)
knitr::kable(idf, digits = 1, caption = 'Curvas IDF')| Tr | 15 | 30 | 45 | 60 | 75 | 90 | 105 | 120 | 135 | 150 | 165 | 180 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 2 | 57.3 | 36.3 | 27.8 | 23.0 | 19.8 | 17.6 | 15.9 | 14.5 | 13.4 | 12.5 | 11.8 | 11.1 |
| 5 | 5 | 67.6 | 42.8 | 32.7 | 27.1 | 23.4 | 20.7 | 18.7 | 17.1 | 15.9 | 14.8 | 13.9 | 13.1 |
| 10 | 10 | 76.6 | 48.5 | 37.1 | 30.7 | 26.5 | 23.5 | 21.2 | 19.4 | 18.0 | 16.8 | 15.7 | 14.9 |
| 25 | 25 | 90.3 | 57.1 | 43.7 | 36.2 | 31.2 | 27.7 | 25.0 | 22.9 | 21.2 | 19.8 | 18.5 | 17.5 |
| 50 | 50 | 102.3 | 64.7 | 49.5 | 41.0 | 35.4 | 31.4 | 28.3 | 25.9 | 24.0 | 22.4 | 21.0 | 19.8 |
| 100 | 100 | 115.9 | 73.3 | 56.1 | 46.4 | 40.1 | 35.5 | 32.1 | 29.4 | 27.2 | 25.4 | 23.8 | 22.5 |
# Reorganiza para plotearlo
idf.tab <- reshape2::melt(idf, id.vars = "Tr", variable.name = "t")
idf.tab$Tr <- as.character(idf.tab$Tr)
idf.plot <- ggplot(idf.tab, aes(t, value, group = Tr, color = Tr)) +
geom_line() + theme_bw() +
labs(title = "Curvas IDF", subtitle = "Método Diaz-Granados & Vargas",
x = "Tiempo (min)", y = "I (mm/hr)") +
theme(plot.title = element_text(hjust = 0.5))
idf.plot
4 Estimación de caudales
Para el presente ejercicio se va a suponer la construcción de una alcantarilla de 0.9 m de diámetro, de acuerdo con el Manual de Drenaje para Carreteras del Invias (2011) para este tipo de obras el periodo de retorno es de 10 años.
4.1 Método Racional
En el método racional se aplica la siguiente ecuación
\[Q=CIA\] Donde, \(C\) es el coeficiente de escorrentía, \(I\) es la intensidad de la lluvia leída de las IDF en un tiempo igual al tiempo de concentración de la cuenca y \(A\) es el área de la cuenca. Las unidades de \(I\) y \(A\) debe ser consistentes de acuerdo a si se desean obtener los caudales en m³/s o en L/s.
Los tiempos de concentración son un tema abierto, para la presente aplicación se van a utilizar las ecuaciones de Kirpich, Témez, Ranser (SCS), V.T. Chow; las cuales se muestran en ese orden a continuación:
\[Tc=0.06628\left(\frac{L}{S_{0}^{0.5}}\right)^{0.77}\]
\[Tc=0.30\left(\frac{L}{(S_{0}*100)^{0.25}}\right)^{0.76}\]
\[Tc=0.947\left(\frac{L^{3}}{H}\right)^{0.385}\]
\[Tc=0.273\left(\frac{L}{S_{0}^{0.5}}\right)^{0.64}\]
Donde, \(Tc\) es el tiempo de concentración, en horas; \(L\) es la longitud del cauce principal, en kilómetros; \(S_{0}\) es la pendiente del cauce principal, en m/m; y \(H\) es la diferencia de cotas entre los extremos del cauce principal, en m.
l <- cnca$L / 1000 # Longitud del cauce principal
s <- cnca$S # Pendiente del cauce principal
h <- cnca$CMax - cnca$CMin # Diferencia de cotas en el cauce principal
# Tiempos de concentración en minutos
tc <- c(0.06628 * (l / s ^ 0.5) ^ 0.77,
0.30 * (l / (s * 100) ^ 0.25) ^ 0.79,
0.947 * (l ^ 3 / h) ^ 0.385,
0.273 * (l / s ^ 0.5) ^ 0.64) * 60
tc.tab <- data.frame("Ecuación" = c("Kirpich", "Témez", "Ranser", "V.T. Chow"),
Tiempo = tc)
knitr::kable(tc.tab, digits = 1, caption = "Tiempos de Concentración")| Ecuación | Tiempo |
|---|---|
| Kirpich | 16.6 |
| Témez | 20.5 |
| Ranser | 14.1 |
| V.T. Chow | 53.8 |
Con estos valores se debe escoger un tiempo de concentración, \(Tc\), para la cuenca. Existen varias posibilidades: 1) el menor \(Tc\) generará caudal mas grande y por consiguiente mayor nivel de seguridad en las obras, a riesgo de sobrestimar; 2) hacer un promedio de los \(Tc\); 3) analizar las condiciones bajo las cuales se propuso cada método y escoger el mas cercano a la cuenca de estudio; entre otros.
Una vez se haya escogido un \(Tc\) se puede calcular la intensidad con la ecuación de Diaz-Granados & Vargas, interpolarlo de las IDF ya calculadas o escoger el del tiempo más cercano. En el presente caso se va a escoger el valor de 15 minutos de las IDF que se acerca a los primeros 3 métodos calculados.
El coeficiente de escorrentía se obtiene ponderando el valor asociado a cada cobertura y tipo de suelo por el área que ocupa dentro de la cuenca.
suelos <- crop(suelos, cnca) # Recorte de las coberturas y suelos al área de la cuenca
coef <- sum(expanse(suelos) * suelos$Ce) / sum(expanse(suelos)) # Suma ponderada de los coeficientes
sueloSpatial <- as(project(suelos, "epsg:4326"), "Spatial")
leaflet() %>% addTiles() %>%
setView(lat = mean(cnca.ext[3:4]), lng = mean(cnca.ext[1:2]), zoom = 15) %>%
addPolygons(data = sueloSpatial, fillOpacity = 0.7, weight = 2, color = "#8B7D7B", fillColor = "#FFE4E1",
popup = ~sprintf("Cobertura: %s\nTextura: %s", nivel_2, Textura)) %>%
addPolygons(data = cncaSpatial, color = "#830000", fillOpacity = 0, popup = "Cuenca") sprintf("El coeficiente de escorrentía de la cuenca es de %.3f", coef)## [1] "El coeficiente de escorrentía de la cuenca es de 0.405"En este punto ya se tienen todos los elementos para aplicar el método racional.
q <- coef * (idf["10", "15"] / 3.6e6) * cnca$Area
sprintf("El caudal de la cuenca para un periodo de retorno de 10 años se estima en %.3f m³/s (%.f L/s)", q, q*1000)## [1] "El caudal de la cuenca para un periodo de retorno de 10 años se estima en 26.311 m³/s (26311 L/s)"4.2 Hidrograma unitario y convolución de caudales
Para el cálculo de caudales por el método del hidrograma unitario primero se debe estimar una tormenta de diseño, en el presente caso se hará mediante el bloque alterno con un tiempo de duración de la tormenta de 3 horas.
A partir de las IDF se puede obtener la precipitación total para esa duración, así como la forma en la que se distribuye en cada intervalo de 15 minutos.
# Precipitación total para cada duración
precip <- (idf["10",-1] * as.numeric(colnames(idf)[-1]) / 60) %>% as.numeric()
hiet <- c(precip[1], diff(precip)) # Precipitación por intervalo
df.plot <- rbind(data.frame(t = seq(15, 180, 15), v = precip, type = "Acumulada"),
data.frame(t = seq(15, 180, 15), v = hiet, type = "Intervalo"))
p <- ggplot(df.plot, aes(t, v, fill = type)) + theme_bw() +
geom_col(position = "dodge") +
scale_fill_manual(name = NULL, values = c("#00008B", "#8EE5EE")) +
labs(title = "Precipitación Acumulada y Producida en cada intervalo",
x = "Tiempo (min)", y = "mm") +
theme(plot.title = element_text(hjust = 0.5))
p
Posteriormente, se reorganiza de acuerdo con el bloque alterno. La mayor precipitación se ubica en la mitad y las demás se empiezan a ubicar a cada lado.
# Reorganización mediante bloque alterno - Hietograma
hiet <- hiet[c(seq(11, 1, -2), seq(2, 12, 2))]
df.plot <- data.frame(t = seq(15, 180, 15), v = hiet)
p <- ggplot(df.plot, aes(t, v)) + theme_bw() +
geom_col(fill = "#00008B") +
labs(title = "Hietograma de Precipitación Total",
x = "Tiempo (min)", y = "mm") +
theme(plot.title = element_text(hjust = 0.5))
p
Al hietograma de precipitación se aplican las perdidas por infiltración dadas por el número de curva, el cual se calcula de manera ponderada como se hizo con el coeficiente de escorrentía. A partir de este, se obtiene la abstracción inicial con \(I_{a}=0.2S\), donde \(S=\frac{25400-254CN}{CN}\) es la retención máxima de la cuenca, en mm. Con esto se puede obtener la precipitación efectiva (\(P_{e}\)), en mm; de acuerdo con:
\[P_{e}=\frac{(P-I_{a})^{2}}{P-I_{a}+S}\] O lo que es equivalente:
\[P_{e}=\frac{(P-0.2S)^{2}}{P+0.8S}\]
Donde, \(P\) es la precipitación total. Al ser precipitación total, se debe acumular el hietograma, estimar la precipitación efectiva en cada intervalo, teniendo en cuenta de que cuando la precipitación sea menor a \(0.2S\) se toma como nula. Finalmente. se obtiene el hietograma efectivo como las diferencias en la precipitación efectiva entre intervalos.
# Número de curva ponderado para la cuenca
cn <- sum(expanse(suelos) * suelos$CN) / sum(expanse(suelos))
# Precipitación acumulada
p.acc <- cumsum(hiet)
# Retención máxima en la cuenca
s <- 25400 / cn - 254
p.ef.acc <- (p.acc - 0.2 * s) ^ 2 / (p.acc + 0.8 * s)
p.ef.acc[p.acc < 0.2 * s] <- 0
hiet.ef <- c(p.ef.acc[1], diff(p.ef.acc))
df.plot$v <- hiet.ef
p <- ggplot(df.plot, aes(t, v)) + theme_bw() +
geom_col(fill = "#00868B") +
labs(title = "Hietograma de Precipitación Efectiva",
x = "Tiempo (min)", y = "mm") +
theme(plot.title = element_text(hjust = 0.5))
p
Una vez se tiene la precipitación efectiva, se procede con la aplicación del hidrograma unitario y la convolución. El hidrograma unitario que se va a usar es el del Soil Conservation Service (SCS). Dentro de este método se debe tener el tiempo de retraso, \(T_{lag}\) y el tiempo al pico, \(T_{p}\), en horas. En algunos casos se acepta que \(T_{lag}=0.6Tc\), sin embargo se suele recomendar la siguiente expresión:
\[T_{lag}=\frac{L_{p}^{0.8}(1000/CN-9)^{0.7}}{1900S_{w}^{0.5}}\] Donde, \(L_{p}\) es la longitud del cauce principal, en pies; y \(S_{w}\) es la pendiente media de la cuenca, en porcentaje. En el caso del tiempo al pico se calcula como \(T_{p}=\frac{D}{2}+t_{lag}\), donde, \(D\) es el intervalo de tiempo en el que se está trabajando, en horas; 0.25 hr (15 min) para este caso.
Adicionalmente, se requiere el caudal pico, \(U_{p}\), obtenido de la siguiente forma:
\[U_{p}=0.208\frac{A}{T_{p}}\]
Donde, \(A\) es el área de la cuenca, en km².
sw <- cnca$Sw
lp <- l * 3280.84
tlag <- (lp ^ 0.8 * (1000 / cn - 9) ^ 0.7 / 1900 / sw ^ 0.5)
tp <- (tlag + 0.25 / 2)
up <- 0.208 * cnca$Area / 1e6 / tpCon esto ya se puede hacer el cálculo del hidrograma unitario y la respectiva convolución de caudal. El hidrograma unitario es adimensional, las abscisas son \(t/t_{p}\) y las ordenadas \(q/U_{p}\). Se debe multiplicar por el \(t_{p}\) y el \(U_{p}\) de la cuenca en estudio. Luego, se multiplica por cada uno de los intervalos del hietograma de precipitación efectiva, así se obtiene la respuesta de cada uno de estos.
La convolución consiste agrupar todos los hidrogramas en una única respuesta de la cuenca. Para los tiempos tienen de todos tienen que ser consistentes. Lo cual se logra interpolandolos a una única escala de tiempo que tenga los mismo pasos de tiempo que el hietograma, es decir de 15 minutos. La última consideración es que en la respuesta de la cuenca, cada hidrograma inicia un paso de tiempo después que el anterior.
Con esto, ya se tienen todas las consideraciones para obtener el hidrograma de la cuenca.
uh <- read.csv("~/Diego/04.InformaciónGeneral/UH-SCS.csv")
# Conversión de las abscisas a los tiempos de la cuenca.
# Se hizo la conversión para trabajar en minutos.
uh[,1] <- uh[,1] * tp * 60
# Conversión de las ordenadas a los caudales de la cuenca, en m³/s.
uh[,2] <- uh[,2] * up
# Serie de tiempo para el hidrograma de la cuenca, 9 horas (540 min) cada 15 minutos.
t <- seq(15, 540, 15)
# Interpolación del hidrograma unitario a la serie de tiempo
uh <- data.frame(t = t, q = approx(uh[,1], uh[,2], t, rule = 2)[["y"]])
qu <- NULL # Para guardar todos los hidrogramas
for (i in 1:length(hiet.ef)){
qi <- uh[,2] * hiet.ef[i] # Hidrograma del paso de tiempo i
qi <- c(rep(0, (i-1)), qi) # Ajuste del tiempo de inicio del hidrograma
qu <- cbind(qu, qi[1:length(t)]) # Una columna por cada intervalo del hietograma.
}
colnames(qu) <- as.character(seq(15, 180, 15))
q <- apply(qu, 1, sum) # Convolución
# Para la gráfica
qu.tab <- qu[,apply(qu, 2, function(x) any(x > 0))]
qu.tab <- data.frame(t, qu) %>%
reshape2::melt(id.vars = "t", variable.name = "Dt")
df.plot <- data.frame(t = t, q = q)
p <- ggplot() + theme_bw() +
geom_line(data = qu.tab, aes(t, value, group = Dt), color = "#708090", size = 0.9) +
geom_line(data = df.plot, aes(t, q), color = "#27408B", size = 1.5) +
labs(title = "Hidrograma de la Cuenca", x = "Tiempo (min)", y = "Q (m³/s)") +
theme(plot.title = element_text(hjust = 0.5))
p
5 Replicabilidad
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] xts_0.12.2 zoo_1.8-11 ggplot2_3.3.6 phylin_2.0.2 leaflet_2.1.1
## [6] sf_1.0-8 raster_3.6-3 sp_1.5-0 terra_1.6-17 dplyr_1.0.10
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.9 lubridate_1.8.0 lattice_0.20-41 class_7.3-17
## [5] assertthat_0.2.1 digest_0.6.30 utf8_1.2.2 R6_2.5.1
## [9] plyr_1.8.7 evaluate_0.17 e1071_1.7-11 highr_0.9
## [13] pillar_1.8.1 rlang_1.0.6 rstudioapi_0.14 jquerylib_0.1.4
## [17] rmarkdown_2.17 labeling_0.4.2 rgdal_1.5-32 stringr_1.4.1
## [21] htmlwidgets_1.5.4 munsell_0.5.0 proxy_0.4-27 compiler_4.0.3
## [25] xfun_0.34 pkgconfig_2.0.3 rgeos_0.5-9 htmltools_0.5.3
## [29] tidyselect_1.2.0 tibble_3.1.8 bookdown_0.29 codetools_0.2-16
## [33] fansi_1.0.3 withr_2.5.0 grid_4.0.3 jsonlite_1.8.2
## [37] gtable_0.3.1 lifecycle_1.0.3 DBI_1.1.3 magrittr_2.0.3
## [41] units_0.8-0 scales_1.2.1 KernSmooth_2.23-17 cli_3.2.0
## [45] stringi_1.7.8 cachem_1.0.6 farver_2.1.1 reshape2_1.4.4
## [49] bslib_0.4.0 generics_0.1.3 vctrs_0.4.2 tools_4.0.3
## [53] glue_1.6.2 crosstalk_1.2.0 fastmap_1.1.0 yaml_2.3.6
## [57] colorspace_2.0-3 classInt_0.4-8 knitr_1.40 sass_0.4.2