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 #Proyecto St
Latitude = -26.076496 #Proyecto St

#Get elevation from google maps
Elevation= 4323.77 #Proyecto St

site.info<-data.frame(Location="Site - Proyecto St",
                      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 - Proyecto St	-68.98	-26.08	4324

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]")

dwi.graph(flow_dly_z,start.year = 1950,
          station.name = flow_idx$nombre_estacion %>% str_wrap(.,12))

daily2monthly.V2(flow_dly_z,FUN=sum) %>%
  window(start=as.Date("1980-01-01")) %>% 
  plot(type="b")

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)

Total runoff mm/mon
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){
  
  #y<-1########################################################ojooooooooo
  
  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

Table continues below
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

Table continues below
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

Table continues below
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

Table continues below
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

Table continues below
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

Table continues below
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

Table continues below
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 St=4544.743 masl(although it would have been less confusing using the expression for the Proyecto St)
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 Proyecto St 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){ #this elevation value is overriden if another elev value is provided when calling the function
        
        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)
        
}


#trial function to test if x is effectively overriden. it is. z can be 

trial_fun<-function(x,z=2){
 
  y<-x*2+z
  
  return(y)
  
}

2.5.1 Rainfall

Values based on range at 3500 to 4800

proyecto_st_elev<-site.info$Elevation

# Rainfall factor function
#rf_fn attempts to bias-correct raw ERA5-Land values based solely on elevation (i.e., ERA5-Land data at a given coord is downloaded and automatically bias corected for the elev of that given coordinate - see Fig 5.18 in Baseline report) 
rf_fn<-function(z){
  
  pmin(pmax(-0.8457+z*0.0002612,0.148),0.367)
  
}

# Rainfall factor function correction for climate change based on Proyecto St
# CC rainfall is a bias corrected synthetic timeseries developed at the Proyecto St coordinate 
# Correction factor rf_fn does not really apply, as it was developed to bias correct ERA5-Land data exclusively
# Assuming that cc precip had been developed from ERA5Land data, this would mean that it would have had to be corrected by rf_fn(proyecto_st_elev)
#this is because cc data was done at the Proyecto St coordinate
#that's where this correction rf_cc_fn comes from

rf_cc_fn<-function(z){
  
  rf_fn(z)/rf_fn(proyecto_st_elev)
  
}

2.5.2 Snowfall

Values based on range at 3500 to 4800

# Snowfall factor function
sf_fn<-function(z){
  
  pmin(pmax(1.54-z*0.000182,0.65),0.85)
}

# Snowfall factor function correction for climate change based on Proyecto St

sf_cc_fn<-function(z){
  
  sf_fn(z)/sf_fn(proyecto_st_elev)
  
}

2.6 Meteo Compiled

2.6.1 Historical

applied to site conditions

x1<-site.info$Longitude
y1<-site.info$Latitude
z1<-site.info$Elevation

#function meteo_info_fn returns a table with baseline dly values of prcp, tavg, tmax, tmin, rainfall, snowfall and pe
meteo_info_fn<-function(x1,y1,z1){
  
  #Use next three lines only for testing the function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  #x1<-site.info$Longitude
  #y1<-site.info$Latitude
  #z1<-site.info$Elevation
  
  data('constants')
  
  #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$tmax[(meteo_info_z$tmax<meteo_info_z$tmin)]<-NA
  
  meteo_info_z<-na.approx(meteo_info_z)
  
  #pe_z must be a zoo time series with date and dly evap
  
  #1. Oudin method for evaporation
  # 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))
  
  #2. Hargreaves-Samani method for evaporation
  
  ##constants argument of HS function
  data('constants') #this uploads a list of constants used in HS method
  constants_hs<-constants
  constants_hs$Elev<-z1
  constants_hs$lat_rad =y1*pi/180

  ##data argument of HS function
  data_hs <- ReadInputs(varnames = c("Tmax","Tmin"),
           meteo_info_z %>% 
             fortify.zoo() %>% 
                   dplyr::rename("Tmax"="tmax",
                                 "Tmin"="tmin",
                                 "date"="Index") %>% 
                   mutate(Year=lubridate::year(date),
                          Month=lubridate::month(date),
                          Day=lubridate::day(date)), 
           stopmissing=c(10,10,3),
           constants = constants_hs,
           timestep = "daily",
           interp_missing_days = FALSE, 
           interp_missing_entries = FALSE, 
           interp_abnormal = FALSE, 
           missing_method = NULL, 
           abnormal_method = NULL)

  
  pe_z <- Evapotranspiration::ET.HargreavesSamani(
    data_hs,
    constants_hs,
    ts = "daily",             # Time scale, e.g., "daily" or "monthly"
    message = "yes",          # Whether to display messages during processing
    AdditionalStats="yes",
    save.csv = "no"           # Whether to save results as a CSV
  ) %>%
    .$ET.Daily
  
  pe_z<-ifelse.zoo(pe_z<0,0,pe_z) #removal of any evap values below 0

  #Updated tk_zoo with pe_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)  
  

}

#function meteo_info_fn applied to site coordinates

#Precipitation without corrections
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")

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")

Corrected precipitation values based on SRK 2024
Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Ann
8.1	7.3	5.2	4.6	12	13.6	12.4	9.5	5.8	1.9	1.2	3.4	85

# #Climate change model
# n_model<-c(2)
# 
# prcp_cc_z<-cc_info_lst[[n_model]] %>% 
#    dplyr::select(prcp=pr_qdm,date) %>% 
#    tk_zoo(silent=T)
# 
# daily2annual.avg(prcp_cc_z,FUN=sum) %>% 
#     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]))

meteo_info_z<-meteo_info_fn(x1,y1,z1)

dwi.graph(meteo_info_z)

plot(meteo_info_z, 
     main="meteorology at site")

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")

site meteorology - mean monthly values 1980 to 2021
Station	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Ann
tmax	13.5	12.7	11.6	8.6	5.5	3.7	3.4	4.9	6.4	8.9	11	12.7	8.6
tavg	4.8	4.2	3.2	0.7	-2	-3.9	-4.5	-3	-1.6	0.6	2.4	4.1	0.4
tmin	-3.7	-3.9	-4.8	-6.5	-8.6	-10.7	-11.5	-10.1	-9	-7.4	-6.4	-4.8	-7.3
rainfall	8.1	7	4.1	0.4	0.1	0	0	0	0.2	0.4	0.8	3	24.1
snowfall	0	0.3	1.1	4.3	11.9	13.6	12.4	9.5	5.6	1.5	0.5	0.3	61
prcp	8.1	7.3	5.2	4.6	12	13.6	12.4	9.5	5.8	1.9	1.2	3.4	85

2.6.2 Under Climate Change

#function meteo_info_fn returns a table with cc dly values of prcp, tavg, tmax, tmin, rainfall, snowfall and pe

#unless specified otherwise when calling the function, the values ts_temp=temp_cc_z and ts_prcp=prcp_cc_z will be used

meteo_info_cc_fn<-function(z1,ts_temp=temp_cc_z,ts_prcp=prcp_cc_z){ 
  
  #Two lines only intended to test function
  #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
  
  #1. Oudin method for evaporation  
  # 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))


  #2. Hargreaves-Samani method for evaporation
  
  ##constants argument of HS function
  data('constants') #this uploads a list of constants used in HS method
  constants_hs<-constants
  constants_hs$Elev<-z1
  constants_hs$lat_rad =y1*pi/180

  ##data argument of HS function
  data_hs <- ReadInputs(varnames = c("Tmax","Tmin"),
           meteo_info_z %>% 
             fortify.zoo() %>% 
                   dplyr::rename("Tmax"="tmax",
                                 "Tmin"="tmin",
                                 "date"="Index") %>% 
                   mutate(Year=lubridate::year(date),
                          Month=lubridate::month(date),
                          Day=lubridate::day(date)), 
           stopmissing=c(10,10,3),
           constants = constants_hs,
           timestep = "daily",
           interp_missing_days = FALSE, 
           interp_missing_entries = FALSE, 
           interp_abnormal = FALSE, 
           missing_method = NULL, 
           abnormal_method = NULL)

  
  pe_z <- Evapotranspiration::ET.HargreavesSamani(
    data_hs,
    constants_hs,
    ts = "daily",             # Time scale, e.g., "daily" or "monthly"
    message = "yes",          # Whether to display messages during processing
    AdditionalStats="yes",
    save.csv = "no"           # Whether to save results as a CSV
  ) %>%
    .$ET.Daily
  
  pe_z<-ifelse.zoo(pe_z<0,0,pe_z) #removal of any evap values below 0
    
  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 Proyecto 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 GR4J - Rio Pulido en Vertedero

id_sel<-wat_ras_idx[2,]

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")

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

  dygraph(flow_dly_z[,as.character(id_sel$id)] %>% 
            window(start=as.Date("1980-01-01"))) %>% 
  dygraphs::dyRangeSelector()

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)

nlayers_project<-3 #number of layers for calibration and validation (must be the same in both)

InputsModel.cal<-CreateInputsModel(
  FUN_MOD = RunModel_CemaNeigeGR4J,
  DatesR = time(Meteo.site.cal) %>% as.POSIXlt(),
  Precip = Meteo.site.cal$prcp%>% coredata, #meteo.site.cal contains data at the centroind of the watershed where calibration is being performed
  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, #mean elevation of the watershed where calibration is being performed 
  PrecipScale = T,
  HypsoData = na.omit(z_hypso$y),
  NLayers=nlayers_project)

#capturing at different elevations  - here meteo.site.cal is corrected by elevation for the mean elevation of each of the nlayers_project
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"), #original 
  #Weights = list(0.4, 0.6)      #original
  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))

#seed for model parameters
start_param_gr4j<-matrix(
      nrow=1,ncol=6,
      c (0.000100,   1.230370, 464.788556,  21.656918,   0.967836,  27.835963)) 

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: values of calibrated parameters
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)+
title(sub = list(paste0('Calibration, period ',year(flow_start),'-',year(flow_end),', ',
                        'NSE=',round(NSE.cal,2),', ',
                        'KGE=',round(KGE.cal,2)),
                 cex = 0.9))

integer(0)

data.frame(NSE=NSE.cal,KGE=KGE.cal) %>% 
  pandoc.table(style="simple",caption="Runoff Model Performance - Daily")

Runoff Model Performance - Daily
NSE	KGE
0.5821	0.7946

3.2.4 Results

#vector of P_eff depending on number of layers in GRXJ model
P_eff_list<-pbapply::pblapply(1:nlayers_project,FUN=function(x){

  y<-if(x<10){
    paste0('OutputsModel.cal$CemaNeigeLayers$Layer0',x,'$PliqAndMelt')
  }else{
    paste0('OutputsModel.cal$CemaNeigeLayers$Layer',x,'$PliqAndMelt')
  }
  
  eval(parse(text=y))
})

  met_gl<-data.frame(
    Date=OutputsModel.cal$DatesR %>% as.Date(),
    P=OutputsModel.cal$Precip,
    P_eff=Reduce(`+`,P_eff_list)/nlayers_project, #averages PliqAndMelt of all layers of GRXJ model
    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")

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	4.35	4.86	3.35	2.1	7.43	3.95	2.72	10.07	29.56	40.34	34.67	51.91	195.3
E	611.6	475.9	422.3	257.6	160	131.3	146.2	193.3	257.7	405.2	524.2	611.3	4197
P net(Pn)	0	0	0	0	0.43	1.47	0	2.65	10.62	11.8	11.27	7.07	45.31
P storage(Ps)	0	0	0	0	0	0	0	0	0.01	0.01	0	0.01	0.02
Pn-Ps	0	0	0	0	0.43	1.47	0	2.65	10.62	11.79	11.26	7.06	45.29
aft.Perc.rout.(Pr)	0	0	0	0	0.43	1.47	0	2.65	10.62	11.79	11.26	7.06	45.29
Exchange (F)	-1.72	-1.18	-1.17	-1.04	-1.02	-1.29	-1.17	-1.05	-1.31	-1.55	-2.13	-1.92	-16.54
Q	4.08	3.09	2.91	2.47	2.27	3.56	2.97	2.5	2.66	3.44	4.48	4.33	38.75
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")

GR4J with snowmelt parameters
Parameter	definition	value
GR4J X1	production store capacity [mm]	0
GR4J X2	intercatchment exchange coefficient [mm/d]	-8.85
GR4J X3	routing store capacity [mm]	859.7
GR4J X4	unit hydrograph time constant [d]	20
CemaNeige X1	weighting coefficient for snow pack thermal state [-]	0.98
CemaNeige X2	degree-day melt coefficient [mm/°C/d]	26.95

  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=nlayers_project)

#capturing at different elevations  (correcting Meteo.site.eval by elevation for the mean elev of nlayers_project)
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)+
title(sub = list(paste0('Validation, period ',year(flow_start),'-',year(flow_end),', ',
                        'NSE=',round(NSE.val,2),', ',
                        'KGE=',round(KGE.val,2)),
                 cex = 0.9))

integer(0)

data.frame(NSE=NSE.val,KGE=KGE.val) %>% 
  pandoc.table(style="simple",caption="Runoff Model Performance - Daily")

Runoff Model Performance - Daily
NSE	KGE
0.5105	0.6327

4 Evaluation on Project (mean daily flows)

The calibrated and validated GR4J model for the Rio Pulido en Vertedero watershed will be evaluated in the regional watershed where the site is located. The model will be evaluated too in the site-specific watershed that envelopes all the following subcatchments of interest:

-Sub C1 -Sub C2 -Sub C3 -Sub C4 -Sub C5 -Sub C6 -Sub C7 -Sub C8 -Bot Norte -Bot Sur -Planta Procesos -Rajo

4.1 Regional watershed

4.1.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)

#For some reason, the hypsometric curve needs 101 points

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 - Regional watershed",xlab="percentage of area [%]",
     ylab="Elevation[masl]")

4.1.2 Meteorology

#Meteorological information evaluated at the centroid of the regional watershed

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 - regional watershed")

start.eval.date<-as.Date("1980-06-01")

end.eval.date<-as.Date("2022-10-30")

Meteo.site.eval<-meteo_info_z %>% 
  window(start=start.eval.date,end=end.eval.date) %>% 
  na.approx(maxgap = 15)

summary(Meteo.site.eval)

4.1.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=nlayers_project)

#capturing at different elevations  (correcting Meteo.site.eval by elevation for the mean elev of nlayers_project)
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.1.4 Results

#vector of P_eff depending on number of layers in GRXJ model
P_eff_list<-pbapply::pblapply(1:nlayers_project,FUN=function(x){

  y<-if(x<10){
    paste0('OutputsModel.eval$CemaNeigeLayers$Layer0',x,'$PliqAndMelt')
  }else{
    paste0('OutputsModel.eval$CemaNeigeLayers$Layer',x,'$PliqAndMelt')
  }
  
  eval(parse(text=y))
})


  met_gl<-data.frame(
    Date=OutputsModel.eval$DatesR %>% as.Date(),
    P=OutputsModel.eval$Precip,
    P_eff=Reduce(`+`,P_eff_list)/nlayers_project, #averages PliqAndMelt of all layers of GRXJ model
    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 - Regional watershed")
  
  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]")

Monthly Runoff [mm/mon]
Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec	Ann
3.2	2.6	2.5	2.2	2	2.2	2	1.7	1.5	1.5	2.1	2.8	26.4

daily2monthly.V2(site_q_mm_d,FUN=sum) %>% 
  zoo2matrix() %>% 
  boxplot(ylab="mm/mon",xlab=NULL,main="Monthly boxplot for runoff at the regional watershed (mm/mon)")

plot(site_q_mm_d,ylab="mm/d",xlab=NULL,
     main="Daily flow estimated for the regional watershed (mm/d)")

#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,
                                               "Botadero Norte",    2.778,
                                                 "Botadero Sur",    1.337,
                                             "Planta Procesos", 0.928,
                                                        "Rajo", 1.494,
                                                      "Sub C1", 1.044,
                                                      "Sub C2", 0.3627,
                                                      "Sub C3", 0.4471,
                                                      "Sub C4", 1.208,
                                                      "Sub C5", 1.929,
                                                      "Sub C6", 2.103,
                                                      "Sub C7", 0.702,
                                                      "Sub C8", 0.819,
  )


max_flow_fa_rw<-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) %>% 
  bind_rows(data.frame(RP=150,max_daily_m3.s=NA)) %>%  #Bind row with missing RP of interest (150 yrs)
  arrange(RP) #arrange by RP

max_flow_fa_rw$max_daily_m3.s<-na.approx(max_flow_fa_rw$max_daily_m3.s) #replace NA with interpolated value - this line works only because you're trying to find 150 and you have 100 and 200 as data, if it were another number it would not be accurate

cross_join(site_watersheds,
           max_flow_fa_rw) %>%
  mutate(Cuenca=factor(Cuenca,site_watersheds$Cuenca,ordered = T)) %>% 
  mutate(value=area_km2*max_daily_m3.s) %>% 
  dplyr::filter(RP%in%c(50,100,150,500,1000)) %>%
  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 50, 100, 150, 500 and 1000 yrs, based on regional watershed")

Maximum daily average values in m3/s for return periods of 50, 100, 150, 500 and 1000 yrs, based on regional watershed
Cuenca	area_km2	50	100	150	500	1000
Canal Derivación Superior	1.89	0.015	0.018	0.02	0.027	0.031
Canal Derivación Inferior	7.12	0.055	0.067	0.074	0.1	0.117
Qda. N°1(Situación Sin P, Sector Botadero Norte)	3.83	0.029	0.036	0.04	0.054	0.063
Qda. N°2(Situación Sin P, Sector Botadero Norte)	0.63	0.005	0.006	0.007	0.009	0.01
Qda. N°3(Situación Sin P, Sector Botadero Norte)	0.48	0.004	0.005	0.005	0.007	0.008
Qda. N°4(Situación Sin P, Sector Botadero Norte)	0.44	0.003	0.004	0.005	0.006	0.007
Qda. N°5(Situación Sin P, Sector Botadero Norte)	0.74	0.006	0.007	0.008	0.01	0.012
Afluente a Qda. N°5(Situación Sin P, Sector Botadero Norte)	0.29	0.002	0.003	0.003	0.004	0.005
Botadero Norte (Situación Con P)	1.94	0.015	0.018	0.02	0.027	0.032
Botadero Norte	2.778	0.021	0.026	0.029	0.039	0.046
Botadero Sur	1.337	0.01	0.013	0.014	0.019	0.022
Planta Procesos	0.928	0.007	0.009	0.01	0.013	0.015
Rajo	1.494	0.011	0.014	0.015	0.021	0.025
Sub C1	1.044	0.008	0.01	0.011	0.015	0.017
Sub C2	0.363	0.003	0.003	0.004	0.005	0.006
Sub C3	0.447	0.003	0.004	0.005	0.006	0.007
Sub C4	1.208	0.009	0.011	0.012	0.017	0.02
Sub C5	1.929	0.015	0.018	0.02	0.027	0.032
Sub C6	2.103	0.016	0.02	0.022	0.03	0.035
Sub C7	0.702	0.005	0.007	0.007	0.01	0.012
Sub C8	0.819	0.006	0.008	0.008	0.012	0.014

4.2 Site-specific watershed

4.2.1 Hypsometric curve

hypso_ss<-read.csv(here::here('data/csv/site_specific_hyps_curve.csv'))

#For some reason, the hypsometric curve needs 101 points

z_hypso_ss<-approx(x=hypso_ss$AREA.PERCENTAGE,
                 y=hypso_ss$BASE_HEIGHT,
                 xout=seq(0.001,1,length.out=101))

plot(z_hypso_ss,main="Hypsometric curve - Site-specific watershed",xlab="percentage of area [%]",
     ylab="Elevation [masl]")

4.2.2 Meteorology

#Meteorological information evaluated at the centroid of the site_specific watershed (for its mean elevation)

x_ss<--68.89322 #centroid of site-specific watershed
y_ss<--26.01595 #centroid of site-specific watershed
z_ss<-4561.3 #mean elevation of site-specific watershed

meteo_info_z_ss<-meteo_info_fn(x = x_ss,
                              y = y_ss,
                              z = z_ss)

plot(meteo_info_z_ss,main="available information")

start.eval.date<-as.Date("1980-06-01")

end.eval.date<-as.Date("2022-10-30")

Meteo.site.eval_ss<-meteo_info_z_ss %>% 
  window(start=start.eval.date,end=end.eval.date) %>% 
  na.approx(maxgap = 15)

summary(Meteo.site.eval_ss)

4.2.3 Runoff

InputsModel.eval<-CreateInputsModel(
  FUN_MOD = RunModel_CemaNeigeGR4J,
  DatesR = time(Meteo.site.eval_ss) %>% as.POSIXlt(),
  Precip = Meteo.site.eval_ss$prcp%>% coredata,
  PotEvap = Meteo.site.eval_ss$pe %>% coredata,
  TempMax = Meteo.site.eval_ss$tmax%>% coredata,
  TempMin = Meteo.site.eval_ss$tmin%>% coredata,
  TempMean = Meteo.site.eval_ss$tavg%>%coredata,
  ZInputs = z_ss,
  PrecipScale = T,
  HypsoData = na.omit(z_hypso_ss$y),
  NLayers=nlayers_project)

#capturing at different elevations  
meteo_info_lst<-lapply(InputsModel.eval$ZLayers,FUN=function(x){
  meteo_info_fn(x = x_ss,
                y = y_ss,
                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_ss)),
                   grep(end.eval.date, time(Meteo.site.eval_ss)))

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.2.4 Results

#vector of P_eff depending on number of layers in GRXJ model
P_eff_list<-pbapply::pblapply(1:nlayers_project,FUN=function(x){

  y<-if(x<10){
    paste0('OutputsModel.eval$CemaNeigeLayers$Layer0',x,'$PliqAndMelt')
  }else{
    paste0('OutputsModel.eval$CemaNeigeLayers$Layer',x,'$PliqAndMelt')
  }
  
  eval(parse(text=y))
})


  met_gl<-data.frame(
    Date=OutputsModel.eval$DatesR %>% as.Date(),
    P=OutputsModel.eval$Precip,
    P_eff=Reduce(`+`,P_eff_list)/nlayers_project, #averages PliqAndMelt of all layers of GRXJ model
    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 - Site-specific watershed")
  
  print(g10)