Runoff Review
Victor Munoz S.
2024-11-04
1 Introduction
SRK Consulting was retained by xxxx to prepare a review of the future meteorological condition to be expected at the site xxxx.
SRK has already prepared a meteorological evaluation at the site including site analogue daily records for total precipitation and mean temperature.
source(here::here("scripts","hydro_support.r"))
#https://developers.google.com/maps/gmp-get-started\
# key_google<-rstudioapi::askForSecret("Google Key")
key_google<-"AIzaSyC2wTuP1tK60AxlCGpPSNetOvyZEgnO6EA"
Longitude= -68.978826
Latitude = -26.076496
#Get elevation from google maps
Elevation=GM_elev(Longitude,Latitude,key_google)
site.info<-data.frame(Location="Site",
Longitude=Longitude,
Latitude=Latitude,
Elevation=Elevation)
#Site in looped with variables
site.info_xy<-site.info %>%
dplyr::rename(x="Longitude",y="Latitude",z="Elevation") %>%
reshape2::melt(id.vars="Location")
pandoc.table(site.info)| Location | Longitude | Latitude | Elevation |
|---|---|---|---|
| Site | -68.98 | -26.08 | 3917 |
2 Input files
2.1 Flow records
loc_flow_files<-dir(here::here("data","csv"),pattern = "q_mm_day.csv",full.names = T,recursive = T)
loc_flow_files_splt<-str_split(loc_flow_files,pattern = "\\/",simplify = T)
#ID table with locations
loc_flow_files_tbl<-data.frame(id=loc_flow_files_splt[,ncol(loc_flow_files_splt)-1] %>%
str_split(,pattern = "_",simplify = T) %>% .[,3] %>% as.numeric(),
file=loc_flow_files)
#ID table with location
flow_idx <- read_excel(here::here("data/csv/Ubicación Estaciones.xlsx")) %>%
dplyr::filter(complete.cases(ID)) %>%
janitor::clean_names() %>%
left_join(.,loc_flow_files_tbl) %>%
mutate(dist_km=geosphere::distVincentyEllipsoid(
p1=cbind(Longitude,Latitude),
p2=cbind(longitud_grados,latitud_grados))/1000) %>%
arrange(dist_km)
pandoc.table(flow_idx %>% dplyr::select(-file),split.table=Inf)| id | nombre_estacion | latitud_grados | longitud_grados | elevacion_en_salida_m_s_n_m | fecha_inicio | fecha_fin | dist_km |
|---|---|---|---|---|---|---|---|
| 3022001 | Rio La Ola En Vertedero | -26.48 | -69.06 | 0 | 1986-06-17 | 2019-08-31 | 45.52 |
| 3414001 | Rio Pulido En Vertedero | -28.09 | -69.94 | NA | 1964-02-04 | 2020-04-07 | 242.2 |
| 2510001 | Rio San Pedro En Cuchabrachi | -22.82 | -68.2 | 2539 | 1947-06-23 | 2015-07-07 | 369 |
| 2105002 | Rio Salado En Sifon Ayquina | -22.29 | -68.34 | 2955 | 1975-03-23 | 2019-02-02 | 424.4 |
| 2103001 | Rio San Pedro En Parshall N"1 | -21.97 | -68.37 | 3764 | 1967-12-28 | 2018-05-15 | 459.5 |
| 2103014 | Rio Siloli Antes B.T. Fcab | -22.01 | -68.03 | 4315 | 2001-02-22 | 2020-04-07 | 460.7 |
| 2101001 | Rio Loa Antes Represa Lequena | -21.66 | -68.66 | 3293 | 1967-07-11 | 2020-04-07 | 490.7 |
flow_db<-lapply(seq_along(flow_idx$id),FUN=function(x){
read.csv(file = flow_idx$file[x]) %>%
mutate(id=flow_idx$id[x]) %>%
dplyr::rename("mm_d"=2) %>%
dplyr::filter(complete.cases(mm_d))
}) %>% rbindlist() %>%
mutate(date=as.Date(date))
flow_dly_z<-reshape2::dcast(data=flow_db,date~id,value.var = "mm_d") %>%
tk_zoo(silent=T)
flow_dly_z<-flow_dly_z[,as.character(flow_idx$id)]
plot(flow_dly_z,main="Daily flows [mm/d]")
daily2monthly.V2(flow_dly_z,FUN=sum) %>%
fortify.zoo() %>%
mutate(mon=lubridate::month(Index,label=T)) %>%
reshape2::melt(id.vars=c("Index","mon")) %>%
dplyr::rename(id="variable") %>%
mutate(id=as.numeric(as.character(id))) %>%
left_join(.,flow_idx %>% dplyr::select(id,nombre_estacion)) %>%
dplyr::filter(complete.cases(value)) %>%
mutate(nombre_estacion=factor(nombre_estacion,
labels = flow_idx$nombre_estacion,ordered = T)) %>%
ggplot(data=.,aes(x=mon,y=value))+
geom_boxplot()+
facet_wrap(~nombre_estacion,scales="free_y")+
theme_bw()+
labs(x="months",y="monthly runoff [mm/mon]",
title="Monthly average runoff")
left_join(flow_idx,
daily2annual.avg(flow_dly_z,FUN=sum) %>%
mutate(Station=as.numeric(Station)),by=c("id"="Station")) %>%
dplyr::select(id,nombre_estacion,month.abb,Ann) %>%
pandoc.table(caption="Total runoff mm/mon",split.table=Inf,round=2)| id | nombre_estacion | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 3022001 | Rio La Ola En Vertedero | 0.13 | 0.12 | 0.13 | 0.13 | 0.14 | 0.13 | 0.13 | 0.13 | 0.14 | 0.15 | 0.13 | 0.13 | 1.59 |
| 3414001 | Rio Pulido En Vertedero | 3.18 | 2.66 | 2.23 | 1.72 | 1.72 | 1.5 | 1.53 | 1.41 | 1.19 | 1.21 | 1.42 | 2.17 | 21.93 |
| 2510001 | Rio San Pedro En Cuchabrachi | 1.67 | 1.68 | 1.71 | 1.44 | 1.6 | 1.66 | 1.8 | 1.58 | 1.47 | 1.41 | 1.31 | 1.39 | 18.72 |
| 2105002 | Rio Salado En Sifon Ayquina | 2.32 | 2.64 | 2.04 | 1.46 | 1.63 | 1.53 | 1.58 | 1.5 | 1.33 | 1.26 | 1.11 | 1.15 | 19.56 |
| 2103001 | Rio San Pedro En Parshall N"1 | 2.45 | 2.21 | 2.46 | 2.42 | 2.48 | 2.36 | 2.36 | 2.28 | 2.28 | 2.44 | 2.48 | 2.42 | 28.65 |
| 2103014 | Rio Siloli Antes B.T. Fcab | 1.87 | 1.73 | 1.9 | 1.83 | 1.87 | 1.86 | 1.9 | 1.89 | 1.8 | 1.97 | 1.81 | 1.85 | 22.27 |
| 2101001 | Rio Loa Antes Represa Lequena | 0.83 | 1.09 | 0.77 | 0.65 | 0.71 | 0.71 | 0.75 | 0.75 | 0.75 | 0.76 | 0.63 | 0.61 | 9.01 |
2.2 Shapefile and Topography
library(sf)
library(terra)
loc_shp_files<-dir(here::here("data","csv"),pattern = "polygon.shp",
full.names = T,recursive = T)
loc_shp_files_splt<-str_split(loc_shp_files,pattern = "\\/",simplify = T)
#ID table with locations
loc_shp_files_tbl<-data.frame(id=loc_flow_files_splt[,ncol(loc_flow_files_splt)-1] %>%
str_split(,pattern = "_",simplify = T) %>% .[,3] %>% as.numeric(),
file=loc_shp_files)
wat_shp<-lapply(loc_shp_files_tbl$file,FUN=function(x){
wat_file_shp<-st_read(x)
}) %>% do.call(rbind,.)
#sorted out based on the distance form the site
wat_shp<-wat_shp[match(flow_idx$id,wat_shp$gauge_id),]# topo_ras<-rast(here::here("data","dem","SRTM_HD.bil"))
#
# plot(topo_ras)
# points(Longitude,Latitude)
#
#
# wat_bb<-ggmap::make_bbox(lon = as.numeric(ext(wat_shp)[1:2]),
# lat= as.numeric(ext(wat_shp)[3:4]),f = 0.1)
# #
# wat_ras<-rast(xmin=ext(wat_bb)[1],
# xmax=ext(wat_bb)[3],
# ymin=ext(wat_bb)[2],
# ymax=ext(wat_bb)[4])
#
# gtopo30_world_ras<-rast("//10.192.16.23/Library/Data/Gtopo30/World-GTOPO30.dem")
# #
# wat_gtopo30_ras<-crop(gtopo30_world_ras,wat_ras)
wat_gtopo30_ras<-rast(here::here("data/dem/wat_salares_norte_Rev2.tif"))
# dataframe
wat_ras_lst<-lapply(1:nrow(wat_shp),FUN=function(y){
x<-wat_shp[y,]
pander::pandoc.header(paste0(x$gauge_id," - ",x$gauge_name),
level=3)
wat_topo_ll<-mask(crop(wat_gtopo30_ras,x),x)
pander::pandoc.p('')
pander::pandoc.p('')
plot(wat_topo_ll)
pander::pandoc.p('')
pander::pandoc.p('')
wat_topo_m<-project(wat_topo_ll,"EPSG:24819")
watershed <- wat_topo_m>0
cell_size <- res(wat_topo_m)[1] * res(wat_topo_m)[2] # cell size in square units
watershed_area <- sum(values(watershed), na.rm=TRUE) * cell_size/1e6
# Get the coordinates of the watershed cells
coords <- xyFromCell(wat_topo_ll, which(values(watershed) == 1))
# Calculate the centroid
x_centroid <- mean(coords[,1], na.rm=TRUE)
y_centroid <- mean(coords[,2], na.rm=TRUE)
z_mean<-tidyterra::fortify( wat_topo_m)[,3] %>%
pull() %>% mean(.,na.rm=T)
z_min<-tidyterra::fortify( wat_topo_m)[,3] %>%
pull() %>% min(.,na.rm=T)
z_max<-tidyterra::fortify( wat_topo_m)[,3] %>%
pull() %>% max(.,na.rm=T)
tab_out<-tibble(id=x$gauge_id,
name=x$gauge_name,
area_topo_km2=watershed_area,
area_gis_km2=x$area_km2,
"diff_%"=100*(area_topo_km2/area_gis_km2-1),
x=x_centroid,y=y_centroid,z_min=z_min,z_mean=z_mean,
z_max=z_max)
tab_out%>%
pandoc.table(style="simple")
results<-list(ras=wat_topo_ll,tbl=tab_out)
return(results)
})2.2.1 3022001 - Rio La Ola En Vertedero
| id | name | area_topo_km2 | area_gis_km2 | diff_% |
|---|---|---|---|---|
| 3022001 | Rio La Ola En Vertedero | 12783 | 12311 | 3.832 |
| x | y | z_min | z_mean | z_max |
|---|---|---|---|---|
| -68.76 | -26.94 | 3543 | 4594 | 6676 |
2.2.2 3414001 - Rio Pulido En Vertedero
| id | name | area_topo_km2 | area_gis_km2 | diff_% |
|---|---|---|---|---|
| 3414001 | Rio Pulido En Vertedero | 2186 | 2022 | 8.106 |
| x | y | z_min | z_mean | z_max |
|---|---|---|---|---|
| -69.72 | -28.17 | 1379 | 3600 | 5669 |
2.2.3 2510001 - Rio San Pedro En Cuchabrachi
| id | name | area_topo_km2 | area_gis_km2 | diff_% |
|---|---|---|---|---|
| 2510001 | Rio San Pedro En Cuchabrachi | 1548 | 1416 | 9.321 |
| x | y | z_min | z_mean | z_max |
|---|---|---|---|---|
| -68.11 | -22.58 | 2565 | 4044 | 5698 |
2.2.4 2105002 - Rio Salado En Sifon Ayquina
| id | name | area_topo_km2 | area_gis_km2 | diff_% |
|---|---|---|---|---|
| 2105002 | Rio Salado En Sifon Ayquina | 899.6 | 806.8 | 11.51 |
| x | y | z_min | z_mean | z_max |
|---|---|---|---|---|
| -68.12 | -22.32 | 2989 | 4148 | 5524 |
2.2.5 2103001 - Rio San Pedro En Parshall N”1
| id | name | area_topo_km2 | area_gis_km2 |
|---|---|---|---|
| 2103001 | Rio San Pedro En Parshall N”1 | 1280 | 1157 |
| diff_% | x | y | z_min | z_mean | z_max |
|---|---|---|---|---|---|
| 10.65 | -68.13 | -22.01 | 3749 | 4482 | 5918 |
2.2.6 2103014 - Rio Siloli Antes B.T. Fcab
| id | name | area_topo_km2 | area_gis_km2 | diff_% |
|---|---|---|---|---|
| 2103014 | Rio Siloli Antes B.T. Fcab | 279.6 | 233.6 | 19.7 |
| x | y | z_min | z_mean | z_max |
|---|---|---|---|---|
| -67.94 | -21.99 | 4280 | 4692 | 5590 |
2.2.7 2101001 - Rio Loa Antes Represa Lequena
| id | name | area_topo_km2 | area_gis_km2 |
|---|---|---|---|
| 2101001 | Rio Loa Antes Represa Lequena | 2170 | 2053 |
| diff_% | x | y | z_min | z_mean | z_max |
|---|---|---|---|---|---|
| 5.703 | -68.65 | -21.33 | 3270 | 4090 | 5968 |
names(wat_ras_lst)<-loc_shp_files_tbl$id
#combine all the rasters
wat_ras_all<-lapply(wat_ras_lst,FUN=function(x){x$ras}) %>%
do.call(c,.)
#produce the index table
wat_ras_idx<-lapply(wat_ras_lst,FUN=function(x){x$tbl}) %>%
do.call(rbind,.)
#ordered based on distance
wat_ras_idx<-wat_ras_idx[match(flow_idx$id,wat_ras_idx$id),]
pandoc.table(wat_ras_idx,split.table=Inf)| id | name | area_topo_km2 | area_gis_km2 | diff_% | x | y | z_min | z_mean | z_max |
|---|---|---|---|---|---|---|---|---|---|
| 3022001 | Rio La Ola En Vertedero | 12783 | 12311 | 3.832 | -68.76 | -26.94 | 3543 | 4594 | 6676 |
| 3414001 | Rio Pulido En Vertedero | 2186 | 2022 | 8.106 | -69.72 | -28.17 | 1379 | 3600 | 5669 |
| 2510001 | Rio San Pedro En Cuchabrachi | 1548 | 1416 | 9.321 | -68.11 | -22.58 | 2565 | 4044 | 5698 |
| 2105002 | Rio Salado En Sifon Ayquina | 899.6 | 806.8 | 11.51 | -68.12 | -22.32 | 2989 | 4148 | 5524 |
| 2103001 | Rio San Pedro En Parshall N”1 | 1280 | 1157 | 10.65 | -68.13 | -22.01 | 3749 | 4482 | 5918 |
| 2103014 | Rio Siloli Antes B.T. Fcab | 279.6 | 233.6 | 19.7 | -67.94 | -21.99 | 4280 | 4692 | 5590 |
| 2101001 | Rio Loa Antes Represa Lequena | 2170 | 2053 | 5.703 | -68.65 | -21.33 | 3270 | 4090 | 5968 |
2.3 Watersheds regional context
library(rnaturalearth)
chile <- ne_countries(country = "Chile", returnclass = "sf",scale = 10)
ggplot()+
tidyterra::geom_spatraster(data=wat_gtopo30_ras)+
geom_sf(data=chile,fill="transparent",lwd=2)+
geom_sf(data=wat_shp,fill="lightgrey",alpha=0.8,col="purple",lwd=1)+
geom_point(aes(x=Longitude,y=Latitude),col="red",size=5)+
geom_label_repel(data=wat_ras_idx,aes(x=x,y=y,label=str_wrap(name,12)),
size=2.5,force=150)+
scale_fill_binned(type="viridis",breaks=seq(0,8000,2000))+
scale_y_continuous(limits=c(-29,-20.5))+
scale_x_continuous(limits=c(-71.6,-66.8))+
theme_bw()+
labs(fill="elevation\nmasl")
2.4 Air Temperature
Temperature was implemented based on the baseline report. Typically should be applicable in range of 3400 to 4600 masl
#read temperature
site_temp_tbl<-arrow::read_parquet(file = here::here("data/rds",
"temp_tdew_rh_dly_rev3.par"))
#Expression to change the temperature based on elevation
#expression tested!
#by default defined at Planta station=4544.743 masl
temp_adj_z<-function(z,ts=site_planta_t2m_dly_z,z0=4544.743){
reg_station_slope<-data.frame(variable=c("tavg","tmax","tmin"),
slope=c(-0.005806119,-0.011188779,0.002942172))
adj_tbl<-reg_station_slope %>%
mutate(elev_dif=z0-z) %>%
mutate(adj=-slope*elev_dif)
temp_z<-merge(
tmax=ts$tmax,
tavg=ts$tavg,
tmin=ts$tmin
) %>%
fortify.zoo() %>%
reshape2::melt(id.var="Index") %>%
left_join(adj_tbl) %>%
mutate(value_adj=value+adj) %>%
dplyr::select(Index,variable,value_adj) %>%
reshape2::dcast(data=.,formula=Index~variable,value.var = "value_adj") %>%
tk_zoo(silent=T)
return(temp_z)
}
#adjustment based on Projecto information
temp_adj_z2<-function(z){
ts_temp<-tk_zoo(site_temp_tbl[,1:4])
names(ts_temp)<-str_remove_all(names(ts_temp),"2m_")
return(temp_adj_z(z=z,ts = ts_temp,z0 = 4323.77))
}2.5 Total Precipitation
era5_land_ras<-terra::rast(here::here("data","netcdf","era5_land_prcp_dly.nc"))
#loading reanalysis
nc2z_fn<-function (Longitude, Latitude, source) {
nc_file <- case_when(source == "ERA5-LAND" ~ "era5_land_ras",
source == "ERA5" ~ "era5_ras", source == "MERRA2" ~ "merra2_ras",
source == "CHIRPS" ~ "chirps_ras", source == "GPM" ~
"gpm_ras", source == "WRF" ~ "wrf_ras",
source ==
"SM2RAIN" ~ "sm2rain_ras")
nc_ras <- eval(parse(text = nc_file))
site_vec <- terra::vect(data.frame(Longitude, Latitude),
geom = c("Longitude", "Latitude"), crs = "epsg:4326")
extract_ts <- terra::extract(x = nc_ras, y = site_vec, method = "simple",
raw = T)
reanalysis_z <- zoo::zoo(unlist(extract_ts[1, -1]), terra::time(nc_ras))
return(reanalysis_z)
}
#reanalysis at any point, where temperature define rainfall / snowfall
nc2z_elev_fn<-function(Longitude=site.info$Longitude,
Latitude=site.info$Latitude,
source="ERA5-LAND",elev=3500){
prcp_z<-nc2z_fn(Longitude = site.info$Longitude,
Latitude = site.info$Latitude,
source = source)
zoo::time(prcp_z)<-as.Date(time(prcp_z))
tavg_z<-temp_adj_z2(elev)
meteo_out_z<-merge(prcp=prcp_z,tavg_z) %>%
window(start="1980-01-01") %>%
fortify.zoo() %>%
mutate(rainfall=if_else(tavg>0,prcp,if_else(is.na(prcp),NA,0)),
snowfall=if_else(tavg<=0,prcp,if_else(is.na(prcp),NA,0))) %>%
tk_zoo(silent=T)
return(meteo_out_z)
}2.5.1 Rainfall
Values based on range at 3500 to 4800
2.6 Meteo Compiled
applied to site conditions
x1<-site.info$Longitude
y1<-site.info$Latitude
z1<-site.info$Elevation
meteo_info_fn<-function(x1,y1,z1){
#Set air temperature
meteo_info_z<-nc2z_elev_fn(Longitude =x1,
Latitude = y1,
elev=z1,
source = "ERA5-LAND")
meteo_info_z$tmin[(meteo_info_z$tmax<meteo_info_z$tmin)]<-NA
meteo_info_z<-na.approx(meteo_info_z)
meteo_info_z$tmax[(meteo_info_z$tmax<meteo_info_z$tmin)]<-NA
pe_z<-zoo(airGR::PE_Oudin(JD = lubridate::yday(meteo_info_z$tavg),
Temp = coredata(meteo_info_z$tavg),
Lat = y1,LatUnit = "deg"),
time(meteo_info_z))
meteo_info_z<-na.approx(meteo_info_z) %>%
fortify.zoo() %>%
mutate(rainfall=rainfall*rf_fn(z1),
snowfall=snowfall*sf_fn(z1),
prcp=rainfall+snowfall) %>%
tk_zoo(silent=T) %>%
merge(.,pe=pe_z)
meteo_info_z$prcp[meteo_info_z$prcp<0]<-0
return(meteo_info_z)
}
meteo_info_z<-meteo_info_fn(x1,y1,z1)
dwi.graph(meteo_info_z)
rbind(
daily2annual.avg(meteo_info_z[,c("tmax","tavg","tmin")],
FUN=mean),
daily2annual.avg(meteo_info_z[,c("rainfall","snowfall","prcp")],
FUN=sum)
) %>%
pandoc.table(round=1,split.table=Inf,
caption="site meteorology - mean monthly values 1980 to 2021")| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| tmax | 18 | 17.3 | 16.2 | 13.2 | 10.1 | 8.3 | 8 | 9.4 | 10.9 | 13.4 | 15.5 | 17.2 | 13.1 |
| tavg | 7.2 | 6.6 | 5.5 | 3 | 0.4 | -1.6 | -2.1 | -0.6 | 0.7 | 2.9 | 4.7 | 6.4 | 2.8 |
| tmin | -4.9 | -5.1 | -6 | -7.7 | -9.8 | -11.9 | -12.7 | -11.3 | -10.2 | -8.6 | -7.6 | -6 | -8.5 |
| rainfall | 5.1 | 4.4 | 2.8 | 0.6 | 0.6 | 0.1 | 0 | 0.2 | 0.4 | 0.5 | 0.6 | 1.9 | 17.2 |
| snowfall | 0 | 0 | 0.2 | 3 | 10.8 | 14.3 | 13.4 | 9.6 | 4.8 | 0.7 | 0.2 | 0.2 | 57.1 |
| prcp | 5.1 | 4.4 | 2.9 | 3.6 | 11.3 | 14.4 | 13.5 | 9.8 | 5.2 | 1.1 | 0.7 | 2.1 | 74.3 |
3 GR5J - Rio Pulido en Vertedero
3.1 Hypsometric Curve
dem_ras<-wat_ras_all[[as.character(id_sel$id)]]
dem_tbl<-tidyterra::fortify(dem_ras) %>%
dplyr::rename("z"=3) %>%
dplyr::filter(complete.cases(z))
plot(dem_ras)
x_ecdf<- dem_tbl$z %>% ecdf()
z_val<-seq(min(dem_tbl$z,na.rm = T),
max(dem_tbl$z,na.rm = T),
length.out = 1000)
z_hypso<-approx(x = x_ecdf(z_val),
y=z_val,
xout=seq(0.001,1,length.out = 101))
z_hypso$y[1]<-min(dem_tbl$z,na.rm = T)
plot(z_hypso,main="Hypsometric curve",xlab="percentage of area [%]",
ylab="Elevation[masl]")
3.2 Calibration 1983 - 1998
3.2.1 Meteorology & Flows
meteo_info_z<-meteo_info_fn(x = id_sel$x,
y = id_sel$y,
z = id_sel$z_mean)
plot(meteo_info_z,main="available information")
#r1 works between 1997-06 to 2008-12
#r2 1990-06 / 2008-12
start.cal.date<-as.Date("1980-06-01")
end.cal.date<-as.Date("1998-05-31")
#flow start 3 years after, that helps on the calibration
flow_start<-start.cal.date+years(3)
flow_end<-end.cal.date
Meteo.site.cal<-merge(meteo_info_z,
local=flow_dly_z[,as.character(id_sel$id)] %>%
window(start=flow_start,end=flow_end)) %>%
window(start=start.cal.date,end=end.cal.date)
plot(Meteo.site.cal)
3.2.2 Kienzle
T_thresh = 1.7
T_range = 14.3
# Tav <- (Tmax_C+Tmin_C)/2
p_rain_fn<-function(Tav){
P_rain <- ifelse(Tav<=T_thresh,
5*((Tav-T_thresh)/(1.4*T_range))^3+
6.76*((Tav-T_thresh)/(1.4*T_range))^2+
3.19*((Tav-T_thresh)/(1.4*T_range))+0.5,
5*((Tav-T_thresh)/(1.4*T_range))^3-
6.76*((Tav-T_thresh)/(1.4*T_range))^2+
3.19*((Tav-T_thresh)/(1.4*T_range))+0.5)
P_rain[P_rain>1] <- 1
P_rain[P_rain<0] <- 0
return(P_rain)
}3.2.3 Runoff Model
library(airGR)
InputsModel.cal<-CreateInputsModel(
FUN_MOD = RunModel_CemaNeigeGR5J,
DatesR = time(Meteo.site.cal) %>% as.POSIXlt(),
Precip = Meteo.site.cal$prcp%>% coredata,
PotEvap = Meteo.site.cal$pe %>% coredata,
TempMax = Meteo.site.cal$tmax%>% coredata,
TempMin = Meteo.site.cal$tmin%>% coredata,
TempMean = Meteo.site.cal$tavg%>%coredata,
ZInputs = id_sel$z_mean,
PrecipScale = T,
HypsoData = na.omit(z_hypso$y),
NLayers=4)
#capturing at different elvations
meteo_info_lst<-lapply(InputsModel.cal$ZLayers,FUN=function(x){
meteo_info_fn(x = id_sel$x,
y = id_sel$y,
z = x) %>%
window(start=start.cal.date,end=end.cal.date)
})
#prcp
prcp_lst<-lapply(meteo_info_lst,FUN=function(x){x$prcp %>% coredata()})
names(prcp_lst)<-paste0("L",1:length(prcp_lst))
InputsModel.cal$LayerPrecip<-prcp_lst
#air temperature
tavg_lst<-lapply(meteo_info_lst,FUN=function(x){x$tavg %>% coredata()})
names(tavg_lst)<-paste0("L",1:length(tavg_lst))
InputsModel.cal$LayerTempMean<-tavg_lst
#Snow Fraction as Kienzle
sf_lst<-lapply(meteo_info_lst,FUN=function(x){
# if_else(x$prcp<=0,0,coredata(x$snowfall/x$prcp))
if_else(x$prcp<=0,0, 1-p_rain_fn(Tav = x$tavg) %>%
coredata())
})
names(sf_lst)<-paste0("L",1:length(sf_lst))
InputsModel.cal$LayerFracSolidPrecip<-sf_lst
Ind_Run.cal <- seq(grep(start.cal.date, time(Meteo.site.cal)),
grep(end.cal.date, time(Meteo.site.cal)))
RunOptions.cal <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR5J,
InputsModel = InputsModel.cal,
IndPeriod_Run = Ind_Run.cal,
IniStates = NULL, IniResLevels = NULL,
IndPeriod_WarmUp =1:60,
IsHyst = F )
###CALIBRATION START##
InputsCrit.cal <- CreateInputsCrit(
FUN_CRIT =list(ErrorCrit_KGE,
ErrorCrit_KGE),
InputsModel = InputsModel.cal,
RunOptions = RunOptions.cal,
VarObs=list("Q", "Q"),
Obs = list(coredata(Meteo.site.cal$local),
coredata(Meteo.site.cal$local)),
transfo = list("0.1", "1.2"),
Weights = list(0.4, 0.6)
)
#limit for x4 extended beyond 20 as default
search_ranges<-matrix(
nrow=2,ncol=7,
c(0,-10903,0,0.5,0,0,0,
21807,10903,21807,30,1,1,109.04),byrow = T)
CalibOptions.cal <- CreateCalibOptions(
SearchRanges = search_ranges,
FUN_MOD = RunModel_CemaNeigeGR5J,
FUN_CALIB = Calibration_Michel,
StartParamList = matrix(
nrow=1,ncol=7,
c(0.0001,-0.110,238,20,0.345,0.996,0.724)),
# SearchRanges =search_ranges,
IsHyst = F)
OutputsCalib.cal <- Calibration_Michel(InputsModel = InputsModel.cal,
RunOptions = RunOptions.cal,
InputsCrit = InputsCrit.cal,
CalibOptions = CalibOptions.cal,
FUN_MOD = RunModel_CemaNeigeGR5J)
Param_eval <- OutputsCalib.cal$ParamFinalR
OutputsModel.cal <- RunModel_CemaNeigeGR5J(InputsModel = InputsModel.cal,
RunOptions = RunOptions.cal,
Param = Param_eval)
NSE.cal<-hydroGOF::NSE(sim=OutputsModel.cal$Qsim,
obs=Meteo.site.cal$local %>% coredata)
KGE.cal<-hydroGOF::KGE(sim=OutputsModel.cal$Qsim,
obs=Meteo.site.cal$local %>% coredata)
plot.OutputsModel(OutputsModel.cal,
which=c(#"Precip", "Temp",
"SnowPack", "Flows", #"Error",
"Regime",
"CumFreq", "CorQQ"),
Qobs = Meteo.site.cal$local%>%
window(start=start.cal.date,
end=end.cal.date) %>%
coredata)
data.frame(NSE=NSE.cal,KGE=KGE.cal) %>%
pandoc.table(style="simple",caption="Runoff Model Performance - Daily")| NSE | KGE |
|---|---|
| 0.456 | 0.7289 |
3.2.4 Results
met_gl<-data.frame(
Date=OutputsModel.cal$DatesR %>% as.Date(),
P=OutputsModel.cal$Precip,
P_eff=(0.25*OutputsModel.cal$CemaNeigeLayers$Layer01$PliqAndMelt+
0.25*OutputsModel.cal$CemaNeigeLayers$Layer02$PliqAndMelt+
0.25*OutputsModel.cal$CemaNeigeLayers$Layer03$PliqAndMelt+
0.25*OutputsModel.cal$CemaNeigeLayers$Layer04$PliqAndMelt),
E=OutputsModel.cal$PotEvap
) %>%
mutate(
# Pn=if_else(P>=E,P-E,0),
# En=if_else(P>=E,0,E-P),
`P net(Pn)`=OutputsModel.cal$Pn,
`P storage(Ps)`=OutputsModel.cal$Ps,
`Pn-Ps`=`P net(Pn)`-`P storage(Ps)`,
# Prod=OutputsModel1$Prod,
# AE=OutputsModel1$AE,
`aft.Perc.rout.(Pr)`=OutputsModel.cal$PR,
`Exchange (F)`=OutputsModel.cal$AExch,
Q=OutputsModel.cal$Qsim
)
wb_gl.z<-tk_zoo(met_gl,silent=T)
wb_gl.syn<-merge(wb_gl.z,
Q_obs=Meteo.site.cal$local%>%
window(start=start.cal.date,
end=end.cal.date)) %>% daily2annual.avg(FUN=sum) %>%
data.frame()
pandoc.table(wb_gl.syn,
round=2,split.table=Inf,missing="-",
caption="Watershed Summary from 1980 to 2020") | Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| P | 3.76 | 3.47 | 2.13 | 1.14 | 9.15 | 17.12 | 12.34 | 11.02 | 5.29 | 1.08 | 0.8 | 2.41 | 69.7 |
| P_eff | 11.53 | 5.53 | 3.24 | 0.82 | 1.44 | 1.66 | 1.44 | 3.36 | 6.6 | 18.47 | 13.42 | 13.62 | 81.12 |
| E | 77.58 | 61.93 | 55.08 | 34.9 | 21.28 | 12.65 | 12.28 | 20.03 | 28.84 | 46.92 | 60.71 | 74.85 | 507.1 |
| P net(Pn) | 1.94 | 0.09 | 0.05 | 0.08 | 0.18 | 0.49 | 0.69 | 0.83 | 2.72 | 8.93 | 5.26 | 3.86 | 25.13 |
| P storage(Ps) | 0 | 0 | 0 | 0 | 0.01 | 0.01 | 0 | 0.01 | 0.01 | 0.01 | 0.01 | 0.01 | 0.07 |
| Pn-Ps | 1.94 | 0.09 | 0.05 | 0.08 | 0.18 | 0.48 | 0.69 | 0.82 | 2.72 | 8.92 | 5.26 | 3.85 | 25.06 |
| aft.Perc.rout.(Pr) | 1.94 | 0.09 | 0.05 | 0.08 | 0.18 | 0.48 | 0.69 | 0.82 | 2.72 | 8.92 | 5.26 | 3.85 | 25.06 |
| Exchange (F) | 0.33 | 0.3 | 0.35 | 0.35 | 0.37 | 0.36 | 0.38 | 0.38 | 0.37 | 0.38 | 0.34 | 0.33 | 4.22 |
| Q | 3.88 | 3.17 | 2.61 | 2.09 | 1.87 | 1.93 | 1.76 | 1.63 | 1.57 | 2.09 | 3.26 | 3.99 | 29.86 |
| Q_obs | 4.7 | 3.11 | 2.74 | 2.08 | 2.18 | 1.8 | 1.9 | 1.92 | 1.69 | 1.9 | 2.35 | 3.2 | 29.56 |
par_gl.tbl<-cbind(
tribble(~Parameter, ~definition,
"GR4J X1","production store capacity [mm]",
"GR4J X2","intercatchment exchange coefficient [mm/d]",
"GR4J X3","routing store capacity [mm]",
"GR4J X4","unit hydrograph time constant [d]",
"GR5J X5","intercatchment exchange threshold [-]",
"CemaNeige X1","weighting coefficient for snow pack thermal state [-]",
"CemaNeige X2","degree-day melt coefficient [mm/°C/d]"
),value=Param_eval)
pandoc.table(par_gl.tbl,round=2,caption="GR4J with snowmelt parameters")| Parameter | definition | value |
|---|---|---|
| GR4J X1 | production store capacity [mm] | 0 |
| GR4J X2 | intercatchment exchange coefficient [mm/d] | -0.03 |
| GR4J X3 | routing store capacity [mm] | 265.1 |
| GR4J X4 | unit hydrograph time constant [d] | 23.67 |
| GR5J X5 | intercatchment exchange threshold [-] | 0.44 |
| CemaNeige X1 | weighting coefficient for snow pack thermal state [-] | 0.98 |
| CemaNeige X2 | degree-day melt coefficient [mm/°C/d] | 6.47 |
g10<-wb_gl.syn|>
mutate(Station=paste0(stringr::str_pad(1:nrow(wb_gl.syn),width=2,
side = "left",pad = 0),". ",Station)) %>%
# window(start=as.Date("1981-01-01"),end=as.Date("2020-12-31")) |>
# daily2annual.avg() |>
reshape2::melt(id.vars="Station") |>
dplyr::filter(variable!="Ann") |>
ggplot(aes(x=variable,y=value))+
#geom_line()+
geom_bar(stat = "identity",position = "dodge",aes(fill=Station),col="black")+
#geom_boxplot()+
geom_text(aes(label=round(value,1)),alpha=0.9,size=3,nudge_y = 4)+
scale_fill_viridis_d(option="E")+
#scale_y_continuous(limits=c(0,280))+
facet_grid(Station~.,scale="free_y")+
theme_bw()+
# scale_y_continuous(limits=c(0,220),breaks = seq(0,250,50))+
theme(legend.position = "none",
strip.text.y = element_text(angle = 0))+
labs(x="",y="Mean monthly values [mm/mon]",
title="Mean monthly values")
print(g10)
3.3 Validation 1998 - 2012
On the validation, the model is tested with another period of the time series.
#r1 works between 1997-06 to 2008-12
#r2 1990-06 / 2008-12
start.val.date<-as.Date("1995-06-01")
end.val.date<-as.Date("2010-06-30")
#flow start 3 years after, that helps on the calibration
flow_start<-start.val.date+years(3)
flow_end<-end.val.date
Meteo.site.val<-merge(meteo_info_z,
local=flow_dly_z[,as.character(id_sel$id)] %>%
window(start=flow_start,end=flow_end)) %>%
window(start=start.val.date,end=end.val.date)
plot(Meteo.site.val)
3.3.1 Runoff Model
library(airGR)
InputsModel.val<-CreateInputsModel(
FUN_MOD = RunModel_CemaNeigeGR5J,
DatesR = time(Meteo.site.val) %>% as.POSIXlt(),
Precip = Meteo.site.val$prcp%>% coredata,
PotEvap = Meteo.site.val$pe %>% coredata,
TempMax = Meteo.site.val$tmax%>% coredata,
TempMin = Meteo.site.val$tmin%>% coredata,
TempMean = Meteo.site.val$tavg%>%coredata,
ZInputs = id_sel$z_mean,
PrecipScale = T,
HypsoData = na.omit(z_hypso$y),
NLayers=4)
#capturing at different elvations
meteo_info_lst<-lapply(InputsModel.val$ZLayers,FUN=function(x){
meteo_info_fn(x = id_sel$x,
y = id_sel$y,
z = x) %>%
window(start=start.val.date,end=end.val.date)
})
#prcp
prcp_lst<-lapply(meteo_info_lst,FUN=function(x){x$prcp %>% coredata()})
names(prcp_lst)<-paste0("L",1:length(prcp_lst))
InputsModel.val$LayerPrecip<-prcp_lst
#air temperature
tavg_lst<-lapply(meteo_info_lst,FUN=function(x){x$tavg %>% coredata()})
names(tavg_lst)<-paste0("L",1:length(tavg_lst))
InputsModel.val$LayerTempMean<-tavg_lst
#Snow Fraction as Kienzle
sf_lst<-lapply(meteo_info_lst,FUN=function(x){
# if_else(x$prcp<=0,0,coredata(x$snowfall/x$prcp))
if_else(x$prcp<=0,0, 1-p_rain_fn(Tav = x$tavg) %>%
coredata())
})
names(sf_lst)<-paste0("L",1:length(sf_lst))
InputsModel.val$LayerFracSolidPrecip<-sf_lst
Ind_Run.val <- seq(grep(start.val.date, time(Meteo.site.val)),
grep(end.val.date, time(Meteo.site.val)))
RunOptions.val <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR5J,
InputsModel = InputsModel.val,
IndPeriod_Run = Ind_Run.val,
IniStates = NULL, IniResLevels = NULL,
IndPeriod_WarmUp =1:60,
IsHyst = F )
OutputsModel.val <- RunModel_CemaNeigeGR5J(InputsModel = InputsModel.val,
RunOptions = RunOptions.val,
Param = Param_eval)
NSE.val<-hydroGOF::NSE(sim=OutputsModel.val$Qsim,
obs=Meteo.site.val$local %>% coredata)
KGE.val<-hydroGOF::KGE(sim=OutputsModel.val$Qsim,
obs=Meteo.site.val$local %>% coredata)
plot.OutputsModel(OutputsModel.val,
which=c(#"Precip", "Temp",
"SnowPack", "Flows", #"Error",
"Regime",
"CumFreq", "CorQQ"),log_scale=FALSE,
Qobs = Meteo.site.val$local%>%
window(start=start.val.date,
end=end.val.date) %>%
coredata)
data.frame(NSE=NSE.val,KGE=KGE.val) %>%
pandoc.table(style="simple",caption="Runoff Model Performance - Daily")| NSE | KGE |
|---|---|
| 0.3359 | 0.5735 |
4 Evaluation on Project
The concept will be to evaluated the GR5J model from the Rio Pulido en Vertedero with the parameters defined from the calibrated and validated model.
The model will be reevaluated with the new hypsometric curve & the local meteorology based on the report SRK 2024.
4.1 Hypsometric Curve
proj_shp<-st_read(here::here("data","csv","Cuenca Proyecto.kml")) %>%
st_transform(.,crs=4326)
dem_ras<-crop(wat_gtopo30_ras,proj_shp) %>% mask(.,proj_shp)
dem_tbl<-tidyterra::fortify(dem_ras) %>%
dplyr::rename("z"=3) %>%
dplyr::filter(complete.cases(z))
plot(dem_ras)
x_ecdf<- dem_tbl$z %>% ecdf()
z_val<-seq(min(dem_tbl$z,na.rm = T),
max(dem_tbl$z,na.rm = T),
length.out = 200)
z_hypso<-approx(x = x_ecdf(z_val),
y=z_val,
xout=seq(0.001,1,length.out = 101))
z_hypso$y[1]<-min(dem_tbl$z,na.rm = T)
z_hypso$y<-na.approx(z_hypso$y)
plot(z_hypso,main="Hypsometric curve",xlab="percentage of area [%]",
ylab="Elevation[masl]")
4.2 Meteorology
meteo_info_z<-meteo_info_fn(x = mean(dem_tbl$x),
y = mean(dem_tbl$y),
z = mean(dem_tbl$z))
plot(meteo_info_z,main="available information")
4.3 Runoff
InputsModel.eval<-CreateInputsModel(
FUN_MOD = RunModel_CemaNeigeGR5J,
DatesR = time(Meteo.site.eval) %>% as.POSIXlt(),
Precip = Meteo.site.eval$prcp%>% coredata,
PotEvap = Meteo.site.eval$pe %>% coredata,
TempMax = Meteo.site.eval$tmax%>% coredata,
TempMin = Meteo.site.eval$tmin%>% coredata,
TempMean = Meteo.site.eval$tavg%>%coredata,
ZInputs = mean(dem_tbl$z),
PrecipScale = T,
HypsoData = na.omit(z_hypso$y),
NLayers=4)
#capturing at different elvations
meteo_info_lst<-lapply(InputsModel.eval$ZLayers,FUN=function(x){
meteo_info_fn(x = mean(dem_tbl$x),
y = mean(dem_tbl$y),
z = x) %>%
window(start=start.eval.date,end=end.eval.date)
})
#prcp
prcp_lst<-lapply(meteo_info_lst,FUN=function(x){x$prcp %>% coredata()})
names(prcp_lst)<-paste0("L",1:length(prcp_lst))
InputsModel.eval$LayerPrecip<-prcp_lst
#air temperature
tavg_lst<-lapply(meteo_info_lst,FUN=function(x){x$tavg %>% coredata()})
names(tavg_lst)<-paste0("L",1:length(tavg_lst))
InputsModel.eval$LayerTempMean<-tavg_lst
#Snow Fraction as Kienzle
sf_lst<-lapply(meteo_info_lst,FUN=function(x){
# if_else(x$prcp<=0,0,coredata(x$snowfall/x$prcp))
if_else(x$prcp<=0,0, 1-p_rain_fn(Tav = x$tavg) %>%
coredata())
})
names(sf_lst)<-paste0("L",1:length(sf_lst))
InputsModel.eval$LayerFracSolidPrecip<-sf_lst
Ind_Run.eval <- seq(grep(start.eval.date, time(Meteo.site.eval)),
grep(end.eval.date, time(Meteo.site.eval)))
RunOptions.eval <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR5J,
InputsModel = InputsModel.eval,
IndPeriod_Run = Ind_Run.eval,
IniStates = NULL, IniResLevels = NULL,
IndPeriod_WarmUp =1:60,
IsHyst = F )
OutputsModel.eval <- RunModel_CemaNeigeGR5J(InputsModel = InputsModel.eval,
RunOptions = RunOptions.eval,
Param = Param_eval)4.4 Results
The results are estimated in mm/d, these values are transformed to m3/s when these values are transformed to m3/s when values are multiplied by their respective watersheds.
site_q_mm_d<-zoo(OutputsModel.eval$Qsim,
as.Date(OutputsModel.eval$DatesR))
daily2annual.avg(site_q_mm_d,FUN=sum) %>%
pandoc.table(round=1,caption="Monthly Runoff [mm/mon]")| Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 4.9 | 4.7 | 4.1 | 3.1 | 2.6 | 2.3 | 2.1 | 1.8 | 1.6 | 1.6 | 2.3 | 3.9 | 35.1 |
daily2monthly.V2(site_q_mm_d,FUN=sum) %>%
zoo2matrix() %>%
boxplot(ylab="mm/mon",xlab=NULL,main="Monthly Boxplot for Runoff at the site mm/mon")
#maximum per year
site_q_mm_d_max_yr<-coredata(IDF(site_q_mm_d,1))[,1]
site_q_m3_s_max_yr<-ud.convert(site_q_mm_d_max_yr*1,"mm/d*km2","m3/s")
site_watersheds<-tibble::tribble(
~Cuenca, ~area_km2,
"Canal Derivación Superior", 1.89,
"Canal Derivación Inferior", 7.12,
"Qda. N°1(Situación Sin P, Sector Botadero Norte)", 3.83,
"Qda. N°2(Situación Sin P, Sector Botadero Norte)", 0.63,
"Qda. N°3(Situación Sin P, Sector Botadero Norte)", 0.48,
"Qda. N°4(Situación Sin P, Sector Botadero Norte)", 0.44,
"Qda. N°5(Situación Sin P, Sector Botadero Norte)", 0.74,
"Afluente a Qda. N°5(Situación Sin P, Sector Botadero Norte)", 0.29,
"Botadero Norte (Situación Con P)", 1.94
)
max_flow_fa<-FA(site_q_m3_s_max_yr)$results %>%
mutate(RP=RP.text(Prob,dry.wet=F)) %>%
dplyr::filter(Prob>=0.5) %>%
arrange(Prob) %>%
dplyr::select(RP,LN) %>%
dplyr::rename("max_daily_m3.s"=2)cross_join(site_watersheds,
max_flow_fa) %>%
mutate(Cuenca=factor(Cuenca,site_watersheds$Cuenca,ordered = T)) %>%
mutate(value=area_km2*max_daily_m3.s) %>%
dplyr::filter(RP%in%c(100,200,500)) %>%
reshape2::dcast(data=.,formula=Cuenca+area_km2~RP,value.var = "value") %>%
pandoc.table(round=3,
caption="maximum daily average values in m3/s for return periods of 100, 200 and 500 yrs.")| Cuenca | area_km2 | 100 | 200 | 500 |
|---|---|---|---|---|
| Canal Derivación Superior | 1.89 | 0.025 | 0.031 | 0.041 |
| Canal Derivación Inferior | 7.12 | 0.094 | 0.117 | 0.153 |
| Qda. N°1(Situación Sin P, Sector Botadero Norte) | 3.83 | 0.051 | 0.063 | 0.082 |
| Qda. N°2(Situación Sin P, Sector Botadero Norte) | 0.63 | 0.008 | 0.01 | 0.014 |
| Qda. N°3(Situación Sin P, Sector Botadero Norte) | 0.48 | 0.006 | 0.008 | 0.01 |
| Qda. N°4(Situación Sin P, Sector Botadero Norte) | 0.44 | 0.006 | 0.007 | 0.009 |
| Qda. N°5(Situación Sin P, Sector Botadero Norte) | 0.74 | 0.01 | 0.012 | 0.016 |
| Afluente a Qda. N°5(Situación Sin P, Sector Botadero Norte) | 0.29 | 0.004 | 0.005 | 0.006 |
| Botadero Norte (Situación Con P) | 1.94 | 0.026 | 0.032 | 0.042 |