Technical note: When best is the enemy of good – critical evaluation of performance criteria in hydrological models
- Guillaume Cinkus 0000-0002-2877-6551
- Naomi Mazzilli 0000-0002-9145-5160
- Hervé Jourde 0000-0001-7124-4879
- Andreas Wunsch 0000-0002-0585-9549
- Tanja Liesch 0000-0001-8648-5333
- Nataša Ravbar 0000-0002-0160-1460
- Zhao Chen 0000-0003-0076-7079
- Nico Goldscheider 0000-0002-8428-5001
This repository contains the following elements:
/script/utils.RCustom functions for KGENP, KGE’’, LME, LCE and NSE/script/run_eval.RFunction to calculate performance criteria values:- KGE and KGE’: hydroGOF R package (Zambrano-Bigiarini, 2020)
- KGENP: supplementary material from Pool et al. (2018)
- DE’: diag-eff Python package (Schwemmle et al., 2021)
- dr: HydroErr Python package (Roberts et al., 2018)
- KGE’’, LME, LCE and NSE: custom code
- KGE_sf: custom code for using KGE with scaling factors
/script/cblind_bw_palette.RColor palette adapted for color vision deficiencies and black/white printing/script/Synthetic time seriesSynthetic time series codesyn_flood_event.RSynthetic time series of a flood event (syn1) or two consecutive flood events (syn2)compare_score_param.RR code to generate Figure 3plot_model_distribution.RR code to generate Figure 4plot_omega_vs_omega2.RR code to generate Figure 5
/script/ANN modelANN model code- Python script to run the ANN model
- Dummy data file to illustrate the structure of input data
/script/Reservoir modelReservoir model code- KarstMod
.propertiesfile - R script to perform the snow routine
- KarstMod
# R packages
library(reticulate)
library(hydroGOF) # for KGE and KGE'
library(dplyr)
library(tidyr)
library(ggplot2)
library(cowplot)
# Python packages
de <- reticulate::import("de") # diaf-eff Python package
hydroErr <- reticulate::import("HydroErr") # HydroEff Python package
# Functions
source("script/utils.R")
source("script/run_eval.R")
source("script/cblind_bw_palette.R")# Import synthetic time series of flood events
## syn1: one event
## syn2: two consecutive events
source("script/Synthetic time series/syn_flood_event.R")
# Define omega bounds
omega_log_min <- -0.36
omega_log_max <- 0.36# Calculate omega bounds (min/max)
omega_min <- 10^(omega_log_min)
omega_max <- 10^(omega_log_max)
# Create min and max transformed time series
model_min <- c(syn2$discharge * omega_min)
model_max <- c(syn2$discharge * omega_max)
# Generate plot
syn2 |>
mutate(model_min, model_max) |>
ggplot(aes(t, discharge)) +
geom_ribbon(aes(ymin = model_min,
ymax = model_max,
fill = "Set of transformed time series")) +
geom_line(aes(color = "Reference time series")) +
scale_color_manual(name = "", values = "black") +
scale_fill_manual(name = "", values = "lightgrey") +
xlab("Time [T]") +
ylab(expression(paste("Discharge [L"^3~T^-1, "]"))) +
coord_cartesian(xlim = c(0, 40)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
theme_bw() +
theme(legend.position = "bottom",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
strip.text = element_text(size = 10),
legend.text = element_text(size = 10))
# Define omega values
omega_1 <- 0.75
omega_2 <- 1.2
# Create two transformed time series
model_bad_bad <- c(syn1$discharge * omega_1, syn1$discharge * omega_2)
model_bad_good <- c(syn1$discharge * omega_1, syn1$discharge)
# Generate plot
syn2 |>
mutate(model_bad_bad, model_bad_good) |>
pivot_longer(c("model_bad_bad", "model_bad_good")) |>
ggplot() +
geom_point(aes(x = t, y = discharge, shape = "Reference time series")) +
geom_line(aes(x = t, y = value, color = "Transformed time series")) +
facet_wrap(~name,
labeller = as_labeller(c("model_bad_bad" = '"Bad-Bad" model (BB model)',
"model_bad_good" = '"Bad-Good" model (BG model)'))) +
coord_cartesian(xlim = c(0, 40)) +
scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
scale_shape_manual(name = "", values = 4) +
scale_color_manual(name = "", values = "#696969") +
xlab("Time [T]") +
ylab(expression(paste("Discharge [L"^3~T^-1, "]"))) +
theme_bw() +
theme(legend.position = "bottom",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
strip.text = element_text(size = 10),
legend.text = element_text(size = 10))
source("script/Synthetic time series/compare_score_param.R")
# Define omega values
omega_1 <- 0.75
omega_2 <- 1.2
# Create two transformed time series
model_bad_bad <- c(syn1$discharge * omega_1, syn1$discharge * omega_2)
model_bad_good <- c(syn1$discharge * omega_1, syn1$discharge)
# Generate plot
compare_score_param(syn2$discharge, model_bad_bad, model_bad_good)
source("script/Synthetic time series/plot_model_distribution.R")
# Vector of criteria to study
criteria <- c("kge", "kge_m", "kge_m2", "kgenp", "de", "lme", "nse", "dr", "lce")
# Generate plot
plot_model_distribution(data = syn1,
perf_c = criteria,
ncol = 3,
omega_log_min = omega_log_min,
omega_log_max = omega_log_max,
step = 0.002)
source("script/Synthetic time series/plot_omega_vs_omega2.R")
# Vector of criteria to study
criteria <- c("kge", "kge_m", "kge_m2", "kgenp", "de", "lme", "nse", "dr", "lce")
# Generate plot
plot_omega_vs_omega2(data = syn1,
perf_c = criteria,
omega_log_min = omega_log_min,
omega_log_max = omega_log_max,
step = 0.002)
Information on the KarstMod platform can be found in Mazzilli et al. (2019b) and the user manual (Mazzilli and Bertin, 2019a). The main workflow is:
- Prepare the input data according to the KarstMod format
- [Optional] Open a KarstMod
.propertiesfile - Import the input data
- Define warm-up/calibration/validation periods
- Define a path for the output directory
- Run calibration
It is possible to modify the model parameters, the objective function,
the number of iterations, the maximum time, and other options. The
Save button allows to save the modifications into a new .properties
file.
Information on 1D-Convolutional Neural Networks can be found in ADD DOI of ANN-lumped study.
Dependencies:
- Python 3.8
- Tensorflow 2.7
- BayesianOptimization 1.2
- Numpy 1.21
- Pandas 1.4
- Scipy 1.7
- Scikit-learn 1.0
- Matplotlib 3.5
Chen, Z., Hartmann, A., Wagener, T., and Goldscheider, N.: Dynamics of water fluxes and storages in an Alpine karst catchment under current and potential future climate conditions, Hydrol. Earth Syst. Sci., 22, 3807–3823, https://doi.org/10.5194/hess-22-3807-2018, 2018.
Hock, R.: A distributed temperature-index ice- and snowmelt model including potential direct solar radiation, J. Glaciol., 45, 101–111, https://doi.org/10.3189/S0022143000003087, 1999.
Mazzilli, N. and Bertin, D.: KarstMod User Guide - version 2.2, hal-01832693, 103927, 2019a.
Mazzilli, N., Guinot, V., Jourde, H., Lecoq, N., Labat, D., Arfib, B., Baudement, C., Danquigny, C., Soglio, L. D., and Bertin, D.: KarstMod: A modelling platform for rainfall - discharge analysis and modelling dedicated to karst systems, Environ. Model. Softw., 122, 103927, https://doi.org/10.1016/j.envsoft.2017.03.015, 2019b.
Pianosi, F., Sarrazin, F., and Wagener, T.: A Matlab toolbox for Global Sensitivity Analysis, Environmental Modelling & Software, 70, 80–85, https://doi.org/10.1016/j.envsoft.2015.04.009, 2015.
Pool, S., Vis, M., and Seibert, J.: Evaluating model performance: towards a non-parametric variant of the Kling-Gupta efficiency, Hydrol. Sci. J., 63, 1941–1953, https://doi.org/10.1080/02626667.2018.1552002, 2018.
Roberts, W., Williams, G., Jackson, E., Nelson, E., Ames, D.: Hydrostats: A Python Package for Characterizing Errors between Observed and Predicted Time Series. Hydrology 5(4) 66, doi:10.3390/hydrology5040066, 2018.
Schwemmle, R., Demand, D., and Weiler, M.: Technical note: Diagnostic efficiency – specific evaluation of model performance, Hydrol. Earth Syst. Sci., 25, 2187–2198, https://doi.org/10.5194/hess-25-2187-2021, 2021.
Zambrano-Bigiarini M.: hydroGOF: Goodness-of-fit functions for comparison of simulated and observed hydrological time series. R package version 0.4-0. https://github.com/hzambran/hydroGOF. DOI:10.5281/zenodo.839854, 2020.