Presentation of the study case

We can see in vignette ‘V02_Calibration_SD_model’ that calibration result for the flow simulated on the Avon at Evesham (Gauging station ‘54002’) and on the Severn at Buildwas (Gauging station ‘54095’) is not satisfactory. These upper basins are heavily influenced by impoundments and inter-basin transfers (Higgs and Petts 1988).

To cope with this influenced flow, we won’t try to simulate it with a hydrological model and we will directly inject in the model the observed flow at these nodes.

Convert a gauging station into a release spot

Modifying the GRiwrm object

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

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"))

To notify the SD model that the provided flow on a node should be directly used instead of an hydrological model, you only need to declare its model as NA:

griwrmV03 <- griwrm
griwrmV03$model[griwrm$id == "54002"] <- NA
griwrmV03$model[griwrm$id == "54095"] <- NA
griwrmV03
#>      id  down length         model    area
#> 1 54057  <NA>     NA RunModel_GR4J 9885.46
#> 2 54032 54057     15 RunModel_GR4J 6864.88
#> 3 54001 54032     45 RunModel_GR4J 4329.90
#> 4 54095 54001     42          <NA> 3722.68
#> 5 54002 54057     43          <NA> 2207.95
#> 6 54029 54032     32 RunModel_GR4J 1483.65

Here, we keep the area of this basin which means that the discharge will be provided in mm per time step. If the discharge is provided in m3/s, then the area should be set to NA and downstream basin areas should be modified subsequently.

The diagram of the network structure is represented below with:

  • in blue the upstream nodes with a GR4J model
  • in green the intermediate nodes with an SD (GR4J + LAG) model
  • in red the node with direct flow injection (no hydrological model)
plot(griwrmV03)
#> NULL

Generate the GRiwrmInputsModel object

The formatting of the input data is described in the vignette “V01_Structure_SD_model.” the following code chunk resumes this formatting procedure:

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)

The GRiwrmInputsModel object can be generated taking into account the new GRiwrm object:

IM_OL <- CreateInputsModel(griwrmV03, DatesR, Precip, PotEvap, Qobs)
#> CreateInputsModel.GRiwrm: Treating sub-basin 54001...
#> CreateInputsModel.GRiwrm: Treating sub-basin 54029...
#> CreateInputsModel.GRiwrm: Treating sub-basin 54032...
#> CreateInputsModel.GRiwrm: Treating sub-basin 54057...

Calibration of the new model

Calibration options and criteria procedures are detailed in vignette “V02_Calibration_SD_model.” The following code chunk resume this procedure:

IndPeriod_Run <- seq(
  which(DatesR == (DatesR[1] + 365*24*60*60)), # Set aside warm-up period
  length(DatesR) # Until the end of the time series
)
IndPeriod_WarmUp = seq(1,IndPeriod_Run[1]-1)
RunOptions <- CreateRunOptions(IM_OL,
                               IndPeriod_WarmUp = IndPeriod_WarmUp,
                               IndPeriod_Run = IndPeriod_Run)
InputsCrit <- CreateInputsCrit(IM_OL,
                               FUN_CRIT = airGR::ErrorCrit_KGE2,
                               RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,])
CalibOptions <- CreateCalibOptions(IM_OL)

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

OC_OL <- suppressWarnings(
  Calibration(IM_OL, RunOptions, InputsCrit, CalibOptions))
#> 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,    2.384
#>       Crit. KGE2[Q]      = 0.9391
#> Steepest-descent local search in progress
#>   Calibration completed (54 iterations, 737 runs)
#>       Param =    2.116, 4865.866, -10903.649,   41.020,    4.538
#>       Crit. KGE2[Q]      = 0.9815
#> 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 54032...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (243 runs)
#>       Param =   11.250,  432.681,   -0.020,   83.096,    2.384
#>       Crit. KGE2[Q]      = 0.9718
#> Steepest-descent local search in progress
#>   Calibration completed (55 iterations, 775 runs)
#>       Param =    1.296,  673.308,   -0.826,  107.608,    1.081
#>       Crit. KGE2[Q]      = 0.9864
#> 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,   42.098,    2.384
#>       Crit. KGE2[Q]      = 0.9723
#> Steepest-descent local search in progress
#>   Calibration completed (171 iterations, 2023 runs)
#>       Param =    0.730,  173.202,   -0.083,    0.031,    4.650
#>       Crit. KGE2[Q]      = 0.9795
ParamV03 <- sapply(griwrm$id, function(x) {OC_OL[[x]]$Param})
save(griwrmV03, ParamV03, file = file.path(tempdir(), "V03.RData"))

Run model with this newly calibrated parameters

OM_OL <- RunModel(
  IM_OL,
  RunOptions = RunOptions,
  Param = ParamV03
)
#> RunModel.GRiwrmInputsModel: Treating sub-basin 54001...
#> RunModel.GRiwrmInputsModel: Treating sub-basin 54029...
#> RunModel.GRiwrmInputsModel: Treating sub-basin 54032...
#> RunModel.GRiwrmInputsModel: Treating sub-basin 54057...

Plot the result

As, we can see below, compare to results of vignette “V02_Calibration_SD_model,” the use of measured flow on upstream influenced basins improves largely the model performance at downstream stations.

htmltools::tagList(lapply(
  names(OM_OL),
  function(x) {
    plot(OM_OL[[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(OM_OL, "Qm3s")
plot(Qm3s[1:150,])

References

Higgs, Gary, and Geoff Petts. 1988. “Hydrological Changes and River Regulation in the UK.” Regulated Rivers: Research & Management 2 (3): 349–68. https://doi.org/10.1002/rrr.3450020312.