Descargar DEM - Modelación Hidrológica

Este cuaderno es una guía para facilitar el desarrollo de las actividades en clase del curso Modelación de procesos hidrológicos de la Facultad de Estudios Ambientales y Rurales de la Pontificia Universidad javeriana, y fue consolidado con la ayuda de Ginna Vergara, estudiante de la Maestría en Hidrosistemas en la PUJ.

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 Link.

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 Link..

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 [@Jeziorska2018].

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 [@Jeziorska2018].

Figura 1: Concepto básico del esquema TOPMODEL [@Jeziorska2018].

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("C://Users//maria.arenasb//OneDrive - Pontificia Universidad Javeriana//Javeriana//Modelacion de Procesos//Topmodel")

Taller 1:

1. Cambie la dirección de su carpeta de trabajo
2. En esta carpeta coloque los archivos del taller enviados

Para poder realizar este modelo se deben descargar e instalar las librerías relacionadas a continuación. 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:

#install.packages("raster")
#install.packages("sp")
#install.packages("magicaxis")
#install.packages("hydroGOF")
#install.packages("moments")
#install.packages("lattice")
#install.packages("fBasics")
#install.packages("ggplot2")

library(topmodel)
library(raster)
library(sp)
library(magicaxis)
library(hydroGOF)
library(moments)
library(lattice)
library(fBasics)
library(ggplot2)

1.3. Cálculo del Índice Topográfico - TWI

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.

Primero iniciamos con la carga del Modelo de Elevación Digital (DEM)

#2. CARGA DEL ARCHIVO RASTER 
DEM<-raster("DEM_30m.tif") #archivo entregado para la clase
plot(DEM) #usamos este comando para ver el DEM

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 30.

Matriz_DEM<-getValues(DEM,format="matrix")
TI<-topidx(Matriz_DEM,resolution=30.9) # Se ingresa el DEM, la resolución y como opción una matriz de ríos. Regresa el indice topografico de cada celda y el ?rea que contribuye

TI_class<-make.classes(TI$atb,30)

TI1<- topidx(Matriz_DEM,resolution=30.9)$atb
image(TI1)

1.4. Cálculo del tiempo de concentración

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

#Se calcula la distancia que existe entre cada celda del DEM y el punto de salida de la cuenca.
len.river<-flowlength(Matriz_DEM)

# 2.Agrupa las celdas que tienen igual distancia promedio hasta la salida de la cuenca
delay<-make.classes(len.river,30) # Clasificarlo en las mismas clases
delay[,1]<-delay[end(delay)[1]:1,1]*30.9
delay[,2]<-cumsum(delay[,2])

2. Correr el modelo lluvia escorrentia

Inicialmente hay que cargar los datos de entrada. Este modelo requiere datos de Precipitación, Evapotranspiración, Caudal y fecha.

# ---- Cargar Datos ----
DATOS <- read.csv("VMM_topmodel11.csv", header = T, sep=",") # Datos en m/d?a
ETO <- data.frame(ET0=DATOS$ET,FECHA=DATOS$Fecha)
PREC <- data.frame(PREC=DATOS$P,FECHA=DATOS$Fecha)
CAUDAL <- data.frame(Q=DATOS$Q,FECHA=DATOS$Fecha)# Datos en m3/s

# Convertir caudal a m/d?a

AREA <- 17105000000 # m2 AREA DE LA CUENCA EL PROFUNDO
CAUDAL$Q <- CAUDAL$Q*86400/AREA # m/d?a

En el siguiente apartado se realiza una exploración inicial de los datos y los parámetros (qs0,lnTe,m,Sr0,Srmax,td,vch,vr,k0,CD,dt) [@Buytaert2018].

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)

Una vez se han cargado los datos de entrada, se genera el vector de secuencia para las agrupaciones temporales que se requieran

# Generar vector de los d?as con la funci?n secuencia seq() by day
H.days<-seq(from=as.Date("01-01-2000",format="%d-%m-%Y"),
            to=as.Date("31-12-2012",format="%d-%m-%Y"),by="day")
# Generar vector de los a?os con la funci?n secuencia seq() by year
H.year<-seq(from=as.Date("01-01-2000",format="%d-%m-%Y"),
            to=as.Date("31-12-2012",format="%d-%m-%Y"),by="year")

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.

Ahora los valores de todos los parámetros se van variar 15 veces, de manera aleatoria.Para esto se utiliza la función uniforme. Los límites mínimos y máximos se definen según la literatura.

# ---- Calibraci?n tipo montecarlo ----
runs<-15
qs0_opt<-runif(runs, min=0, max=0.00015)  # Crea aleatoriamente un vector entre el valor min y max
lnTe_opt<-runif(runs, min=1, max=3)
m_opt<-runif(runs, min=0.15, max=0.22)
Sr0_opt<-runif(runs, min=0, max=0.1)
Srmax_opt<-runif(runs, min=0, max=0.1)
td_opt<-runif(runs, min=60, max=100)
vch_opt<-runif(runs, min = 1500, max = 2500)
vr_opt<-runif(runs, min=1000, max=1500)
k0_opt<-runif(runs, min=7, max=10)
CD_opt<-runif(runs, min=3, max=5)
dt<-24 # Time step siempre en horas
parameters<-cbind(qs0_opt,lnTe_opt,m_opt,Sr0_opt,Srmax_opt,td_opt,vch_opt,vr_opt,k0_opt,CD_opt,dt)

Después de tener 15 conjuntos de parámetros, con valores aleatorios, se define el periodo de calibración. En el protocolo de modelación esta establecido 1/4 del tiempo total para la calibración.

#Periodo de calibraci?n 01/01/2000-31/12/2012
Start.cal<-which(H.days=="2000-01-01"); End.Cal<-which(H.days=="2012-12-31")
Period.cal<-c(Start.cal:End.Cal)

Iniciamos las corridas del proceso de calibración

result<-matrix(0,length(Period.cal),runs)
M.NSE_opt<-c()
pb <- winProgressBar(title="Test Progress Bar", min=0, max=runs, width = 300)
for(i in 1:runs){
  result[,i]<- topmodel(parameters[i,], TI_class, delay, PREC$P[Period.cal], ETO$ET0[Period.cal])
  M.NSE_opt[i]<-NSeff(CAUDAL$Q[Period.cal], result[,i])
  bias_cal<-pbias(result[,i],CAUDAL$Q[Period.cal])
 
  Sys.sleep(0.1) # slow down the code for illustration purposes
  setWinProgressBar(pb, i, title=paste( round(i/runs*100, 0),
                                        "% Avance"))
}
close(pb)

## NULL

Ahora hacemos el análisis de sensibilidad comparando con el NSE y generamos tablas y gráficas

pos<-which.max(M.NSE_opt)
NSEop_opt<-M.NSE_opt[pos]; qs0_op<-qs0_opt[pos];lnTe_op<-lnTe_opt[pos];m_op<-m_opt[pos];Sr0_op<-Sr0_opt[pos];
Srmax_op<-Srmax_opt[pos];td_op<-td_opt[pos];vch_op<-vch_opt[pos];vr_op<-vr_opt[pos];k0_op<-k0_opt[pos];CD_op<-CD_opt[pos]

Qop_opt<-result[,pos]
Qop_opt1<-Qop_opt*AREA/86400
# write.csv(Qop_opt1, file="Qop_feb2.csv") #guarda los archivos en formato csv

#save(M.NSE,result,file = "TOPMODEL_R_Result.RData")

ind<-which.max(M.NSE_opt)
par(mfrow=c(5,2))
for(k in c(1:7,8:10)){
  par(mgp=c(1.6,0.3,0))
  par(mar=c(2.5,3,0.5,1.5)+0.1)
  plot(parameters[,k],M.NSE_opt,ylab="NSE",xlab=colnames(parameters)[k],xaxt="n",yaxt="n",bty="n",cex=0.8,
       ylim=c(-2,1))
  points(parameters[ind,k],M.NSE_opt[ind],col="red",pch=19,cex=1.4)
  magaxis(2,las=2,cex.axis=0.7,col="gray45",lty=16,lwd=0.3)
  magaxis(1,cex.axis=0.8,col="gray45",lty=1,lwd=0.3,usepar=TRUE)
  box(col="gray45",lwd=0.3)
}

Ahora hacemos la gráfica de los caudales observados vs los simulados por TopModel

par(mfrow=c(1,1))
par(mar=c(2.5,3,4.5,1.5)+0.1)
Caudal_m3s<-zoo(CAUDAL$Q*AREA/86400,H.days)
caudal1<- Caudal_m3s
Cau_m3s_sim<-zoo(Qop_opt*AREA/86400,H.days)

plot(caudal1, type="l", col="steelblue",
     main="CAUDALES SIMULADOS Y OBSERVADOS VMM", 
     xlab="AÑO", ylab="CAUDAL(m3/s)", cex.main=0.8)
lines(Cau_m3s_sim, type="l", col="red",cex.main=0.8)

Taller 2:

1. Aumnte el número de corridas y analice si esto cambia el comportamiento del modelo y del NSE.
2. Varíe el periodo de calibración y ajustelo a 1/4 del tiempo total del estudio.

Finalmente se hace la validación con el mejor resutlado de la calibración. Definimos cada parámetro con el valor óptimo de la calibración:

# ----  montecarlo con datos de validacion----
runs<-1
qs0_val<-qs0_op  # trabaja con el valor optimo de la calibracion
lnTe_val<-lnTe_op
m_val<-m_op
Sr0_val<-Sr0_op
Srmax_val<-Srmax_op
td_val<-td_op
vch_val<-vch_op
vr_val<-vr_op
k0_val<-k0_op
CD_val<-CD_op
dt<-24 # Time step siempre en horas
parameters_val<-cbind(qs0_val,lnTe_val,m_val,Sr0_val,Srmax_val,td_val,vch_val,vr_val,k0_val,CD_val,dt)

Iniciamos la corrida de la Validación definiendo la ventana de tiempo

#Periodo de validacion 01/01/2000-31/12/2012
Start.cal<-which(H.days=="2009-01-01"); End.Cal<-which(H.days=="2012-12-31")
Period.cal<-c(Start.cal:End.Cal)

Iniciamos las corridas del proceso de validación

result1<-matrix(0,length(Period.cal),runs)
M.NSE_val<-c()

for(i in 1:runs){
  result1[,i]<- topmodel(parameters_val[i,], TI_class,delay, PREC$P[Period.cal], ETO$ET0[Period.cal])
  M.NSE_val[i]<-NSeff(result1[,i],CAUDAL$Q[Period.cal])
 }

pos<-which.max(M.NSE_val)
NSEop_val<-M.NSE_val[pos]; qs0_op_val<-qs0_val[pos];lnTe_op_val<-lnTe_val[pos];m_op_val<-m_val[pos];Sr0_op_val<-Sr0_val[pos];
Srmax_op_val<-Srmax_val[pos];td_op_val<-td_val[pos];vch_op_val<-vch_val[pos];vr_op_val<-vr_val[pos];k0_op_val<-k0_val[pos];CD_op_val<-CD_val[pos]

#write.table(c(NSEop_val,qs0_op_val,lnTe_op_val,m_op_val,Sr0_op_val,Srmax_op_val,td_op_val,vch_op_val,vr_op_val,k0_op_val,CD_op_val),"parametros_validacion_5.txt")
Qop_val<-result[,pos]

#write.table(Qop_val*AREA/86400,"TQsim_vali_5.txt") # en m3/s

# se grafica el resultado
Caudal_m3s<-zoo(CAUDAL$Q*AREA/86400,H.days)
plot(Caudal_m3s, type="l", col="steelblue",
     main="CAUDALES SIMULADOS Y OBSERVADOS VMM", 
     xlab="AÑO", ylab="CAUDAL(m3/s)", cex.main=0.8)
lines(H.days, Qop_val*AREA/86400, type="l", col="red",cex.main=0.8)

#axis.Date(side = 1, at = H.year, format="%Y",labels = TRUE,
#          tck=-0.018,las=2,cex.axis=0.63,col="gray45",lwd=0.3)

Taller 3:

1. Varíe el periodo de validación y ajustelo a 3/4 del tiempo total del estudio. Taller
2. Repita el calculo del TWI usando al menos otro DEM (por ejemplo el Aster y el hydrosheds). Responda: nota diferencias en el proceso de computo del modelo (menor tiempo, mayor NSE)