In inrae/airGRiwrm: 'airGR' Integrated Water Resource Management

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(airGRiwrm)

As in vignette V05, this vignette proposes an example of calibration of influenced flow with the Marne reservoir but here with the GR6J model instead of GR4J. It will use influenced observation flows directly measured at gauging stations and flows recorded at reservoir inlets and outlets.

Set the data

Loading naturalized data and influenced flow configuration:

#load("_cache/V01.RData")
load("_cache/V04.RData")

We remove extra items from a complete configuration to keep only the Marne system:

selectedNodes <- c("MARNE_P23", "STDIZ_04", "LOUVE_19", "VITRY_25", "MARNE_P28", "MARNE_R25", "CHALO_21", "MONTR_18", "NOISI_17")
griwrm3 <- griwrm2[griwrm2$id %in% selectedNodes,]
griwrm3$model[!is.na(griwrm3$model)] <- "RunModel_GR6J"
griwrm3[griwrm3$id == "NOISI_17", c("down", "length")] = NA # Downstream station instead of PARIS_05
plot(griwrm3)

We can now generate the new GRiwrmInputsModel object:

library(seinebasin)
data(QOBS)
iEnd <- which(DatesR == as.POSIXct("2008-07-31", tz = "UTC"))

data(Qreservoirs)
QresMarne <- Qreservoirs[1:iEnd, grep("MARNE", colnames(Qreservoirs))]
id_GR_nodes <- griwrm3$id[!is.na(griwrm3$model)]
InputsModel3 <- CreateInputsModel(griwrm3,
                                  DatesR[1:iEnd],
                                  Precip[1:iEnd, id_GR_nodes],
                                  PotEvap[1:iEnd, id_GR_nodes],
                                  QresMarne)

GriwmRunOptions object

We first define the run period:

IndPeriod_Run <- seq.int(
  which(DatesR == (DatesR[1] + 365 * 24 * 60 * 60)), # Set aside warm-up period
  iEnd # Until the end of the time series
)

We define the (optional but recommended) warm up period as a one-year period before the run period:

IndPeriod_WarmUp <- seq.int(1,IndPeriod_Run[1] - 1)

RunOptions <- CreateRunOptions(
  InputsModel3,
  IndPeriod_WarmUp = IndPeriod_WarmUp,
  IndPeriod_Run = IndPeriod_Run
)

InputsCrit object

We define the objective function for the calibration:

InputsCrit <- CreateInputsCrit(
  InputsModel = InputsModel3,
  FUN_CRIT = ErrorCrit_KGE2,
  RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,]
)

GRiwrmCalibOptions object

CalibOptions <- CreateCalibOptions(InputsModel3)
str(CalibOptions)

Calibration

The optimization (i.e. calibration) of parameters can now be performed:

OutputsCalib <- Calibration(InputsModel3, RunOptions, InputsCrit, CalibOptions)

Run model with Michel calibration

Now that the model is calibrated, we can run it with the optimized parameter values:

Param5 <- extractParam(OutputsCalib)

OutputsModels3 <- RunModel(
  InputsModel3,
  RunOptions = RunOptions,
  Param = Param5
)

Comparison with simulated flows

We can compare these simulated flows with influenced discharge measurements:

htmltools::tagList(lapply(
  griwrm3$id[!is.na(griwrm3$model)],
  function(x) {
    Q3 <- Qobs[RunOptions[[1]]$IndPeriod_Run, x]
    iQ3 <- which(!is.na(Q3))
    IndPeriod_Obs <- iQ3[1]:tail(iQ3, 1)
    OutputsModels <- ReduceOutputsModel(OutputsModels3[[x]], IndPeriod_Obs)
    plot(OutputsModels, Qobs = Q3[IndPeriod_Obs], main = x)
  }
))