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 used influenced observation flows directly measured at gauging stations and flows recorded at reservoir inlets and outlet.

Set the data

Loading naturalised data and influenced flow configuration:

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

Remove extra items from global configuration

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)
#> NULL

Generate GRiwrmInputsModel object

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

data(Qreservoirs)
Qobs3 <- cbind(Qobs, Qreservoirs[1:iEnd,])
InputsModel3 <- CreateInputsModel(griwrm3, DatesR[1:iEnd], Precip[1:iEnd,], PotEvap[1:iEnd,], Qobs3)
#> CreateInputsModel.GRiwrm: Treating sub-basin STDIZ_04...
#> CreateInputsModel.GRiwrm: Treating sub-basin MONTR_18...
#> CreateInputsModel.GRiwrm: Treating sub-basin LOUVE_19...
#> CreateInputsModel.GRiwrm: Treating sub-basin VITRY_25...
#> CreateInputsModel.GRiwrm: Treating sub-basin CHALO_21...
#> CreateInputsModel.GRiwrm: Treating sub-basin NOISI_17...

GriwmRunOptions object

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
)

The warmup period could also be defined as is:

IndPeriod_WarmUp = seq.int(1,IndPeriod_Run[1]-1)
RunOptions <- CreateRunOptions(
  InputsModel = InputsModel3,
  IndPeriod_WarmUp = IndPeriod_WarmUp,
  IndPeriod_Run = IndPeriod_Run
)

InputsCrit object

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

GRiwrmCalibOptions object

CalibOptions <- CreateCalibOptions(InputsModel3)
str(CalibOptions)
#> List of 6
#>  $ STDIZ_04:List of 4
#>   ..$ FixedParam       : logi [1:7] NA NA NA NA NA NA ...
#>   ..$ SearchRanges     : num [1:2, 1:7] 1.00e-02 2.00e+01 4.59e-05 2.18e+04 -1.09e+04 ...
#>   ..$ FUN_TRANSFO      :function (ParamIn, Direction)  
#>   ..$ StartParamDistrib: num [1:3, 1:7] 11.2 12.5 15 36.6 49.4 ...
#>   ..- attr(*, "class")= chr [1:4] "CalibOptions" "GR6J" "SD" "HBAN"
#>  $ MONTR_18:List of 4
#>   ..$ FixedParam       : logi [1:6] NA NA NA NA NA NA
#>   ..$ SearchRanges     : num [1:2, 1:6] 4.59e-05 2.18e+04 -1.09e+04 1.09e+04 4.59e-05 ...
#>   ..$ FUN_TRANSFO      :function (ParamIn, Direction)  
#>   ..$ StartParamDistrib: num [1:3, 1:6] 36.598 49.402 90.017 -1.175 -0.521 ...
#>   ..- attr(*, "class")= chr [1:3] "CalibOptions" "GR6J" "HBAN"
#>  $ LOUVE_19:List of 4
#>   ..$ FixedParam       : logi [1:6] NA NA NA NA NA NA
#>   ..$ SearchRanges     : num [1:2, 1:6] 4.59e-05 2.18e+04 -1.09e+04 1.09e+04 4.59e-05 ...
#>   ..$ FUN_TRANSFO      :function (ParamIn, Direction)  
#>   ..$ StartParamDistrib: num [1:3, 1:6] 36.598 49.402 90.017 -1.175 -0.521 ...
#>   ..- attr(*, "class")= chr [1:3] "CalibOptions" "GR6J" "HBAN"
#>  $ VITRY_25:List of 4
#>   ..$ FixedParam       : logi [1:6] NA NA NA NA NA NA
#>   ..$ SearchRanges     : num [1:2, 1:6] 4.59e-05 2.18e+04 -1.09e+04 1.09e+04 4.59e-05 ...
#>   ..$ FUN_TRANSFO      :function (ParamIn, Direction)  
#>   ..$ StartParamDistrib: num [1:3, 1:6] 36.598 49.402 90.017 -1.175 -0.521 ...
#>   ..- attr(*, "class")= chr [1:3] "CalibOptions" "GR6J" "HBAN"
#>  $ CHALO_21:List of 4
#>   ..$ FixedParam       : logi [1:7] NA NA NA NA NA NA ...
#>   ..$ SearchRanges     : num [1:2, 1:7] 1.00e-02 2.00e+01 4.59e-05 2.18e+04 -1.09e+04 ...
#>   ..$ FUN_TRANSFO      :function (ParamIn, Direction)  
#>   ..$ StartParamDistrib: num [1:3, 1:7] 11.2 12.5 15 36.6 49.4 ...
#>   ..- attr(*, "class")= chr [1:4] "CalibOptions" "GR6J" "SD" "HBAN"
#>  $ NOISI_17:List of 4
#>   ..$ FixedParam       : logi [1:7] NA NA NA NA NA NA ...
#>   ..$ SearchRanges     : num [1:2, 1:7] 1.00e-02 2.00e+01 4.59e-05 2.18e+04 -1.09e+04 ...
#>   ..$ FUN_TRANSFO      :function (ParamIn, Direction)  
#>   ..$ StartParamDistrib: num [1:3, 1:7] 11.2 12.5 15 36.6 49.4 ...
#>   ..- attr(*, "class")= chr [1:4] "CalibOptions" "GR6J" "SD" "HBAN"

Calibration

OutputsCalib <- Calibration(InputsModel3, RunOptions, InputsCrit, CalibOptions)
#> Calibration.GRiwrmInputsModel: Treating sub-basin STDIZ_04...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (2187 runs)
#>       Param =   15.000,   36.598,   -0.521,   60.340,    2.345,    0.220,   20.086
#>       Crit. KGE2[Q]      = 0.7948
#> Steepest-descent local search in progress
#>   Calibration completed (184 iterations, 4770 runs)
#>       Param =   19.990,  161.075,   -0.105,   49.258,    3.700,    0.302,    3.513
#>       Crit. KGE2[Q]      = 0.9180
#> Calibration.GRiwrmInputsModel: Treating sub-basin MONTR_18...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (729 runs)
#>       Param =   90.017,   -0.521,   60.340,    2.345,    0.220,   20.086
#>       Crit. KGE2[Q]      = 0.7904
#> Steepest-descent local search in progress
#>   Calibration completed (79 iterations, 1699 runs)
#>       Param =  175.611,   -0.318,   46.474,    2.412,    0.272,    6.006
#>       Crit. KGE2[Q]      = 0.8382
#> Calibration.GRiwrmInputsModel: Treating sub-basin LOUVE_19...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (729 runs)
#>       Param =   49.402,   -0.521,   60.340,    2.345,    0.020,   20.086
#>       Crit. KGE2[Q]      = 0.9127
#> Steepest-descent local search in progress
#>   Calibration completed (33 iterations, 1099 runs)
#>       Param =   58.518,   -0.521,   94.421,    2.221,    0.026,   14.515
#>       Crit. KGE2[Q]      = 0.9290
#> Calibration.GRiwrmInputsModel: Treating sub-basin VITRY_25...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (729 runs)
#>       Param =   36.598,   -0.521,  148.413,    2.345,    0.020,   20.086
#>       Crit. KGE2[Q]      = 0.9197
#> Steepest-descent local search in progress
#>   Calibration completed (56 iterations, 1394 runs)
#>       Param =   33.665,   -0.521,  163.534,    3.691,    0.043,   16.471
#>       Crit. KGE2[Q]      = 0.9473
#> Calibration.GRiwrmInputsModel: Treating sub-basin CHALO_21...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (2187 runs)
#>       Param =   11.250,   90.017,   -1.175,  148.413,    1.369,    0.020,  148.413
#>       Crit. KGE2[Q]      = 0.9016
#> Steepest-descent local search in progress
#>   Calibration completed (135 iterations, 4052 runs)
#>       Param =    0.410,  584.058,   -1.160,  487.846,    2.374,    0.022,    3.896
#>       Crit. KGE2[Q]      = 0.9553
#> Calibration.GRiwrmInputsModel: Treating sub-basin NOISI_17...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (2187 runs)
#>       Param =   11.250,   90.017,   -1.175,  148.413,    2.345,    0.220,  148.413
#>       Crit. KGE2[Q]      = 0.8327
#> Steepest-descent local search in progress
#>   Calibration completed (57 iterations, 2955 runs)
#>       Param =    1.680,  962.949,   -1.175,  181.272,    2.159,    0.216,   21.542
#>       Crit. KGE2[Q]      = 0.9482

Run model with Michel calibration

Param5 <- sapply(griwrm3$id, function(x) {OutputsCalib[[x]]$Param})

OutputsModels3 <- RunModel(
  InputsModel3,
  RunOptions = RunOptions,
  Param = Param5
)
#> RunModel.GRiwrmInputsModel: Treating sub-basin STDIZ_04...
#> Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 120
#> time steps with negative flow, set to zero.
#> RunModel.GRiwrmInputsModel: Treating sub-basin MONTR_18...
#> RunModel.GRiwrmInputsModel: Treating sub-basin LOUVE_19...
#> RunModel.GRiwrmInputsModel: Treating sub-basin VITRY_25...
#> RunModel.GRiwrmInputsModel: Treating sub-basin CHALO_21...
#> Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 14
#> time steps with negative flow, set to zero.
#> RunModel.GRiwrmInputsModel: Treating sub-basin NOISI_17...
#> Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 5 time
#> steps with negative flow, set to zero.

Comparison with simulated flows

htmltools::tagList(lapply(
  griwrm3$id[!is.na(griwrm3$model)],
  function(x) {
    Q3 <- Qobs3[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)
  }
))
#> Warning in plot.OutputsModel(OutputsModels, Qobs = Q3[IndPeriod_Obs], main = x):
#> zeroes detected in 'Qsim': some plots in the log space will not be created using
#> all time-steps