In emf-creaf/medfateland: Mediterranean Landscape Simulation

library(kableExtra)
library(sf)
library(knitr)
library(medfate)
library(ggplot2)
library(cowplot)
library(hydroGOF)
evalDIR = "~/OneDrive/mcaceres_work/model_development/medfate_evaluation/WatershedLevelEvaluation/"

medfateland_versions <- list.files(file.path(evalDIR, "data", "output"))
medfateland_version_vec <- character(0)
model_vec <- character(0)
watershed_vec <- character(0)
for(medfateland_version in medfateland_versions) {
  models <- list.files(file.path(evalDIR, "data", "output", medfateland_version))
  for(model in models) {
    watersheds <- list.files(file.path(evalDIR, "data", "output", medfateland_version, model))
    for(watershed in watersheds) {
      medfateland_version_vec <- c(medfateland_version_vec, medfateland_version)
      model_vec <- c(model_vec, models)
      watershed_vec <- c(watershed_vec, watershed)
    }
  }
}
df <- tibble::tibble(medfateland = medfateland_version_vec, 
                     medfate = NA,
                     model = model_vec, watershed = watershed_vec) 
df$before <- vector("list", nrow(df))
df$after <- vector("list", nrow(df))
df$Q_daily <- vector("list", nrow(df))
df$Q_monthly <- vector("list", nrow(df))
for(r in 1:nrow(df)) {
  df$before[[r]] <- yaml::read_yaml(file.path(evalDIR, "data", "output", df$medfateland[r], df$model[r], df$watershed[r], "uncalibrated.yaml"))
  df$after[[r]] <- yaml::read_yaml(file.path(evalDIR, "data", "output", df$medfateland[r], df$model[r], df$watershed[r], "calibrated.yaml"))
  df$Q_daily[[r]] <- read.csv2(file.path(evalDIR, "data", "output", df$medfateland[r], df$model[r], df$watershed[r], "Q_daily.csv"))
  df$Q_monthly[[r]] <- read.csv2(file.path(evalDIR, "data", "output", df$medfateland[r], df$model[r], df$watershed[r], "Q_monthly.csv"))
}

Introduction

Goal

The aim of this article is to provide an assessment of the performance of spwb_land() for the prediction of watershed outflow. To this aim, we simulate hydrological processes in a set of benchmark watersheds and compare the model predictions of watershed outflow against measurements obtained using river gauges at watershed outlets.

Simulation procedure

For each watershed, the following procedure has been conducted:

Initial warm-up simulation for a specified number of years
Simulation for the period with observed data before calibrating watershed parameters
Manual calibration of watershed parameters (to be replaced with automatized calibration)
Final simulation for the period with observed data after calibrating watershed parameters.

Goodness-of-fit statistics

The following goodness of fit statistics are calculated using package hydroGOF:

Nash-Sutcliffe Efficiency (NSE): This coefficient is sensitive to extreme values and might yield sub-optimal results when the dataset contains large outliers.
Kling–Gupta efficiency (KGE): Provides a decomposition of the Nash-Sutcliffe efficiency, which facilitates the analysis of the importance of different components (bias, correlation and variability).
Index of agreement (d): Initially proposed by Willmott (1981) to overcome the drawbacks of the R2, such as the differences in observed and predicted means and variances (Legates and McCabe, 1999). d is also dimensionless and bounded between 0 and 1 and can be interpreted similarly to R2.
Volumetric efficiency index (VE): Originally proposed by Criss and Winston (2008) to circumvent some of the NSE flaws. VE values are also bounded [0, 1] and represent the fraction of water delivered at the proper time.
Root mean squared error (RMSE): The usual estimation of average model error (i.e. the square root of mean squared errors).

Hydrological analysis

Watershed (TETIS) parameters

The following table contains the set of TETIS parameters employed in spwb_land() simulations on all watersheds, before and after calibration:

df_params <- df |> 
  dplyr::select(-Q_daily, -Q_monthly) |>
  tidyr::pivot_longer(cols = c(before, after), names_to = "Calibration", values_to ="Params")
df_params$R_localflow <- NA
df_params$R_interflow <- NA
df_params$R_baseflow <- NA
df_params$n_interflow <- NA
df_params$n_baseflow <- NA
df_params$num_daily_substeps <- NA
df_params$rock_max_infiltration <- NA
df_params$deep_aquifer_loss <- NA
for(i in 1:nrow(df_params)) {
  df_params$medfateland[i] <- df_params$Params[[i]]$medfateland_version
  df_params$medfate[i] <- df_params$Params[[i]]$medfate_version
  df_params$R_localflow[i] <- df_params$Params[[i]]$tetis_parameters$R_localflow
  df_params$R_interflow[i] <- df_params$Params[[i]]$tetis_parameters$R_interflow
  df_params$R_baseflow[i] <- df_params$Params[[i]]$tetis_parameters$R_baseflow
  df_params$n_interflow[i] <- df_params$Params[[i]]$tetis_parameters$n_interflow
  df_params$n_baseflow[i] <- df_params$Params[[i]]$tetis_parameters$n_baseflow
  df_params$num_daily_substeps[i] <- df_params$Params[[i]]$tetis_parameters$num_daily_substeps
  df_params$rock_max_infiltration[i] <- df_params$Params[[i]]$tetis_parameters$rock_max_infiltration
  df_params$deep_aquifer_loss[i] <- df_params$Params[[i]]$tetis_parameters$deep_aquifer_loss
}
df_params |>
  dplyr::select(-Params) |>
  kableExtra::kbl() |>
  kableExtra::kable_styling()

Evaluation results

for(i in 1:nrow(df)) {
   medfateland_version <- df$medfateland[i] 
   model <- df$model[i]
   watershed <- df$watershed[i]
   Q_daily <- df$Q_daily[[i]] |>
    dplyr::mutate(Date = as.Date(Date))
  Q_monthly <- df$Q_monthly[[i]] |>
    dplyr::mutate(Date = as.Date(Date))
  res <- knitr::knit_child('_WLEvaluation_child.Rmd', quiet = TRUE)
  cat(res, sep = '\n')
}