V05 - Hydrologic model calibration

library(tidyverse)
library(airGRccia)

Network definition: GRiwrm object

Load network from V01:

nodes <- openxlsx2::read_xlsx(getCachePath("V01", "network", "xlsx"))
# Remove unecessary columns
nodes <- nodes |>
  select(id, down, length, tot_area, libelle_site, outlet_x, outlet_y)

Adding withdrawals as Direct Injection nodes:

nodes <- rbind(
  nodes,
  data.frame(
    id = paste0(nodes$id, "_WD"),
    down = nodes$id,
    length = 0,
    tot_area = NA_real_,
    libelle_site = paste("Withdrawals on", nodes$libelle_site),
    outlet_x = nodes$outlet_x,
    outlet_y = nodes$outlet_y
  )
)
nodes
#>             id     down   length   tot_area
#> 2     Y2230010 Y2300020 20.61114   78.76875
#> 3     Y2210010 Y2300020 32.87108  182.41750
#> 4     Y2100020 Y2140020 29.98713  910.30438
#> 5     Y2140020 Y2140010 17.98615 1226.60000
#> 6     Y2140010 Y2300020 13.69147 1407.18875
#> 7     Y2300020     <NA>       NA 1954.85188
#> 1  Y2230010_WD Y2230010  0.00000         NA
#> 21 Y2210010_WD Y2210010  0.00000         NA
#> 31 Y2100020_WD Y2100020  0.00000         NA
#> 41 Y2140020_WD Y2140020  0.00000         NA
#> 51 Y2140010_WD Y2140010  0.00000         NA
#> 61 Y2300020_WD Y2300020  0.00000         NA
#>                                                              libelle_site
#> 2                            Le Salagou à Clermont-l'Hérault [Mas Audran]
#> 3                                                      La Lergue à Lodève
#> 4                                                     L'Hérault à Laroque
#> 5                 L'Hérault à Saint-Guilhem-le-Désert [La Combe du Corps]
#> 6                                                      L'Hérault à Gignac
#> 7                                                     L'Hérault à Aspiran
#> 1             Withdrawals on Le Salagou à Clermont-l'Hérault [Mas Audran]
#> 21                                      Withdrawals on La Lergue à Lodève
#> 31                                     Withdrawals on L'Hérault à Laroque
#> 41 Withdrawals on L'Hérault à Saint-Guilhem-le-Désert [La Combe du Corps]
#> 51                                      Withdrawals on L'Hérault à Gignac
#> 61                                     Withdrawals on L'Hérault à Aspiran
#>    outlet_x outlet_y
#> 2  733408.0  6285655
#> 3  726110.0  6292216
#> 4  759105.0  6313209
#> 5  748274.3  6294640
#> 6  743142.0  6283755
#> 7  737989.0  6274704
#> 1  733408.0  6285655
#> 21 726110.0  6292216
#> 31 759105.0  6313209
#> 41 748274.3  6294640
#> 51 743142.0  6283755
#> 61 737989.0  6274704

Load the observations:

BasinsObs <- readRDS(getCachePath("V02", "BasinsObs_safran_historical", "RDS"))
str(BasinsObs)
#> List of 5
#>  $ DatesR  : POSIXlt[1:4536], format: "2012-08-01" "2012-08-02" ...
#>  $ Precip  :'data.frame':    4536 obs. of  6 variables:
#>   ..$ Y2230010: num [1:4536] 0 0 0 0.324 3.24 ...
#>   ..$ Y2210010: num [1:4536] 0 0 0 1.79 8.5 ...
#>   ..$ Y2100020: num [1:4536] 0 0 0 6.67 13.53 ...
#>   ..$ Y2140020: num [1:4536] 0 0 0 0.945 3.338 ...
#>   ..$ Y2140010: num [1:4536] 0 0 0 0.474 3.602 ...
#>   ..$ Y2300020: num [1:4536] 0 0 0 0.467 3.78 ...
#>  $ PotEvap :'data.frame':    4536 obs. of  6 variables:
#>   ..$ Y2230010: num [1:4536] 4.86 9 6.84 6.77 5.41 ...
#>   ..$ Y2210010: num [1:4536] 5.76 8.43 7.04 6.76 5.71 ...
#>   ..$ Y2100020: num [1:4536] 6.44 7.43 7.21 5.83 4.34 ...
#>   ..$ Y2140020: num [1:4536] 6.24 7.1 6.27 6.21 5.29 ...
#>   ..$ Y2140010: num [1:4536] 5.22 8.57 6.9 6.56 5.42 ...
#>   ..$ Y2300020: num [1:4536] 4.99 8.84 7.18 6.66 5.43 ...
#>  $ TempMean:'data.frame':    4536 obs. of  6 variables:
#>   ..$ Y2230010: num [1:4536] 23.3 26 25.1 25.1 23 ...
#>   ..$ Y2210010: num [1:4536] 21.5 23.1 22.2 22.2 20.3 ...
#>   ..$ Y2100020: num [1:4536] 22.4 23.4 22.6 22.4 19.9 ...
#>   ..$ Y2140020: num [1:4536] 23.1 25 25.5 24.8 23.1 ...
#>   ..$ Y2140010: num [1:4536] 23.3 25.9 25.5 25.1 23.2 ...
#>   ..$ Y2300020: num [1:4536] 23.4 26.3 25.6 25.3 23.3 ...
#>  $ Qobs    : tibble [4,536 × 5] (S3: tbl_df/tbl/data.frame)
#>   ..$ Y2210010: num [1:4536] 0.769 0.666 0.661 0.71 0.702 0.726 0.722 0.697 0.7 0.669 ...
#>   ..$ Y2100020: num [1:4536] 2.91 2.9 2.88 2.86 2.92 ...
#>   ..$ Y2140020: num [1:4536] 1.61 1.55 1.55 1.5 1.52 ...
#>   ..$ Y2140010: num [1:4536] 2.24 2.42 2.51 2.46 2.44 ...
#>   ..$ Y2300020: num [1:4536] 3.67 3.67 3.56 3.39 3.4 ...
#>   ..- attr(*, "codes_entite_station")='data.frame':  8 obs. of  2 variables:
#>   .. ..$ code_entite : chr [1:8] "Y2230010" "Y2210010" "Y2210010" "Y2100020" ...
#>   .. ..$ code_station: chr [1:8] "Y223001001" "Y221001001" "Y221001002" "Y210002001" ...
#>   ..- attr(*, "DatesR")= POSIXlt[1:4536], format: "2012-08-01" "2012-08-02" ...
#>  - attr(*, "class")= chr [1:2] "BasinsObs" "list"

No observed flow is available at station Y2230010 which corresponds to location of the Salagou dam.

So we need to:

  • Add a reservoir in the model corresponding to the Salagou dam
  • model the reservoir inflows using an ungauged model with parameters transferred from a nearby basin.
# Add the dam
nodes <- rbind(
  nodes,
  data.frame(
    id = "DAM",
    down = "Y2300020",
    length = 48.50,
    tot_area = as.numeric(NA),
    libelle_site = "Barrage du Salagou",
    outlet_x = nodes$outlet_x[nodes$id == "Y2230010"],
    outlet_y = nodes$outlet_y[nodes$id == "Y2230010"]
  )
)
# Adjust Dam inflow basin
nodes$down[nodes$id == "Y2230010"] <- "DAM"
nodes$length[nodes$id == "Y2230010"] <- 0

Choice of the models used:

nodes$model <- "RunModel_GR4J"
nodes$model[nodes$id == "Y2230010"] <- "Ungauged"
nodes$model[nodes$id == "DAM"] <- "RunModel_Reservoir"
nodes$model[grep("_WD", nodes$id)] <- NA_character_

Definition of the basin given its parameters for the ungauged basin Y2230010. The basin “La Lergue à Lodève (Y2210010)” is the closest one.

nodes$donor <- as.character(NA)
nodes$donor[nodes$id == "Y2230010"] <- "Y2210010"

Creation of the GRiwrm object:

griwrm <- CreateGRiwrm(nodes, cols = list(area = "tot_area"))
plot(griwrm)

Inputs of the model: createInputsModel

Format all the inputs:

Qinf <- openxlsx2::read_xlsx(getCachePath(
  "V03",
  "Qinf_historical",
  ".xlsx"
)) %>%
  mutate(DatesR = as.Date(DatesR))
Qinf <- Qinf %>%
  filter(DatesR >= min(BasinsObs$DatesR) & DatesR <= max(BasinsObs$DatesR))
names(Qinf)[-1] <- paste0(names(Qinf)[-1], "_WD")
str(Qinf)
#> 'data.frame':    4536 obs. of  7 variables:
#>  $ DatesR     : Date, format: "2012-08-01" "2012-08-02" ...
#>  $ Y2100020_WD: num  5258 5258 5258 5258 5258 ...
#>  $ Y2140010_WD: num  95933 95933 95933 95933 95933 ...
#>  $ Y2140020_WD: num  329 329 329 329 329 ...
#>  $ Y2210010_WD: num  810 810 810 810 810 ...
#>  $ Y2230010_WD: num  1289 1289 1289 1289 1289 ...
#>  $ Y2300020_WD: num  12949 12949 12949 12949 12949 ...
Qrelease <- openxlsx2::read_xlsx(getCachePath(
  "V04",
  "Qrelease_historical",
  "xlsx"
))
str(Qrelease)
#> 'data.frame':    4536 obs. of  2 variables:
#>  $ DatesR: POSIXct, format: "2012-08-01" "2012-08-02" ...
#>  $ DAM   : num  0.34 0.34 0.34 0.34 0.34 0.34 0.34 0.34 0.34 0.34 ...

Remove extra columns and convert flows into m3/day:

Qinf$DatesR <- NULL
Qrelease$DatesR <- NULL
Qinf <- -Qinf # withdrawals are negative flows
Qrelease <- Qrelease * 86400 # m3/s -> m3/day

Call CreateInputsModel:

IM <- CreateInputsModel(
  griwrm,
  DatesR = BasinsObs$DatesR,
  Precip = BasinsObs$Precip,
  PotEvap = BasinsObs$PotEvap,
  Qinf = Qinf,
  Qrelease = Qrelease
)
#> CreateInputsModel.GRiwrm: Processing sub-basin Y2210010...
#> CreateInputsModel.GRiwrm: Processing sub-basin Y2100020...
#> CreateInputsModel.GRiwrm: Processing sub-basin Y2230010...
#> CreateInputsModel.GRiwrm: Processing sub-basin Y2140020...
#> CreateInputsModel.GRiwrm: Processing sub-basin Y2140010...
#> CreateInputsModel.GRiwrm: Processing sub-basin DAM...
#> CreateInputsModel.GRiwrm: Processing sub-basin Y2300020...

CreateRunoptions

Definition of the calibration period

Ind_Run <- seq(
  which(BasinsObs$DatesR == as.POSIXlt("2013-08-01", tz = "UTC")),
  which(BasinsObs$DatesR == as.POSIXlt("2023-07-31", tz = "UTC"))
)
RO <- CreateRunOptions(IM, IndPeriod_Run = Ind_Run)
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, : model warm up period not defined: default configuration used
#>   the year preceding the run period is used
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, : model warm up period not defined: default configuration used
#>   the year preceding the run period is used
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, : model warm up period not defined: default configuration used
#>   the year preceding the run period is used
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, : model warm up period not defined: default configuration used
#>   the year preceding the run period is used
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, : model warm up period not defined: default configuration used
#>   the year preceding the run period is used
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, : model warm up period not defined: default configuration used
#>   the year preceding the run period is used
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL,
#> IndPeriod_Run, : model states initialisation not defined: default configuration
#> used
#> Warning in (function (FUN_MOD, InputsModel, IndPeriod_WarmUp = NULL, IndPeriod_Run, : model warm up period not defined: default configuration used
#>   the year preceding the run period is used

It seems that in 2012, the dam was barely full.

RO[["DAM"]]$IniStates <- c("Reservoir.V" = 100E6)

CreateInputsCrit

Observed flows need to be converted from m3/s to mm/day.

Qobs <- BasinsObs$Qobs[Ind_Run, ]
areas <- sapply(names(Qobs), function(x) griwrm$area[griwrm$id == x])
k <- 86400 * 1E3 / (areas * 1E6) # seconds in a day * m3->L / area in m2 = L/m2/day = mm/day
Qobs <- t(replicate(nrow(Qobs), k)) * Qobs
summary(Qobs)
#>     Y2210010          Y2100020          Y2140020           Y2140010      
#>  Min.   : 0.1994   Min.   : 0.1391   Min.   : 0.07424   Min.   : 0.1040  
#>  1st Qu.: 0.4751   1st Qu.: 0.3179   1st Qu.: 0.25978   1st Qu.: 0.2257  
#>  Median : 0.8047   Median : 0.7637   Median : 0.69185   Median : 0.6563  
#>  Mean   : 1.6277   Mean   : 1.6148   Mean   : 1.48510   Mean   : 1.4231  
#>  3rd Qu.: 1.5952   3rd Qu.: 1.5964   3rd Qu.: 1.60466   3rd Qu.: 1.4193  
#>  Max.   :57.6793   Max.   :45.0272   Max.   :43.20078   Max.   :50.6621  
#>  NAs    :129                         NAs    :475                         
#>     Y2300020       
#>  Min.   : 0.09175  
#>  1st Qu.: 0.26368  
#>  Median : 0.59601  
#>  Mean   : 1.29645  
#>  3rd Qu.: 1.19281  
#>  Max.   :44.69826  
#>  NAs    :11
IC <- CreateInputsCrit(
  InputsModel = IM,
  FUN_CRIT = ErrorCrit_KGE2,
  RunOptions = RO,
  Obs = Qobs
)

CreateCalibOptions

CO <- CreateCalibOptions(
  IM,
  FixedParam = list(DAM = c(Vmax = 102.2E6, celerity = 1))
)

Calibration

OC <- Calibration(
  InputsModel = IM,
  RunOptions = RO,
  InputsCrit = IC,
  CalibOptions = CO
)
#> Calibration.GRiwrmInputsModel: Processing sub-basin 'Y2210010'...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (243 runs)
#>       Param =    1.250,  432.681,   -0.020,   83.096,    1.417
#>       Crit. KGE2[Q]      = 0.6498
#> Steepest-descent local search in progress
#>   Calibration completed (44 iterations, 664 runs)
#>       Param =    1.250,  702.786,    1.278,   44.305,    1.379
#>       Crit. KGE2[Q]      = 0.8785
#> Calibration.GRiwrmInputsModel: Processing sub-basin 'Y2100020'...
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (243 runs)
#>       Param =    1.250,  432.681,   -0.020,   83.096,    1.417
#>       Crit. KGE2[Q]      = 0.7841
#> Steepest-descent local search in progress
#>   Calibration completed (37 iterations, 605 runs)
#>       Param =    1.250,  633.937,    0.294,   29.363,    1.447
#>       Crit. KGE2[Q]      = 0.8028
#> Calibration.GRiwrmInputsModel: Processing sub-basin 'Y2230010'...
#> Tranferring parameters from node 'Y2210010' to node 'Y2230010'
#>       Param =    1.250,  702.786,    1.278,   44.305,    1.072
#> 
#> Calibration.GRiwrmInputsModel: Processing sub-basin 'Y2140020'...
#> Parameter regularization: get a priori parameters from node Y2100020: 1.25, 633.937, 0.294, 29.363, 1.054
#> Crit. KGE2[Q] = 0.8280
#>  SubCrit. KGE2[Q] cor(sim, obs, "pearson") = 0.8339 
#>  SubCrit. KGE2[Q] cv(sim)/cv(obs)          = 0.9971 
#>  SubCrit. KGE2[Q] mean(sim)/mean(obs)      = 0.9553 
#> 
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (243 runs)
#>       Param =    1.250,  169.017,   -0.020,   83.096,    1.417
#>       Crit. Composite    = 0.9060
#> Steepest-descent local search in progress
#>   Calibration completed (26 iterations, 481 runs)
#>       Param =    0.750,   66.686,    0.555,  123.965,    1.349
#>       Crit. Composite    = 0.9165
#>  Formula: sum(0.87 * KGE2[sqrt(Q)], 0.13 * GAPX[ParamT])
#> Calibration.GRiwrmInputsModel: Processing sub-basin 'Y2140010'...
#> Parameter regularization: get a priori parameters from node Y2140020: 0.75, 66.686, 0.555, 123.965, 1.14
#> Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 12
#> time steps with negative flow, set to zero.
#> Crit. KGE2[Q] = 0.8663
#>  SubCrit. KGE2[Q] cor(sim, obs, "pearson") = 0.8847 
#>  SubCrit. KGE2[Q] cv(sim)/cv(obs)          = 0.9370 
#>  SubCrit. KGE2[Q] mean(sim)/mean(obs)      = 1.0249 
#> 
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (243 runs)
#>       Param =    1.250,  169.017,   -0.020,   83.096,    1.417
#>       Crit. Composite    = 0.9145
#> Steepest-descent local search in progress
#>   Calibration completed (31 iterations, 531 runs)
#>       Param =    0.800,   47.942,    0.797,  198.343,    1.252
#>       Crit. Composite    = 0.9323
#>  Formula: sum(0.87 * KGE2[sqrt(Q)], 0.13 * GAPX[ParamT])
#> Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 4
#> time steps with negative flow, set to zero.
#> Calibration.GRiwrmInputsModel: Processing sub-basin 'DAM'...
#> Parameters already fixed - no need for calibration
#>       Param = 102200000.000,    1.000
#> Calibration.GRiwrmInputsModel: Processing sub-basin 'Y2300020'...
#> Parameter regularization: test a priori parameters from node Y2210010: 1.25, 702.786, 1.278, 44.305, 1.579
#> Crit. KGE2[Q] = 0.8373
#>  SubCrit. KGE2[Q] cor(sim, obs, "pearson") = 0.8995 
#>  SubCrit. KGE2[Q] cv(sim)/cv(obs)          = 0.8731 
#>  SubCrit. KGE2[Q] mean(sim)/mean(obs)      = 1.0166 
#> 
#> Parameter regularization: test a priori parameters from node Y2140010: 0.8, 47.942, 0.797, 198.343, 1.437
#> Crit. KGE2[Q] = 0.8492
#>  SubCrit. KGE2[Q] cor(sim, obs, "pearson") = 0.9049 
#>  SubCrit. KGE2[Q] cv(sim)/cv(obs)          = 0.9009 
#>  SubCrit. KGE2[Q] mean(sim)/mean(obs)      = 1.0622 
#> 
#> Parameter regularization: test a priori parameters from node Y2230010: 1.25, 702.786, 1.278, 44.305, 1.579
#> Crit. KGE2[Q] = 0.8373
#>  SubCrit. KGE2[Q] cor(sim, obs, "pearson") = 0.8995 
#>  SubCrit. KGE2[Q] cv(sim)/cv(obs)          = 0.8731 
#>  SubCrit. KGE2[Q] mean(sim)/mean(obs)      = 1.0166 
#> 
#> Parameter regularization: set a priori parameters from node Y2140010: 0.8, 47.942, 0.797, 198.343, 1.437
#> Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
#>   Screening completed (243 runs)
#>       Param =    1.250,  169.017,   -0.020,   83.096,    1.417
#>       Crit. Composite    = 0.9309
#> Steepest-descent local search in progress
#>   Calibration completed (24 iterations, 461 runs)
#>       Param =    0.860,   49.899,    0.721,  179.469,    1.417
#>       Crit. Composite    = 0.9406
#>  Formula: sum(0.87 * KGE2[sqrt(Q)], 0.13 * GAPX[ParamT])
params <- extractParam(OC)
str(params)
#> List of 7
#>  $ Y2210010: num [1:5] 1.25 702.79 1.28 44.3 1.38
#>  $ Y2100020: num [1:5] 1.25 633.937 0.294 29.363 1.447
#>  $ Y2230010: num [1:5] 1.25 702.79 1.28 44.3 1.07
#>  $ Y2140020: num [1:5] 0.75 66.686 0.555 123.965 1.349
#>  $ Y2140010: num [1:5] 0.8 47.942 0.797 198.343 1.252
#>  $ DAM     : Named num [1:2] 1.02e+08 1.00
#>   ..- attr(*, "names")= chr [1:2] "Vmax" "celerity"
#>  $ Y2300020: num [1:5] 0.86 49.899 0.721 179.469 1.417

Check calibration result

OM <- RunModel(
  IM,
  RunOptions = RO,
  Param = params
)
#> RunModel.GRiwrmInputsModel: Processing sub-basin Y2210010...
#> RunModel.GRiwrmInputsModel: Processing sub-basin Y2100020...
#> RunModel.GRiwrmInputsModel: Processing sub-basin Y2230010...
#> RunModel.GRiwrmInputsModel: Processing sub-basin Y2140020...
#> RunModel.GRiwrmInputsModel: Processing sub-basin Y2140010...
#> Warning in RunModel_Lag(InputsModel, RunOptions, Param[1], OutputsModel): 4
#> time steps with negative flow, set to zero.
#> RunModel.GRiwrmInputsModel: Processing sub-basin DAM...
#> RunModel.GRiwrmInputsModel: Processing sub-basin Y2300020...
plot(OM, Qobs = Qobs[Ind_Run, ])

#> Warning in FUN(X[[i]], ...): Column "Y2230010" not found in Qobs
#> Warning in plot.OutputsModel(x[[id]], Qobs = Qobs_id, Vobs = Vobs_id, BasinArea
#> = BasinArea, : the 'Error' and 'CorQQ' plot(s) cannot be drawn if there is no
#> 'Qobs'

#> Warning in plot.OutputsModel(x[[id]], Qobs = Qobs_id, Vobs = Vobs_id, BasinArea
#> = BasinArea, : zeroes detected in 'Qsim': some plots in the log space will not
#> be created using all time-steps

#> Warning in FUN(X[[i]], ...): Column "DAM" not found in Qobs