1. Modelo Topmodel

Modelo lluvia escorrentía semidistribuido, que representa una versión simplificada de la variabilidad espacial de la respuesta hidrológica, en donde la topográfica se considera espacialmente. Este modelo usa un número reducido de parámetros, pero mantiene su interpretación física. Tiene una estructura simplificada, que reduce los requisitos de datos. Así, este modelo integra la capacidad de simular la distribución espacial de sus resultados en cualquier paso de tiempo, con una buena eficiencia computacional (Jeziorska & Niedzielski, 2018).

Los factores predominantes que influyen en la generación de caudales en el modelo son la topografía de la cuenca y las características del suelo (Franchini et al. 1996). La topografía se expresa cuantitativamente mediante el índice topográfico (también conocido como índice de humedad topográfica, TWI). Su valor se calcula a partir de la topografía; donde αi es el área que contribuye a la pendiente ascendente de la enésima cuenca y tan βi es la pendiente de la superficie del suelo de esta cuenca (Jeziorska & Niedzielski, 2018).

El área contribuyente aguas arriba (Figura 1) representa el área que potencialmente puede producir escorrentía hacia la salida de la cuenca. En la representación ráster del terreno, debe ser reemplazada por el área de drenaje pendiente arriba por unidad de longitud de contorno, que es equivalente al tamaño de celda de la cuadrícula DEM. Las áreas asociadas con valores altos de TWI tienden a saturarse primero y, por lo tanto, constituirán áreas contribuyentes potenciales del subsuelo o la superficie (Jeziorska & Niedzielski, 2018).

TWI se refiere al concepto de área de fuente variable de generación de escorrentía y se basa en las siguientes tres suposiciones simplificadas con respecto al sistema hidrológico:

  • la dinámica de la zona saturada se puede aproximar mediante sucesivas representaciones de estado estacionario.

  • el gradiente hidráulico de la zona saturada se puede aproximar por la pendiente topográfica de la superficie local (tan βi), y, por lo tanto, el nivel freático y el flujo saturado son paralelos a la pendiente de la superficie local.

  • la distribución de la transmisividad pendiente abajo con la profundidad es una función exponencial del déficit de almacenamiento o la profundidad del nivel freático.

Este enfoque implica que todos los puntos con el mismo valor de TWI responden de la misma manera. Los cálculos deben realizarse solo para valores representativos del índice, lo que simplifica enormemente el procedimiento y reduce el costo computacional al tiempo que mantiene la capacidad de identificación de los niveles freáticos y la humedad del suelo dentro de la cuenca (Jeziorska & Niedzielski, 2018).

Figura 1: Concepto básico del esquema TOPMODEL (Jeziorska & Niedzielski, 2018).

1.2. Procesos preliminares

Antes de iniciar, se realiza la limpieza del espacio de trabajo.

rm(list=ls(all=TRUE))

Se direcciona a la carpeta de trabajo, en donde se encuentra la información a usar y en donde se guardará la información generada.

setwd("D:/Maestria Hidrosistemas/3 Semestre/Modelacion de Procesos Hidrologicos/Parcial")

Para poder realizar este modelo se deben descargar e instalar las siguientes librerías:

  • topmodel.
  • Hmisc.

Estas deben ser instaladas usando el comando en la consola de install.packages(“Nombre de la librería”) Y deben ser llamadas en el código:

library("topmodel")
library("Hmisc")

1.3. Cálculo del índice topográfico

Topmodel requiere dos tipos de información topográfica, la distribución del índice topográfico (topidx), y el tiempo de concentración (delay). Al ser un modelo semidistribuido, no usa estos mapas directamente, si no usa los histogramas.

Para realizar el cálculo del índice topográfico, debe ingresarse el DEM de forma matricial, y se debe tener en cuenta la resolución del DEM, en este caso es 25. 
data(huagrahuma.dem)
DEM <- sinkfill(huagrahuma.dem, degree=0.1, res = 25) 
## number of iterations =  79 
## number of sinks left =  0
topindex <- topidx(DEM, resolution=25)
topindex = as.data.frame(topindex)

1.4. Cálculo del tiempo de concentración

Para este ejemplo el tiempo de concentración, fue generado con la función flowlength.

n <- 7 # un numero alto de clases daría un histograma suavizado.
delay <- flowlength(huagrahuma.dem)*25 
delay <- make.classes(delay, n)
delay <- delay[n:1,]
delay[,2] <- c(0, cumsum(delay[1:(n-1),2]))

2. Correr el modelo lluvia escorrentia

Inicialmente hay que cargar los datos que viene por defecto para probar este modelo, los datos son llamados huagrahuma y se realiza una exploración inicial de los datos, de entrada como lo son, índice topográfico, precipitación, evapotranspiración y los parámetros (qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,CD,dt) (Buytaert, 2018).
Parametro Definición Unidades
qs0 Flujo subsuperficial inicial por unidad de área (m)
lnTe Log del promedio de área del T0 (m2/h)
m Controla la tasa de disminución de la transmisividad
en el perfil del suelo
Sr0 Déficit de almacenamiento inicial en la zona radicular (m)
Sr max Déficit de almacenamiento máximo en zona radicular (m)
td Tiempo de retardo de la zona no saturada por unidad de (h/m)
déficit de almacenamiento
vr Caudal del canal dentro de la cuenca (m/h)
k0 Conductividad hidráulica superficial (m/h)
CD Acondicionamiento Capilar
dt Delta de tiempo (h) 
data(huagrahuma)
attach(huagrahuma)

str(huagrahuma)
## List of 6
##  $ ETp       : num [1:10000] 4.5e-06 4.5e-06 4.5e-06 4.5e-06 4.5e-06 4.5e-06 4.5e-06 4.5e-06 4.4e-06 4.4e-06 ...
##  $ Qobs      : num [1:10000] 3.34e-05 NA 3.45e-05 NA 3.34e-05 ...
##  $ delay     :'data.frame':  4 obs. of  2 variables:
##   ..$ delay1: num [1:4] 0 1000 3000 5000
##   ..$ delay2: num [1:4] 0 0.1 0.5 1
##  $ parameters: Named num [1:11] 3.17e-05 -5.99e-01 2.13e-02 2.63e-03 8.68e-01 ...
##   ..- attr(*, "names")= chr [1:11] "qs0" "lnTe" "m" "Sr0" ...
##  $ rain      : num [1:10000] 0 0 0 0 0 0 0 0 0 0 ...
##  $ topidx    :'data.frame':  16 obs. of  2 variables:
##   ..$ breaks: num [1:16] 16.2 15.3 14.3 13.4 12.4 ...
##   ..$ area  : num [1:16] 0.000619 0.000464 0.001857 0.004488 0.003869 ...
head(topidx)
##     breaks         area
## 1 16.21528 0.0006190034
## 2 15.27187 0.0004642526
## 3 14.32846 0.0018570102
## 4 13.38504 0.0044877747
## 5 12.44163 0.0038687713
## 6 11.49822 0.0095945528
head(parameters)
##           qs0          lnTe             m           Sr0         Srmax 
##  3.167914e-05 -5.990615e-01  2.129723e-02  2.626373e-03  8.683245e-01 
##            td 
##  2.850000e+00
head(rain)
## [1] 0 0 0 0 0 0
head(ETp)
## [1] 4.5e-06 4.5e-06 4.5e-06 4.5e-06 4.5e-06 4.5e-06
plot(rain, type="h")
Se realiza la corrida del modelo y se grafican los resultados, comparando los caudales simulados con los observados, haciendo uso de la librería topmodel. 
Qsim <- topmodel(parameters, topidx, delay, rain, ETp)
plot(Qsim, type="l", col="red")
points(Qobs)

Se evalúa el modelo con una métrica de desempeño (Nash - Sutcliffe).

NSeff(Qobs, Qsim)
## [1] 0.6046691

3. Análisis de Sensibilidad

Se realiza un análisis de sensibilidad, para saber que parámetros en el modelo son sensibles. En este caso solo se realizará para una sola variable m, la función runif utiliza valores aleatorios con distribución normal. 
parameters["m"] <- runif(1, min = 0, max = 0.1)
parameters["m"]
##          m 
## 0.08497109
Correr nuevamente el modelo, usando el nuevo valor de m, y se vuelve a realizar la métrica de desempeño, que con este valor da 0.3997425, y se grafica nuevamente, para verificar si es una buena simulación. 
Qsim <- topmodel(parameters, topidx, delay, rain, ETp)
NSeff(Qobs, Qsim)
## [1] 0.1620278
plot(Qsim, type="l", col="red")
points(Qobs)

Ahora los valores de todos los parámetros se van variar 100 veces, de manera aleatoria.

n <- 100

qs0 <- runif(n, min = 0.0001, max = 0.00025)
lnTe <- runif(n, min = -2, max = 3)
m <- runif(n, min = 0, max = 0.1)
Sr0 <- runif(n, min = 0, max = 0.2)
Srmax <- runif(n, min = 0, max = 0.1)
td <- runif(n, min = 0, max = 3)
vch <- runif(n, min = 100, max = 2500)
vr <- runif(n, min = 100, max = 2500)
k0 <- runif(n, min = 0, max = 10)
CD <- runif(n, min = 0, max = 5)
dt <- 0.25

parameters <- cbind(qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,CD,dt)
Después de tener 100 conjuntos de parámetros, con valores aleatorios, se genera una nueva corrida, evaluando el valor de Nash - Sutcliffe. La función topmodel recibe la función acepta una tabla con juegos de parámetros un juego de parámetros por filas. Como ejemplo se ve la visualización de la sensibilidad para la variable lnTe. 
NS <- topmodel(parameters, topidx, delay, rain, ETp, Qobs = Qobs)
max(NS)
## [1] 0.8493528
plot(lnTe, NS, ylim = c(0,1))

4. Análisis de Incertidumbre - GLUE

Ahora se elegirá el umbral de comportamiento y se eliminará el set de parámetros “malo.”

parameters <- parameters[NS > 0.3,]
NS <- NS[NS > 0.3]
Y se realiza una nueva corrida, sin los sets de parámetros malos y se visualiza los caudales simulados para el primer paso de tiempo. 
Qsim <- topmodel(parameters,topidx,delay,rain,ETp)
hist(Qsim[1,], xlab= "Qsim")
Se construye ponderaciones basadas en las medidas de desempeño, y se generan los límites de predicción mediante el cálculo del cuantil ponderado, para esto se necesita la libreria Hmisc. 
weights <- NS - 0.3
weights <- weights / sum(weights)


limits <- apply(Qsim, 1, "wtd.quantile", weights = weights,
                probs = c(0.05,0.95), normwt=T)

plot(limits[2,], type="l", ylab='Limits')
points(limits[1,], type="l")
points(Qobs, col="red")
Por último, se revisa cuántas medidas caen fuera de los límites de predicción, y se calcula el ancho de los límites de predicción. 
outside <- (Qobs > limits[2,]) | (Qobs < limits[1,])
summary(outside)
##    Mode   FALSE    TRUE    NA's 
## logical    5007    1765    3228
mean(limits[2,] - limits[1,]) / mean(Qobs, na.rm=T)
## [1] 0.7413749
Buytaert, W. (2018). Package ’topmodel’ - Implementation of the Hydrological Model TOPMODEL in R. Repository CRAN.
Jeziorska, J., & Niedzielski, T. (2018). Applicability of TOPMODEL in the mountainous catchments in the upper Nysa Kłodzka river basin (SW Poland). Acta Geophysica, 66(2), 203–222. https://doi.org/10.1007/S11600-018-0121-6/FIGURES/8