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

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.asp = 0.68,
  out.width = "70%",
  fig.align = "center"
)

library(airGRiwrm)
data(Severn)

The study case

In this vignette, we are still using the study case of the vignette #1 and #2, and we are going to simulate a diversion at the node "54001" to provide flow to the node "54029" in order to support low flows in the Severn river.

This objective is to maintain a minimum flow equals to the 3th decile of the flow at station "54029" without remaining less than the first decile of the flow downstream the station "54001".

Let's compute the flow quantiles at each station in m3/s:

Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
Qobs <- Qobs[, Severn$BasinsInfo$gauge_id]
Qobs_m3s <- t(apply(Qobs, 1, function(r) r * Severn$BasinsInfo$area * 1E3 / 86400))
apply(Qobs_m3s[, c("54029", "54001")], 2, quantile, probs = seq(0,1,0.1), na.rm = TRUE)

The rule to apply is expressed as follow:

The diversion from "54001" is computed to maintain a minimum flow of 5.3 m3/s at the gauging station "54029". The diversion is allowed as far as the flow at the gauging station "54001" is above 11.5 m3/s.

Network

nodes_div <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
nodes_div$model <- "RunModel_GR4J"
nodes_div <- rbind(nodes_div, data.frame(gauge_id = "54001",
                                         downstream_id = "54029",
                                         distance_downstream = 25,
                                         model = "Diversion",
                                         area = NA))
renameCols <- list(id = "gauge_id", down = "downstream_id", length = "distance_downstream")
griwrmV06 <- CreateGRiwrm(nodes_div, renameCols)
plot(griwrmV06)

GRiwrmInputsModel object

We produce below the same operations as in the vignette "V02_Calibration_SD_model" to prepare the input data:

data(Severn)
nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
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}))
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)

The main difference here is that we need to provide a diverted flow and a minimum remaining flow at station "54029". As, the diverted flow will be calculated during simulation, we provide a initial diverted flow equals to zero for all the time steps.

Qdiv <- matrix(rep(0, length(DatesR)), ncol = 1)
colnames(Qdiv) <- "54001"
Qmin <- matrix(rep(11.5 * 86400, length(DatesR)), ncol = 1)
colnames(Qmin) <- "54001"
IM_div <- CreateInputsModel(griwrmV06, DatesR, Precip, PotEvap, Qinf = Qdiv, Qmin = Qmin)

Implementation of the regulation controller

The supervisor

The simulation is piloted through a Supervisor that can contain one or more Controller. The supervision time step will be here the same as the simulation time step: 1 day. Each day, the decision is taken for the current day from the measurement simulated the previous time step.

sv <- CreateSupervisor(IM_div, TimeStep = 1L)

The control logic function

On a Diversion node, the minimum remaining flow is provided with the inputs and is automatically taken into account by the model. So the control logic function has only to take care about how much water to divert to the gauging station "54002".

We need to enclose the logic function in a function factory providing the Supervisor in the environment of the function:

#' @param sv the Supervisor environment
logicFunFactory <- function(sv) {
  #' @param Y Flow measured at "54002" the previous time step
  function(Y) {
    Qnat <- Y
    #  We need to remove the diverted flow to compute the natural flow at "54002"
    lastU <- sv$controllers[[sv$controller.id]]$U
    if (length(lastU) > 0) {
      Qnat <- max(0, Y + lastU)
    }
    return(-max(5.3 * 86400 - Qnat, 0))
  }
}

The controller

We declare the controller which defines where is the measurement Y , where to apply the decision U with which logic function:

CreateController(sv,
                 ctrl.id = "Low flow support",
                 Y = "54029",
                 U = "54001",
                 FUN = logicFunFactory(sv))

Running the simulation

First we need to create a GRiwrmRunOptions object and load the parameters calibrated in the vignette "V02_Calibration_SD_model":

# Running simulation on year 2003
IndPeriod_Run <- which(
  DatesR >= as.POSIXct("2003-03-01", tz = "UTC") &
    DatesR <= as.POSIXct("2004-01-01", tz = "UTC")
)
IndPeriod_WarmUp = seq(IndPeriod_Run[1] - 366,IndPeriod_Run[1] - 1)
RunOptions <- CreateRunOptions(IM_div,
                               IndPeriod_WarmUp = IndPeriod_WarmUp,
                               IndPeriod_Run = IndPeriod_Run)
ParamV02 <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))

The node "54029" was initially an upstream node. As it receives water from "54001" it is no longer an upstream node and needs a parameter for routing its upstream flows. We arbitrary say that the velocity in the channel between "54001" and "54029" is 1 m/s.

ParamV02$`54029` <- c(1, ParamV02$`54029`)

And we run the supervised model:

OM_div <- RunModel(sv, RunOptions = RunOptions, Param = ParamV02)

To compare results with diversion and without diversion, we run the model without the supervision (remember that we have set the diverted flow at zero in the inputs):

OM_nat <- RunModel(IM_div, RunOptions = RunOptions, Param = ParamV02)

Exploring results

Let's plot the diverted flow, and compare the flow at stations 54029 and 54001, with and without low-flow support at station 54001:

dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR,
                     Diverted_flow = OM_div$`54001`$Qdiv_m3 / 86400)

oldpar <- par(mfrow=c(3,1), mar = c(2.5,4,1,1))
plot.Qm3s(dfQdiv)

# Plot natural and influenced flow at station "54001" and "54029"
thresholds <- c("54001" = 11.5, "54029" = 5.3)
lapply(names(thresholds), function(id) {
  df <- cbind(attr(OM_div, "Qm3s")[, c("DatesR", id)],
                   attr(OM_nat, "Qm3s")[, id])
  names(df) <- c("DatesR",
                      paste(id, "with low-flow support"),
                      paste(id, "natural flow"))
  plot.Qm3s(df, ylim = c(0,50), lty = c("solid", "dashed"))
  abline(h = thresholds[id], col = "green", lty = "dotted")
})
par(oldpar)