Runoff Review
Victor Munoz S.
2024-11-06
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_")
output_z<-temp_adj_z(z=z,ts = ts_temp,z0 = 4323.77)
return(output_z)
}
#as Adjustment based on Project information AND climate change where
#elevation gradients are included
temp_adj_hist_cc_z2<-function(z,ts=NULL){
if(is.null(ts)){
#case that it is historical
ts_temp<-tk_zoo(site_temp_tbl[,1:4])
names(ts_temp)<-str_remove_all(names(ts_temp),"2m_")
output_z<-temp_adj_z(z=z,ts = ts_temp,z0 = 4323.77)
}else{
#case when TS is climate change info, that will adjusted
#based on elevation gradients
output_z<-temp_adj_z(z=z,ts = ts,z0 = 4323.77)
}
return(output_z)
}
# Test to check that function works
# temp_hist_z<-ts_temp
#
# temp_cc_z<-cc_info_lst[[5]] %>%
# dplyr::select(tmin=tasmin_qdm,tavg=tas_qdm,tmax=tasmax_qdm,date) %>%
# tk_zoo(silent=T)
#
# temp_adj_hist_cc_z2(z = 4500,ts = temp_cc_z) %>%
# daily2annual.avg(FUN=mean)
#
# temp_adj_hist_cc_z2(z = 4300,ts = temp_hist_z) %>%
# daily2annual.avg(FUN=mean)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
#includes elevation gradients on air temperature
nc2z_elev_fn<-function(Longitude=site.info$Longitude,
Latitude=site.info$Latitude,
source="ERA5-LAND",elev=3500){
prcp_z<-nc2z_fn(Longitude = Longitude,
Latitude = 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
2.6.1 Historical
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)
}#Precipitation without corrections
meteo_info_z<-nc2z_elev_fn(Longitude =x1,
Latitude = y1,
elev=z1,
source = "ERA5-LAND")
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 |
2.6.2 Under Climate Change
meteo_info_cc_fn<-function(z1,ts_temp=temp_cc_z,ts_prcp=prcp_cc_z){
#Set air temperature
meteo_info_z<-temp_adj_hist_cc_z2(z = z1,ts = ts_temp)
meteo_info_z$prcp<-prcp_cc_z
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<-fortify(meteo_info_z) %>%
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))) %>%
#elevation correction based on Projecto Bias Correction
mutate(rainfall=rainfall*rf_cc_fn(z1),
snowfall=snowfall*sf_cc_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)
}
#Test!
# n_model<-c(3)
#
# prcp_cc_z<-cc_info_lst[[n_model]] %>%
# dplyr::select(prcp=pr_qdm,date) %>%
# tk_zoo(silent=T)
#
# temp_cc_z<-cc_info_lst[[n_model]] %>%
# dplyr::select(tmin=tasmin_qdm,tavg=tas_qdm,tmax=tasmax_qdm,date) %>%
# tk_zoo(silent=T)
#
#
# meteo_info_cc_z<-meteo_info_cc_fn(z1 = proyecto_st_elev,
# ts_temp = temp_cc_z,ts_prcp = prcp_cc_z) %>%
# fortify.zoo() %>%
# dplyr::select(tavg,prcp,Index) %>%
# tk_zoo(silent=T)
#
# rbind(
# meteo_info_cc_z$tavg%>%
# daily2annual.avg(FUN=mean),
#
# meteo_info_cc_z$prcp%>%
# daily2annual.avg(FUN=sum)
# )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 - 1992
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("1992-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)
rbind(
daily2annual.avg(Meteo.site.cal[,c("tmax","tavg","tmin")],
FUN=mean),
daily2annual.avg(Meteo.site.cal[,c("rainfall","snowfall","prcp")],
FUN=sum)
) %>%
pandoc.table(round=1,split.table=Inf,
caption=
"Rio Pulido en Vertedero - mean monthly values")| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| tmax | 21.4 | 20.7 | 19.7 | 16.7 | 13.5 | 11.7 | 11.2 | 12.7 | 13.9 | 17 | 19.2 | 20.8 | 16.5 |
| tavg | 8.8 | 8.3 | 7.3 | 4.8 | 2.1 | 0 | -0.7 | 0.9 | 2 | 4.7 | 6.6 | 8.2 | 4.4 |
| tmin | -6.2 | -6.2 | -7 | -8.7 | -10.8 | -13.1 | -14.2 | -12.5 | -11.7 | -9.5 | -8.3 | -6.9 | -9.6 |
| rainfall | 2.4 | 3.1 | 2.3 | 1.1 | 2 | 0.4 | 0.7 | 1.7 | 1.7 | 1.6 | 1.3 | 2.7 | 21 |
| snowfall | 0 | 0 | 0 | 4.3 | 21.3 | 26.1 | 60 | 27.9 | 13.8 | 1.9 | 0.2 | 0 | 155.4 |
| prcp | 2.4 | 3.1 | 2.3 | 5.4 | 23.3 | 26.4 | 60.7 | 29.5 | 15.5 | 3.5 | 1.5 | 2.7 | 176.3 |
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_CemaNeigeGR4J,
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_CemaNeigeGR4J,
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
# start_par_list_gr4j<-matrix(
# nrow=1,ncol=6,
# c(0.0001,-0.110,238,20,0.996,0.724))
start_param_gr4j<-matrix(
nrow=1,ncol=6,
c(0.0001,-0.110,238,20,0.996,0.724))
search_ranges_gr4j<-matrix(
nrow=2,ncol=6,
c(0,-10903,0,0.5,0,0,
21807,10903,21807,30,1,109.04),byrow = T)
# start_param_gr5j<-matrix(
# nrow=1,ncol=7,
# c(0.0001,-0.110,238,20,0.345,0.996,0.724))
#
#
# search_ranges_gr5j<-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_gr4j,
FUN_MOD = RunModel_CemaNeigeGR4J,
FUN_CALIB = Calibration_Michel,
StartParamList = start_param_gr4j,
# SearchRanges =search_ranges,
IsHyst = F)
OutputsCalib.cal <- Calibration_Michel(InputsModel = InputsModel.cal,
RunOptions = RunOptions.cal,
InputsCrit = InputsCrit.cal,
CalibOptions = CalibOptions.cal,
FUN_MOD = RunModel_CemaNeigeGR4J)
Param_eval <- OutputsCalib.cal$ParamFinalR
OutputsModel.cal <- RunModel_CemaNeigeGR4J(
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.6764 | 0.8399 |
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 | 2.35 | 3.1 | 2.33 | 5.38 | 23.3 | 26.44 | 60.68 | 29.54 | 15.52 | 3.47 | 1.49 | 2.73 | 176.3 |
| P_eff | 24.45 | 14.25 | 8.45 | 2.83 | 4.5 | 2.65 | 4.51 | 9.8 | 13.21 | 34.11 | 38.01 | 36.67 | 193.4 |
| E | 76.23 | 62.41 | 54.98 | 34.49 | 20.67 | 12.5 | 11.8 | 19.74 | 27.66 | 47.64 | 61 | 74.43 | 503.5 |
| P net(Pn) | 0.01 | 0 | 0.15 | 0.13 | 0.52 | 0.23 | 1.13 | 1.23 | 0.68 | 3.93 | 2.79 | 1.24 | 12.03 |
| P storage(Ps) | 0 | 0 | 0 | 0 | 0.01 | 0.01 | 0.02 | 0.01 | 0.02 | 0.02 | 0.01 | 0.01 | 0.1 |
| Pn-Ps | 0.01 | 0 | 0.14 | 0.13 | 0.51 | 0.22 | 1.11 | 1.22 | 0.66 | 3.91 | 2.78 | 1.23 | 11.93 |
| aft.Perc.rout.(Pr) | 0.01 | 0 | 0.14 | 0.13 | 0.51 | 0.22 | 1.11 | 1.22 | 0.66 | 3.91 | 2.79 | 1.23 | 11.93 |
| Exchange (F) | 2.18 | 1.69 | 1.64 | 1.47 | 1.42 | 1.58 | 1.53 | 1.54 | 1.51 | 1.63 | 1.99 | 2.38 | 20.55 |
| Q | 3.84 | 2.71 | 2.5 | 2.16 | 2.07 | 2.42 | 2.3 | 2.34 | 2.31 | 2.61 | 3.63 | 4.62 | 33.51 |
| Q_obs | 4.24 | 3.15 | 2.92 | 2.27 | 2.54 | 2.11 | 2.3 | 2.26 | 2.01 | 2.48 | 3.07 | 3.46 | 32.79 |
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] | 2.62 |
| GR4J X3 | routing store capacity [mm] | 141.5 |
| GR4J X4 | unit hydrograph time constant [d] | 20.34 |
| CemaNeige X1 | weighting coefficient for snow pack thermal state [-] | 1 |
| CemaNeige X2 | degree-day melt coefficient [mm/°C/d] | 0.72 |
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 1995 - 2010
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")
start.val.date<-as.Date("1992-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_CemaNeigeGR4J,
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_CemaNeigeGR4J,
InputsModel = InputsModel.val,
IndPeriod_Run = Ind_Run.val,
IniStates = NULL, IniResLevels = NULL,
IndPeriod_WarmUp =1:60,
IsHyst = F )
OutputsModel.val <- RunModel_CemaNeigeGR4J(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.5055 | 0.7284 |
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_CemaNeigeGR4J,
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_CemaNeigeGR4J,
InputsModel = InputsModel.eval,
IndPeriod_Run = Ind_Run.eval,
IniStates = NULL, IniResLevels = NULL,
IndPeriod_WarmUp =1:60,
IsHyst = F )
OutputsModel.eval <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel.eval,
RunOptions = RunOptions.eval,
Param = Param_eval)4.4 Results
met_gl<-data.frame(
Date=OutputsModel.eval$DatesR %>% as.Date(),
P=OutputsModel.eval$Precip,
P_eff=(0.25*OutputsModel.eval$CemaNeigeLayers$Layer01$PliqAndMelt+
0.25*OutputsModel.eval$CemaNeigeLayers$Layer02$PliqAndMelt+
0.25*OutputsModel.eval$CemaNeigeLayers$Layer03$PliqAndMelt+
0.25*OutputsModel.eval$CemaNeigeLayers$Layer04$PliqAndMelt),
E=OutputsModel.eval$PotEvap
) %>%
mutate(
# Pn=if_else(P>=E,P-E,0),
# En=if_else(P>=E,0,E-P),
`P net(Pn)`=OutputsModel.eval$Pn,
`P storage(Ps)`=OutputsModel.eval$Ps,
`Pn-Ps`=`P net(Pn)`-`P storage(Ps)`,
# Prod=OutputsModel1$Prod,
# AE=OutputsModel1$AE,
`aft.Perc.rout.(Pr)`=OutputsModel.eval$PR,
`Exchange (F)`=OutputsModel.eval$AExch,
Q=OutputsModel.eval$Qsim
)
wb_gl.z<-tk_zoo(met_gl,silent=T)
wb_gl.syn<-merge(wb_gl.z
) %>% daily2annual.avg(FUN=sum) %>%
data.frame()
pandoc.table(wb_gl.syn,
round=2,split.table=Inf,missing="-",
caption="Watershed Summary from 1980 to 2020")
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")
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)
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 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2.3 | 2.3 | 2.4 | 2.1 | 1.9 | 1.9 | 1.8 | 1.7 | 1.6 | 1.7 | 1.8 | 2.1 | 23.7 |
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,2) %>%
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.02 | 0.028 | 0.043 |
| Canal Derivación Inferior | 7.12 | 0.074 | 0.104 | 0.163 |
| Qda. N°1(Situación Sin P, Sector Botadero Norte) | 3.83 | 0.04 | 0.056 | 0.088 |
| Qda. N°2(Situación Sin P, Sector Botadero Norte) | 0.63 | 0.007 | 0.009 | 0.014 |
| Qda. N°3(Situación Sin P, Sector Botadero Norte) | 0.48 | 0.005 | 0.007 | 0.011 |
| Qda. N°4(Situación Sin P, Sector Botadero Norte) | 0.44 | 0.005 | 0.006 | 0.01 |
| Qda. N°5(Situación Sin P, Sector Botadero Norte) | 0.74 | 0.008 | 0.011 | 0.017 |
| Afluente a Qda. N°5(Situación Sin P, Sector Botadero Norte) | 0.29 | 0.003 | 0.004 | 0.007 |
| Botadero Norte (Situación Con P) | 1.94 | 0.02 | 0.028 | 0.044 |
5 Climate Change
climate_change_files<-dir(
here::here("data","csv","GCM"),
recursive = T,full.names = T)
# climate_change_file<-climate_change_files[[1]]
#division of the models
cc_info_tbl<-lapply(climate_change_files,
FUN=function(climate_change_file){
gcm_records_tbl <- read_csv(
climate_change_file,
col_types = cols(date = col_date(
format = "%Y-%m-%d"))) %>%
data.table()
return(gcm_records_tbl)
}) %>% data.table::rbindlist()
cc_info_lst<-split(cc_info_tbl,cc_info_tbl$id)5.1 Runoff models
meteo_info_z<-nc2z_elev_fn(Longitude =x1,
Latitude = y1,
elev=z1,
source = "ERA5-LAND")
meteo_info_z$prcp %>%
daily2annual.avg(.,FUN=sum) %>%
pandoc.table(style="simple",split.table=Inf,round=1,
caption="Raw precipitation values from ERA5-Land")| Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 28.7 | 25 | 15.8 | 7 | 16.2 | 18.1 | 16.5 | 12.7 | 8.1 | 3.4 | 3.3 | 11.1 | 165.8 |
#Precipitation WITH corrections
meteo_info_z<-meteo_info_fn(x1,y1,z1)
meteo_info_z$prcp %>%
daily2annual.avg(.,FUN=sum)%>%
pandoc.table(style="simple",split.table=Inf,round=1,
caption="Corrected precipitation values based on SRK 2024")| Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 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 |
# #Climate change model
n_model<-c(2)
#
temp_cc_z<-cc_info_lst[[n_model]] %>%
dplyr::select(tmin=tasmin_qdm,
tavg=tas_qdm,
tmax=tasmax_qdm,date) %>%
tk_zoo(silent=T) %>%
na.approx(maxgap = 15)
prcp_cc_z<-cc_info_lst[[n_model]] %>%
dplyr::select(prcp=pr_qdm,date) %>%
tk_zoo(silent=T)
#
rbind(
daily2annual.avg(prcp_cc_z,FUN=sum),
daily2annual.avg(temp_cc_z,FUN=mean)
) %>%
pandoc.table(style="simple",split.table=Inf,round=1,
caption=paste0("Climate change prcp bias corrected by SRK 2024 / Model: ",cc_info_lst[[n_model]]$id[1]))| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| prcp | 21.8 | 21 | 14.7 | 0.7 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 3.1 | 14.5 | 76.2 |
| tmin | -1.5 | -1.4 | -2.3 | -4.5 | -6.5 | -9.8 | -10.5 | -9 | -6.6 | -5.6 | -3.9 | -2.1 | -5.3 |
| tavg | 6.2 | 6 | 5 | 2.9 | 0.4 | -2.9 | -3.8 | -1.4 | 0.6 | 1.7 | 3.7 | 6 | 2 |
| tmax | 13.4 | 12.5 | 11.9 | 10.5 | 8.1 | 5.7 | 5.2 | 7.5 | 8.8 | 9.8 | 11.6 | 14 | 9.9 |
# cc_out_tbl<-pbapply::pblapply(
# seq_along(cc_info_lst),
# FUN=function(n_model){
#
# # n_model<-53
# print(n_model)
#
# #Meteorology
# prcp_cc_z<-cc_info_lst[[n_model]] %>%
# dplyr::select(prcp=pr_qdm,date) %>%
# tk_zoo(silent=T)%>%
# na.approx(maxgap = 15)
#
# prcp_cc_z<<-prcp_cc_z
#
# temp_cc_z<-cc_info_lst[[n_model]] %>%
# dplyr::select(tmin=tasmin_qdm,
# tavg=tas_qdm,
# tmax=tasmax_qdm,date) %>%
# tk_zoo(silent=T) %>%
# na.approx(maxgap = 15)
#
# temp_cc_z<<-temp_cc_z
#
# #Start - End
# start_cc_date<-head(zoo::index(prcp_cc_z),1)
#
# end_cc_date<-tail(zoo::index(prcp_cc_z),1)
#
# #Implementation fo meteorology at elevation of station Proyecto
# meteo_info_cc_z<-meteo_info_cc_fn(z1 = z1,
# ts_temp = temp_cc_z,
# ts_prcp = prcp_cc_z) %>%
# na.approx(maxgap = 15)
#
# dwi_tmax_tmin<-sum(dwi(meteo_info_cc_z$tmax),dwi(meteo_info_cc_z$tmin))
#
# if(dwi_tmax_tmin>0){
#
# InputsModel.cc<-CreateInputsModel(
# FUN_MOD = RunModel_CemaNeigeGR4J,
# DatesR = time(meteo_info_cc_z) %>% as.POSIXlt(),
# Precip = meteo_info_cc_z$prcp%>% coredata,
# PotEvap = meteo_info_cc_z$pe %>% coredata,
# TempMax = meteo_info_cc_z$tmax%>% coredata,
# TempMin = meteo_info_cc_z$tmin%>% coredata,
# TempMean = meteo_info_cc_z$tavg%>%coredata,
# ZInputs = proyecto_st_elev,
# PrecipScale = T,
# HypsoData = na.omit(z_hypso$y),
# NLayers=4)
#
# #capturing at different elvations
# meteo_info_cc_lst<-lapply(InputsModel.cc$ZLayers,FUN=function(x){
#
# meteo_info_cc_fn(z1 = x,
# ts_temp = temp_cc_z,
# ts_prcp = prcp_cc_z) %>%
# window(start=start_cc_date,end=end_cc_date)
#
# })
#
# #prcp
# prcp_lst<-lapply(meteo_info_cc_lst,
# FUN=function(x){x$prcp %>% coredata()})
#
# names(prcp_lst)<-paste0("L",1:length(prcp_lst))
#
# InputsModel.cc$LayerPrecip<-prcp_lst
#
# #air temperature
# tavg_lst<-lapply(meteo_info_cc_lst,
# FUN=function(x){x$tavg %>% coredata()})
#
# names(tavg_lst)<-paste0("L",1:length(tavg_lst))
#
# InputsModel.cc$LayerTempMean<-tavg_lst
#
# #Snow Fraction as Kienzle
# sf_lst<-lapply(meteo_info_cc_lst,
# FUN=function(x){
#
# #Zero degree partition
# # if_else(x$prcp<=0,0,coredata(x$snowfall/x$prcp))
#
# #Kienzle
# if_else(x$prcp<=0,0, 1-p_rain_fn(Tav = x$tavg) %>%
# coredata())
#
# })
#
# names(sf_lst)<-paste0("L",1:length(sf_lst))
#
# InputsModel.cc$LayerFracSolidPrecip<-sf_lst
#
#
# Ind_Run.cc <- seq(grep(start_cc_date, time(meteo_info_cc_z)),
# grep(end_cc_date, time(meteo_info_cc_z)))
#
# RunOptions.cc <- CreateRunOptions(
# FUN_MOD = RunModel_CemaNeigeGR4J,
# InputsModel = InputsModel.cc,
# IndPeriod_Run = Ind_Run.cc,
# IniStates = NULL, IniResLevels = NULL,
# IndPeriod_WarmUp =1:60,
# IsHyst = F )
#
# OutputsModel.cc <- RunModel_CemaNeigeGR4J(
# InputsModel = InputsModel.cc,
# RunOptions = RunOptions.cc,
# Param = Param_eval)
#
# out_tbl<-data.frame(date=OutputsModel.cc$DatesR,
# q_mm.d=OutputsModel.cc$Qsim,
# start_date=start_cc_date,
# end_date=end_cc_date,
# period=cc_info_lst[[n_model]]$period[1],
# scenario=cc_info_lst[[n_model]]$scenario[1],
# model=cc_info_lst[[n_model]]$model[1]) %>%
# as.data.table()
#
# return(out_tbl)
#
# }else{
#
# return()
# }
#
# }) %>%
# rbindlist()
# arrow::write_parquet(x = cc_out_tbl,
# sink = here::here("outputs","cc_flow_results.par"))5.2 Resutts
library("nsRFA")
gev_lmom<-function(a,rp=100){
ll <- Lmoments(a)
par <- par.GEV(ll[1], ll[2], ll[4])
result <- invF.GEV(1-(1/rp),
xi = par$xi,alfa = par$alfa, k = par$k) %>%
as.vector()
return(result)
}
cc_out_tbl<-arrow::read_parquet(here::here("outputs",
"cc_flow_results.par"))
cc_out_res<-cc_out_tbl %>%
mutate(yr=lubridate::year(date)) %>%
group_by(model,scenario,period,start_date,end_date,yr) %>%
summarize(q_mm.d=max(q_mm.d)) %>%
dplyr::filter(yr>=lubridate::year(start_date)+2) %>%
group_by(model,scenario,period) %>%
summarize(q_mm.d_100y=gev_lmom(q_mm.d,rp=100),
q_mm.d_150y=gev_lmom(q_mm.d,rp=150),
q_mm.d_500y=gev_lmom(q_mm.d,rp=500)) %>%
mutate(q_m3_s_100y=ud.convert(q_mm.d_100y*1,"mm/d*km2","m3/s"),
q_m3_s_150y=ud.convert(q_mm.d_150y*1,"mm/d*km2","m3/s"),
q_m3_s_500y=ud.convert(q_mm.d_500y*1,"mm/d*km2","m3/s")) %>%
dplyr::select(-q_mm.d_100y,-q_mm.d_150y,-q_mm.d_500y) %>%
reshape2::melt(id.vars=c("model","scenario","period")) %>%
mutate(variable=str_remove_all(variable,"q_m3_s_") %>%
str_remove_all(.,"y") %>% as.numeric())
cc_out_hist<-cc_out_res%>%
group_by(scenario,period,variable) %>%
summarize(value=median(value)) %>%
dplyr::filter(period=="historical") %>%
rename(historical="value") %>%
ungroup() %>%
dplyr::select(-scenario,-period)
cc_tbl_rp<-cc_out_res%>%
group_by(scenario,period,variable) %>%
summarize(value=median(value)) %>%
dplyr::filter(period!="historical") %>%
rename(projected="value") %>%
left_join(.,cc_out_hist)%>%
mutate(rc=(projected-historical)/historical) %>%
rename(rp=variable) %>%
mutate(rp=paste0("1:",rp," yrs."))
ggplot(data=cc_tbl_rp,aes(x=period,y=rc))+
geom_bar(stat = "Identity",position = "dodge")+
facet_grid(scenario~rp)+
geom_text(aes(label=paste0(round(100*rc,0),"%")),nudge_y = 0.1,
size=2.5)+
theme_bw()+
labs(x="period",
y="Rate of Change with respect\n to historical Conditions [%]",
title="Median of GCM models for different return periods",
caption="Note: probabilistic distribution GEV with L-moments")+
scale_y_continuous(labels=scales::percent)