airGRiwrm automates the execution of airGR semi-distributive models. The steps for running or calibrating the model are the same as the ones of ‘airGR’.

Load library

Preparation of functions inputs

To run a model, as for airGR, the functions of the airGRiwrm package (e.g. the models, calibration and criteria calculation functions) require data and options with specific formats.

To facilitate the use of the package, there are several functions dedicated to the creation of these objects:

  • CreateInputsModel(): prepares the inputs for the different hydrological models (times series of dates, precipitation, observed discharge, etc.)
  • CreateRunOptions(): prepares the options for the hydrological model run (warm up period, calibration period, etc.)
  • CreateInputsCrit(): prepares the options in order to compute the efficiency criterion (choice of the criterion, choice of the transformation on discharge: “log”, “sqrt”, etc.)
  • CreateCalibOptions(): prepares the options for the hydrological model calibration algorithm (choice of parameters to optimize, predefined values for uncalibrated parameters, etc.)

GRiwrmInputsModel object

The production method of the GRiwrmInputsModel object is detailed in the vignette “V01_Structure_SD_model” of the package. The following code chunk resumes all the steps of this vignette:

data(Severn)
nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
nodes$distance_downstream <- nodes$distance_downstream
nodes$model <- "RunModel_GR4J"
griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
BasinsObs <- Severn$BasinsObs
DatesR <- BasinsObs[[1]]$DatesR
PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
## CreateInputsModel.GRiwrm: Treating sub-basin 54095...
## CreateInputsModel.GRiwrm: Treating sub-basin 54002...
## CreateInputsModel.GRiwrm: Treating sub-basin 54029...
## CreateInputsModel.GRiwrm: Treating sub-basin 54001...
## CreateInputsModel.GRiwrm: Treating sub-basin 54032...
## CreateInputsModel.GRiwrm: Treating sub-basin 54057...
str(InputsModel)
## List of 6
##  $ 54095:List of 7
##   ..$ DatesR    : POSIXlt[1:11536], format: "1984-02-29 23:00:00" "1984-03-01 23:00:00" ...
##   ..$ Precip    : num [1:11536] 4.01 0.57 2 0.37 0.01 0.3 0.12 0 0.12 0.07 ...
##   ..$ PotEvap   : num [1:11536] 0.59 1.64 1.42 0.31 0.59 0.73 0.59 0.66 0.57 0.61 ...
##   ..$ id        : chr "54095"
##   ..$ down      : chr "54001"
##   ..$ BasinAreas: num 3723
##   ..$ FUN_MOD   : chr "RunModel_GR4J"
##   ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR"
##  $ 54002:List of 7
##   ..$ DatesR    : POSIXlt[1:11536], format: "1984-02-29 23:00:00" "1984-03-01 23:00:00" ...
##   ..$ Precip    : num [1:11536] 1.58 0.47 1.35 1.92 0.06 0 0.01 0 0.08 0.34 ...
##   ..$ PotEvap   : num [1:11536] 0.61 1.7 1.61 0.3 0.44 0.69 0.52 0.71 0.73 0.57 ...
##   ..$ id        : chr "54002"
##   ..$ down      : chr "54057"
##   ..$ BasinAreas: num 2208
##   ..$ FUN_MOD   : chr "RunModel_GR4J"
##   ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR"
##  $ 54029:List of 7
##   ..$ DatesR    : POSIXlt[1:11536], format: "1984-02-29 23:00:00" "1984-03-01 23:00:00" ...
##   ..$ Precip    : num [1:11536] 2.38 0.33 2.16 0.38 0.01 0.12 0.08 0.05 0.05 0.29 ...
##   ..$ PotEvap   : num [1:11536] 0.58 1.64 1.49 0.23 0.56 0.72 0.63 0.72 0.62 0.64 ...
##   ..$ id        : chr "54029"
##   ..$ down      : chr "54032"
##   ..$ BasinAreas: num 1484
##   ..$ FUN_MOD   : chr "RunModel_GR4J"
##   ..- attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR"
##  $ 54001:List of 11
##   ..$ DatesR          : POSIXlt[1:11536], format: "1984-02-29 23:00:00" "1984-03-01 23:00:00" ...
##   ..$ Precip          : num [1:11536] 1.3 0.427 2.642 0.441 0.01 ...
##   ..$ PotEvap         : num [1:11536] 0.59 1.711 1.563 0.239 0.519 ...
##   ..$ Qupstream       : num [1:11536, 1] 3350412 3350412 3499319 3238732 3201505 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr "54095"
##   ..$ LengthHydro     : Named num 42
##   .. ..- attr(*, "names")= chr "54095"
##   ..$ BasinAreas      : Named num [1:2] 3723 607
##   .. ..- attr(*, "names")= chr [1:2] "54095" "54001"
##   ..$ id              : chr "54001"
##   ..$ down            : chr "54032"
##   ..$ UpstreamNodes   : chr "54095"
##   ..$ UpstreamIsRunoff: logi TRUE
##   ..$ FUN_MOD         : chr "RunModel_GR4J"
##   ..- attr(*, "class")= chr [1:4] "InputsModel" "daily" "GR" "SD"
##  $ 54032:List of 11
##   ..$ DatesR          : POSIXlt[1:11536], format: "1984-02-29 23:00:00" "1984-03-01 23:00:00" ...
##   ..$ Precip          : num [1:11536] 1.737 0.469 2.187 1.229 0.01 ...
##   ..$ PotEvap         : num [1:11536] 0.604 1.599 1.565 0.203 0.543 ...
##   ..$ Qupstream       : num [1:11536, 1:2] 3334023 3334023 3290724 3204126 3117528 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "54001" "54029"
##   ..$ LengthHydro     : Named num [1:2] 45 32
##   .. ..- attr(*, "names")= chr [1:2] "54001" "54029"
##   ..$ BasinAreas      : Named num [1:3] 4330 1484 1051
##   .. ..- attr(*, "names")= chr [1:3] "54001" "54029" "54032"
##   ..$ id              : chr "54032"
##   ..$ down            : chr "54057"
##   ..$ UpstreamNodes   : chr [1:2] "54001" "54029"
##   ..$ UpstreamIsRunoff: logi [1:2] TRUE TRUE
##   ..$ FUN_MOD         : chr "RunModel_GR4J"
##   ..- attr(*, "class")= chr [1:4] "InputsModel" "daily" "GR" "SD"
##  $ 54057:List of 11
##   ..$ DatesR          : POSIXlt[1:11536], format: "1984-02-29 23:00:00" "1984-03-01 23:00:00" ...
##   ..$ Precip          : num [1:11536] 1.523 0.179 1.536 1.545 0 ...
##   ..$ PotEvap         : num [1:11536] 0.536 1.599 1.576 0.31 0.437 ...
##   ..$ Qupstream       : num [1:11536, 1:2] 5766499 5697850 5560553 5423255 5354606 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "54032" "54002"
##   ..$ LengthHydro     : Named num [1:2] 15 43
##   .. ..- attr(*, "names")= chr [1:2] "54032" "54002"
##   ..$ BasinAreas      : Named num [1:3] 6865 2208 813
##   .. ..- attr(*, "names")= chr [1:3] "54032" "54002" "54057"
##   ..$ id              : chr "54057"
##   ..$ down            : chr NA
##   ..$ UpstreamNodes   : chr [1:2] "54032" "54002"
##   ..$ UpstreamIsRunoff: logi [1:2] TRUE TRUE
##   ..$ FUN_MOD         : chr "RunModel_GR4J"
##   ..- attr(*, "class")= chr [1:4] "InputsModel" "daily" "GR" "SD"
##  - attr(*, "class")= chr [1:2] "GRiwrmInputsModel" "list"
##  - attr(*, "GRiwrm")=Classes 'GRiwrm' and 'data.frame':  6 obs. of  5 variables:
##   ..$ id    : chr [1:6] "54057" "54032" "54001" "54095" ...
##   ..$ down  : chr [1:6] NA "54057" "54032" "54001" ...
##   ..$ length: num [1:6] NA 15 45 42 43 32
##   ..$ model : chr [1:6] "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" "RunModel_GR4J" ...
##   ..$ area  : num [1:6] 9885 6865 4330 3723 2208 ...
##  - attr(*, "TimeStep")= num 86400

GriwmRunOptions object

The CreateRunOptions() function allows to prepare the options required for the RunModel() function.

The user must at least define the following arguments:

  • InputsModel: the associated input data
  • IndPeriod_Run: the period on which the model is run

Here below, we start the period to run one year after the beginning of the time series.

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

So the warmup period could be defined as:

IndPeriod_WarmUp = seq(1,IndPeriod_Run[1]-1)

Parameters of CreateRunOptions for airGRiwrm are the same as the function in airGR and are copied for each node running a rainfall-runoff model.

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

GRiwrmInputsCrit object

The CreateInputsCrit() function allows to prepare the input in order to calculate a criterion. It takes the following arguments:

  • InputsModel: the inputs of the GRiwrm network previously prepared by the CreateInputsModel() function
  • FUN_CRIT: the name of the error criterion function (see the available functions description in the airGR package)
  • RunOptions: the options of the GRiwrm network previously prepared by the CreateRunOptions() function
  • Qobs: the observed variable time serie (e.g. the discharge expressed in mm/time step)
InputsCrit <- CreateInputsCrit(
  InputsModel = InputsModel,
  FUN_CRIT = airGR::ErrorCrit_KGE2,
  RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,]
)
str(InputsCrit)
## List of 6
##  $ 54095:List of 8
##   ..$ FUN_CRIT: chr "ErrorCrit_KGE2"
##   ..$ Obs     : num [1:11171] 1.19 1.17 1.25 1.51 2.05 1.61 1.34 1.34 1.17 1.06 ...
##   ..$ VarObs  : chr "Q"
##   ..$ BoolCrit: logi [1:11171] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   ..$ idLayer : logi NA
##   ..$ transfo : chr ""
##   ..$ epsilon : NULL
##   ..$ Weights : NULL
##   ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit"
##  $ 54002:List of 8
##   ..$ FUN_CRIT: chr "ErrorCrit_KGE2"
##   ..$ Obs     : num [1:11171] 1.09 1.13 1.07 1.04 0.82 0.68 0.87 0.88 0.77 0.68 ...
##   ..$ VarObs  : chr "Q"
##   ..$ BoolCrit: logi [1:11171] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   ..$ idLayer : logi NA
##   ..$ transfo : chr ""
##   ..$ epsilon : NULL
##   ..$ Weights : NULL
##   ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit"
##  $ 54029:List of 8
##   ..$ FUN_CRIT: chr "ErrorCrit_KGE2"
##   ..$ Obs     : num [1:11171] 1.37 1.31 1.52 1.6 1.29 1.14 1.45 1.47 1.23 1.12 ...
##   ..$ VarObs  : chr "Q"
##   ..$ BoolCrit: logi [1:11171] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   ..$ idLayer : logi NA
##   ..$ transfo : chr ""
##   ..$ epsilon : NULL
##   ..$ Weights : NULL
##   ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit"
##  $ 54001:List of 8
##   ..$ FUN_CRIT: chr "ErrorCrit_KGE2"
##   ..$ Obs     : num [1:11171] 1.1 1.05 1.19 1.29 1.77 1.49 1.31 1.27 1.1 0.97 ...
##   ..$ VarObs  : chr "Q"
##   ..$ BoolCrit: logi [1:11171] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   ..$ idLayer : logi NA
##   ..$ transfo : chr ""
##   ..$ epsilon : NULL
##   ..$ Weights : NULL
##   ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit"
##  $ 54032:List of 8
##   ..$ FUN_CRIT: chr "ErrorCrit_KGE2"
##   ..$ Obs     : num [1:11171] 1.19 1.19 1.23 1.35 1.43 1.37 1.32 1.35 1.21 1.07 ...
##   ..$ VarObs  : chr "Q"
##   ..$ BoolCrit: logi [1:11171] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   ..$ idLayer : logi NA
##   ..$ transfo : chr ""
##   ..$ epsilon : NULL
##   ..$ Weights : NULL
##   ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit"
##  $ 54057:List of 8
##   ..$ FUN_CRIT: chr "ErrorCrit_KGE2"
##   ..$ Obs     : num [1:11171] 1.06 1.1 1.09 1.22 1.23 1.21 1.18 1.25 1.13 0.98 ...
##   ..$ VarObs  : chr "Q"
##   ..$ BoolCrit: logi [1:11171] TRUE TRUE TRUE TRUE TRUE TRUE ...
##   ..$ idLayer : logi NA
##   ..$ transfo : chr ""
##   ..$ epsilon : NULL
##   ..$ Weights : NULL
##   ..- attr(*, "class")= chr [1:2] "Single" "InputsCrit"
##  - attr(*, "class")= chr [1:2] "list" "GRiwrmInputsCrit"

GRiwrmCalibOptions object

Before using the automatic calibration tool, the user needs to prepare the calibration options with the CreateCalibOptions() function. The GRiwrmInputsModel argument contains all the necessary information:

CalibOptions <- CreateCalibOptions(InputsModel)

Calibration

The airGR calibration process is applied on each node of the GRiwrm network from upstream nodes to downstream nodes.

OutputsCalib <- suppressWarnings(
  Calibration(InputsModel, RunOptions, InputsCrit, CalibOptions))
## Calibration.GRiwrmInputsModel: Treating sub-basin 54095...
## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
##   Screening completed (81 runs)
##       Param =  247.151,   -0.020,   83.096,    2.384
##       Crit. KGE2[Q]      = 0.9413
## Steepest-descent local search in progress
##   Calibration completed (26 iterations, 274 runs)
##       Param =  252.524,    0.031,   55.098,    3.293
##       Crit. KGE2[Q]      = 0.9587
## Calibration.GRiwrmInputsModel: Treating sub-basin 54002...
## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
##   Screening completed (81 runs)
##       Param =  247.151,   -0.020,   20.697,    1.944
##       Crit. KGE2[Q]      = 0.9396
## Steepest-descent local search in progress
##   Calibration completed (16 iterations, 193 runs)
##       Param =  223.632,   -0.020,   18.916,    2.198
##       Crit. KGE2[Q]      = 0.9517
## Calibration.GRiwrmInputsModel: Treating sub-basin 54029...
## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
##   Screening completed (81 runs)
##       Param =  247.151,   -0.020,   42.098,    1.944
##       Crit. KGE2[Q]      = 0.9455
## Steepest-descent local search in progress
##   Calibration completed (41 iterations, 409 runs)
##       Param =  220.601,   -0.084,   37.744,    2.116
##       Crit. KGE2[Q]      = 0.9615
## Calibration.GRiwrmInputsModel: Treating sub-basin 54001...
## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
##   Screening completed (243 runs)
##       Param =   11.250,  432.681,   -2.376,   20.697,    1.417
##       Crit. KGE2[Q]      = 0.9262
## Steepest-descent local search in progress
##   Calibration completed (52 iterations, 718 runs)
##       Param =    2.348, 4402.818, -10903.649,   35.163,   17.500
##       Crit. KGE2[Q]      = 0.9541
## Calibration.GRiwrmInputsModel: Treating sub-basin 54032...
## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
##   Screening completed (243 runs)
##       Param =   11.250,  432.681,   -0.020,   83.096,    1.417
##       Crit. KGE2[Q]      = 0.9493
## Steepest-descent local search in progress
##   Calibration completed (76 iterations, 995 runs)
##       Param =    1.105, 1925.349,   -0.145,    6.321,    1.817
##       Crit. KGE2[Q]      = 0.9633
## Calibration.GRiwrmInputsModel: Treating sub-basin 54057...
## Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
##   Screening completed (243 runs)
##       Param =   11.250,  169.017,   -0.649,   83.096,    2.384
##       Crit. KGE2[Q]      = 0.9552
## Steepest-descent local search in progress
##   Calibration completed (218 iterations, 2522 runs)
##       Param =    0.777,  146.867,   -0.101,    0.089,    8.615
##       Crit. KGE2[Q]      = 0.9616
ParamMichel <- sapply(griwrm$id, function(x) {OutputsCalib[[x]]$Param})

Run model with Michel calibration

OutputsModels <- RunModel(
  InputsModel,
  RunOptions = RunOptions,
  Param = ParamMichel
)
## RunModel.GRiwrmInputsModel: Treating sub-basin 54095...
## RunModel.GRiwrmInputsModel: Treating sub-basin 54002...
## RunModel.GRiwrmInputsModel: Treating sub-basin 54029...
## RunModel.GRiwrmInputsModel: Treating sub-basin 54001...
## RunModel.GRiwrmInputsModel: Treating sub-basin 54032...
## RunModel.GRiwrmInputsModel: Treating sub-basin 54057...

Plot the result for each basin

htmltools::tagList(lapply(
  names(OutputsModels),
  function(x) {
    plot(OutputsModels[[x]], Qobs = Qobs[IndPeriod_Run,x] , main = x)
  }
))

The resulting flow of each node in m3/s is directly available and can be plotted with these commands:

Qm3s <- attr(OutputsModels, "Qm3s")
plot(Qm3s[1:150,])