diff --git a/DESCRIPTION b/DESCRIPTION
index dbb0c5f3b2558bb9c2620634cbc3bb8d8fc4c062..c10112e14b3b68b011a6807d2b04d16b4e8dfab5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -25,6 +25,7 @@ Imports:
     dplyr,
     graphics,
     grDevices,
+    rlang,
     utils
 Suggests:
     htmltools,
diff --git a/NAMESPACE b/NAMESPACE
index fdbec6d4c899b152ef69bf5ecda1410d4f5a49cf..a53032b9529e6f6ef859b576060be08a7564b061 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -21,7 +21,6 @@ S3method(CreateRunOptions,character)
 S3method(RunModel,GR)
 S3method(RunModel,GRiwrmInputsModel)
 S3method(RunModel,InputsModel)
-S3method(RunModel,SD)
 S3method(RunModel,Supervisor)
 S3method(isNodeDownstream,GRiwrm)
 S3method(isNodeDownstream,GRiwrmInputsModel)
@@ -51,11 +50,11 @@ export(plot.Qm3s)
 import(airGR)
 import(dplyr)
 importFrom(dplyr,"%>%")
-importFrom(grDevices,rainbow)
 importFrom(graphics,matplot)
 importFrom(graphics,par)
 importFrom(graphics,plot)
 importFrom(graphics,title)
+importFrom(rlang,.data)
 importFrom(stats,setNames)
 importFrom(utils,read.table)
 importFrom(utils,tail)
diff --git a/R/Calibration.GRiwrmInputsModel.R b/R/Calibration.GRiwrmInputsModel.R
index af7080f0cf1fe6d3cefede01b2743e173a137d5c..9d65235cb4eccb33aaddca27864fb9db39216bf8 100644
--- a/R/Calibration.GRiwrmInputsModel.R
+++ b/R/Calibration.GRiwrmInputsModel.R
@@ -34,7 +34,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
   b <- sapply(InputsModel, function(IM) !IM$isUngauged)
   gaugedIds <- names(b[b])
 
-  for(id in gaugedIds) {
+  for (id in gaugedIds) {
     IM <- InputsModel[[id]]
     message("Calibration.GRiwrmInputsModel: Processing sub-basin ", id, "...")
 
@@ -79,18 +79,18 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
       # Select nodes with model in the sub-network
       g <- attr(IM, "GRiwrm")
       Ids <- g$id[!is.na(g$donor) & g$donor == id]
-      if(IM[[id]]$model$hasX4) {
+      if (IM[[id]]$model$hasX4) {
         # Extract the X4 calibrated for the whole intermediate basin
         X4 <- OutputsCalib[[id]]$ParamFinalR[IM[[id]]$model$iX4] # Global parameter
       }
       for (uId in Ids) {
-        if(!IM[[uId]]$isReservoir) {
+        if (!IM[[uId]]$isReservoir) {
           # Add OutputsCalib for ungauged nodes
           OutputsCalib[[uId]] <- OutputsCalib[[id]]
           # Copy parameters and transform X4 relatively to the sub-basin area
           OutputsCalib[[uId]]$ParamFinalR <-
             OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$indexParamUngauged]
-          if(IM[[id]]$model$hasX4) {
+          if (IM[[id]]$model$hasX4) {
             subBasinAreas <- calcSubBasinAreas(IM)
             OutputsCalib[[uId]]$ParamFinalR[IM[[uId]]$model$iX4] <- max(
               X4 * (subBasinAreas[uId] / subBasinAreas[id]) ^ 0.3,
@@ -107,7 +107,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
           )
         }
       }
-      if(useUpstreamQsim) {
+      if (useUpstreamQsim) {
         OM_subnet <- RunModel_Ungauged(IM,
                                        RunOptions[[id]],
                                        OutputsCalib[[id]]$ParamFinalR,
@@ -115,7 +115,7 @@ Calibration.GRiwrmInputsModel <- function(InputsModel,
         OutputsModel <- c(OutputsModel, OM_subnet)
       }
       IM <- IM[[id]]
-    } else if(useUpstreamQsim) {
+    } else if (useUpstreamQsim) {
       # Run the model for the sub-basin
       OutputsModel[[id]] <- RunModel(
         x = IM,
@@ -150,7 +150,7 @@ getInputsCrit_Lavenne <- function(id, OutputsModel, InputsCrit) {
   AprCelerity <- attr(InputsCrit[[id]], "AprCelerity")
   Lavenne_FUN <- attr(InputsCrit[[id]], "Lavenne_FUN")
   AprParamR <- OutputsModel[[AprioriId]]$RunOptions$Param
-  if(!inherits(OutputsModel[[AprioriId]], "SD")) {
+  if (!inherits(OutputsModel[[AprioriId]], "SD")) {
     # Add default velocity parameter for a priori upstream catchment
     AprParamR <- c(AprCelerity, AprParamR)
   }
@@ -199,6 +199,7 @@ reduceGRiwrmObj4Ungauged <- function(griwrm, obj) {
 #' - `RunOptions`: a *GRiwrmRunOptions* of the reduced network
 #' @noRd
 #' @importFrom dplyr "%>%"
+#' @importFrom rlang .data
 #'
 updateParameters4Ungauged <- function(GaugedId,
                                       InputsModel,
@@ -211,12 +212,12 @@ updateParameters4Ungauged <- function(GaugedId,
   # Select nodes identified with the current node as donor gauged node
   griwrm <- attr(InputsModel, "GRiwrm")
   donorIds <- griwrm$id[!is.na(griwrm$donor) & griwrm$donor == GaugedId]
-  gDonor <- griwrm %>% dplyr::filter(id %in% donorIds)
+  gDonor <- griwrm %>% dplyr::filter(.data$id %in% donorIds)
   # Add upstream nodes for routing upstream flows
   upNodes <- griwrm %>%
-    dplyr::filter(down %in% gDonor$id,
-                  !id %in% gDonor$id) %>%
-    dplyr::mutate(model = ifelse(!is.na(model) & model != "Diversion", NA, model))
+    dplyr::filter(.data$down %in% gDonor$id,
+                  !.data$id %in% gDonor$id) %>%
+    dplyr::mutate(model = ifelse(!is.na(.data$model) & .data$model != "Diversion", NA, .data$model))
   upIds <- upNodes$id
   g <- rbind(upNodes, gDonor)
   # Set downstream nodes
@@ -236,7 +237,7 @@ updateParameters4Ungauged <- function(GaugedId,
   # Update Qupstream already modeled in the reduced network upstream nodes
   idIM <- unique(g$down[g$id %in% upIds])
   for (id in idIM) {
-    if(useUpstreamQsim && any(InputsModel[[id]]$UpstreamIsModeled)) {
+    if (useUpstreamQsim && any(InputsModel[[id]]$UpstreamIsModeled)) {
       # Temporarily switch off upstream nodes belonging to the donor basin
       UpIsModeledBackUp <- InputsModel[[id]]$UpstreamIsModeled
       ImUpIds <- InputsModel[[id]]$UpstreamNodes
@@ -269,7 +270,7 @@ updateParameters4Ungauged <- function(GaugedId,
 calcSubBasinAreas <- function(IM) {
   unlist(
     sapply(IM, function(x) {
-      if(is.list(x)) as.numeric(x$BasinAreas[length(x$BasinAreas)])})
+      if (is.list(x)) as.numeric(x$BasinAreas[length(x$BasinAreas)])})
   )
 }
 
@@ -293,7 +294,7 @@ calcSubBasinAreas <- function(IM) {
 #' @references Lobligeois, Florent. Mieux connaître la distribution spatiale des
 #' pluies améliore-t-il la modélisation des crues ? Diagnostic sur 181 bassins
 #' versants français. Phdthesis, AgroParisTech, 2014.
-#' <https://pastel.archives-ouvertes.fr/tel-01134990/document>
+#' <https://pastel.hal.science/tel-01134990/document>
 #'
 #' @inheritParams airGR::RunModel
 #' @param ouput.all [logical] if `TRUE` returns the output of [RunModel.GRiwrm],
@@ -311,7 +312,7 @@ RunModel_Ungauged <- function(InputsModel, RunOptions, Param, output.all = FALSE
       return(IM$FixedParam)
     }
     p <- Param[IM$model$indexParamUngauged]
-    if(IM$model$hasX4) {
+    if (IM$model$hasX4) {
       p[IM$model$iX4] <- max(
         Param[InputsModel[[donor]]$model$iX4] *
           (IM$BasinAreas[length(IM$BasinAreas)] / donorArea) ^ 0.3,
diff --git a/R/ConvertMeteoSD.R b/R/ConvertMeteoSD.R
index 9dc8a6caafcbbbcb8a4212f3f68636bdfa6eb82c..e37351214aa586183d1953a5ad4515a025b09f51 100644
--- a/R/ConvertMeteoSD.R
+++ b/R/ConvertMeteoSD.R
@@ -34,7 +34,7 @@ ConvertMeteoSD.GRiwrm <- function(x, meteo, ...) {
 ConvertMeteoSD.character <- function(x, griwrm, meteo, ...) {
   griwrm <- griwrm[getDiversionRows(griwrm, inverse = TRUE), ]
   upperIDs <- getUpstreamRunOffIds(x, griwrm)
-  if(length(upperIDs) == 1) {
+  if (length(upperIDs) == 1) {
     return(meteo[,x])
   }
   areas <- griwrm$area[match(upperIDs, griwrm$id)]
@@ -52,29 +52,29 @@ ConvertMeteoSD.character <- function(x, griwrm, meteo, ...) {
 #' @rdname ConvertMeteoSD
 ConvertMeteoSD.matrix <- function(x, areas, temperature = FALSE, ...) {
   # Check arguments
-  if(nrow(x) < 2) {
+  if (nrow(x) < 2) {
     stop("Meteorological data matrix should contain more than one row")
   }
-  if(length(areas) != ncol(x)) {
+  if (length(areas) != ncol(x)) {
     stop("'areas' length and meteo data matrix number of columns should be equal")
   }
-  if(areas[1] <= sum(areas[-1])) {
+  if (areas[1] <= sum(areas[-1])) {
     stop("Basin area 'areas[1]' should be greater than the sum of the upstream sub-basin areas")
   }
-  if(ncol(x) == 1) {
+  if (ncol(x) == 1) {
     return(x)
   }
   # Convert mm to 1E3 m3
   V <- x * rep(areas, rep(nrow(x), length(areas)))
   # Sum upstream data
-  if(ncol(x) > 2) {
+  if (ncol(x) > 2) {
     Vup <- rowSums(V[,-1])
   } else {
     Vup <- V[,2]
   }
   # Remove to basin to get downstream data
   Vdown <- V[,1] - Vup
-  if(!temperature) Vdown[Vdown < 0] <- 0
+  if (!temperature) Vdown[Vdown < 0] <- 0
   # Convert to mm
   meteoDown <- Vdown / (areas[1] - sum(areas[-1]))
   return(as.matrix(meteoDown, ncol = 1))
@@ -85,7 +85,7 @@ getUpstreamRunOffIds <- function(id, griwrm) {
   upstreamNodeIds <- griwrm$id[griwrm$down == id & !is.na(griwrm$down)]
   upstreamRunOffIds <- griwrm$id[griwrm$id %in% upstreamNodeIds & !is.na(griwrm$area)]
   upstreamNaAreaIds <- upstreamNodeIds[!upstreamNodeIds %in% upstreamRunOffIds]
-  if(length(upstreamNaAreaIds) > 0) {
+  if (length(upstreamNaAreaIds) > 0) {
     upstreamRunOffIds <-  c(
       upstreamRunOffIds,
       unlist(sapply(upstreamNaAreaIds, getUpstreamRunOffIds, griwrm = griwrm))
diff --git a/R/CreateCalibOptions.GRiwrmInputsModel.R b/R/CreateCalibOptions.GRiwrmInputsModel.R
index 8d06d706da31f21c159308696824e74af281292a..7367dfe646b3e54644eed1cccaf01fffc6542dc8 100644
--- a/R/CreateCalibOptions.GRiwrmInputsModel.R
+++ b/R/CreateCalibOptions.GRiwrmInputsModel.R
@@ -43,7 +43,7 @@ CreateCalibOptions.GRiwrmInputsModel <- function(x, FixedParam = NULL, ...) {
   CalibOptions <- list()
   class(CalibOptions) <- c("GRiwrmCalibOptions", class(CalibOptions))
 
-  for(id in gaugedIds) {
+  for (id in gaugedIds) {
     IM <- x[[id]]
     FP <- NULL
     if (!is.null(FixedParam)) {
diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index 42352c618bc36abd19115655d367b5252a5340b3..2e9cb6d6d0b5a6e7fb43baef5047c16587dfe05a 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -42,7 +42,7 @@ CreateCalibOptions.InputsModel <- function(x, FixedParam = NULL, ...) {
   dots <- list(...)
   # Add FUN_MOD in parameters if carried by InputsModel
   if (!"FUN_MOD" %in% names(dots)) {
-    if(!is.null(x$FUN_MOD)) {
+    if (!is.null(x$FUN_MOD)) {
       dots$FUN_MOD <- x$FUN_MOD
     } else {
       stop(" The parameter `FUN_MOD` must be defined")
diff --git a/R/CreateController.R b/R/CreateController.R
index 6cad87fe8e7706a1fb528f15972c7936cc445cf3..fe1fb7dace016e5177dc0f7729ccf64165bacf28 100644
--- a/R/CreateController.R
+++ b/R/CreateController.R
@@ -47,10 +47,10 @@ CreateController <- function(supervisor, ctrl.id, Y, U, FUN){
 
   # Function called from Supervisor environment
   #environment(ctrlr$FUN) <- supervisor
-  if(!is.null(supervisor$controllers[[ctrl.id]])) {
+  if (!is.null(supervisor$controllers[[ctrl.id]])) {
     warning("The existing controller '", ctrl.id, "' has been overwritten in the supervisor")
   } else {
-    message("The controller has been added to the supervisor")
+    message("The controller '", ctrl.id, "' has been added to the supervisor")
   }
   supervisor$controllers[[ctrl.id]] <- ctrlr
   invisible(ctrlr)
diff --git a/R/CreateGRiwrm.R b/R/CreateGRiwrm.R
index 6ff7945830ecb56f856eec834b2f6a56de19b9a3..0ceb06b6e34bd2c8ae9413d1f03470ed40dc7dcd 100644
--- a/R/CreateGRiwrm.R
+++ b/R/CreateGRiwrm.R
@@ -150,7 +150,7 @@ getNodeRanking <- function(griwrm) {
       upId <- upIds[1]
       #Browse the ungauged sub-network until the donor
       upDonor <- g$donor[g$id == upId]
-      g2 <- g %>% filter(donor == upDonor)
+      g2 <- g[g$donor == upDonor, ]
       g2$donor <- g2$id
       ungaugedIds <- getNodeRanking(g2)
       upIds <- upIds[!upIds %in% ungaugedIds]
@@ -187,7 +187,7 @@ checkNetworkConsistency <- function(db) {
   })
   id_reservoirs <- db3$id[db3$model == "RunModel_Reservoir"]
   sapply(id_reservoirs, function(id) {
-    if(length(db$id[!is.na(db$down) & db$down == id]) == 0) {
+    if (length(db$id[!is.na(db$down) & db$down == id]) == 0) {
       stop("The reservoir ", id,
            " must have at least one upstream node as inflows.")
     }
@@ -246,7 +246,7 @@ getGaugedId <- function(id, griwrm) {
     id_down <- g2$down[g2$id == id]
     if (!is.na(id_down)) {
       return(getGaugedId(id_down, griwrm))
-    } else if(length(getDiversionRows(griwrm)) > 0) {
+    } else if (length(getDiversionRows(griwrm)) > 0) {
       # Search on Diversion
       g3 <- griwrm[getDiversionRows(griwrm), ]
       id_down <- g3$down[g3$id == id]
@@ -263,7 +263,7 @@ getDiversionRows <- function(griwrm, inverse = FALSE) {
 
   rows <- which(!is.na(griwrm$model) & griwrm$model == "Diversion")
   if (inverse) {
-    if(length(rows) == 0) {
+    if (length(rows) == 0) {
       rows <- seq.int(nrow(griwrm))
     } else {
       rows <- setdiff(seq.int(nrow(griwrm)), rows)
@@ -276,7 +276,7 @@ setDonor <- function(griwrm) {
   sapply(seq(nrow(griwrm)), function(i) {
     id <- griwrm$id[i]
     model <- griwrm$model[i]
-    if(model == "RunModel_Reservoir" && is.na(griwrm$down[i])){
+    if (model == "RunModel_Reservoir" && is.na(griwrm$down[i])){
       # RunModel_Reservoir needs to be its own "donor" only if at downstream
       # Otherwise we search the first gauged station downstream to allow
       # calibration with ungauged upstream nodes
diff --git a/R/CreateInputsCrit.GRiwrmInputsModel.R b/R/CreateInputsCrit.GRiwrmInputsModel.R
index f78a11c5105eb34c5ba3e4a3f12485386f3163eb..e161a15dd6ac25dca01e9bdd56e10cc50f8db355 100644
--- a/R/CreateInputsCrit.GRiwrmInputsModel.R
+++ b/R/CreateInputsCrit.GRiwrmInputsModel.R
@@ -74,7 +74,7 @@ CreateInputsCrit.GRiwrmInputsModel <- function(InputsModel,
 
   np <- getAllNodesProperties(attr(InputsModel, "GRiwrm"))
   gaugedIds <- np$id[np$calibration == "Gauged"]
-  for(id in gaugedIds) {
+  for (id in gaugedIds) {
     if (id %in% colnames(Obs)) {
       IM <- InputsModel[[id]]
       InputsCrit[[IM$id]] <- CreateInputsCrit.InputsModel(
diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 23306b4b47de50993514511b8c5283d60ab578cb..b2dbeec36ae5fcab30f9163481a0c73d07dea867 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -175,7 +175,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
 
   InputsModel <- CreateEmptyGRiwrmInputsModel(x)
 
-  for(id in getNodeRanking(x)) {
+  for (id in getNodeRanking(x)) {
     message("CreateInputsModel.GRiwrm: Processing sub-basin ", id, "...")
 
     InputsModel[[id]] <-
@@ -248,7 +248,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, ..., Qobs, Qmin, Qrel
   LengthHydro <- NULL
   BasinAreas <- NULL
 
-  if(length(UpstreamNodeRows) > 0) {
+  if (length(UpstreamNodeRows) > 0) {
     # Sub-basin with hydraulic routing
     Qupstream <- NULL
     Qupstream <- as.matrix(cbind(
@@ -309,7 +309,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, ..., Qobs, Qmin, Qrel
   # Add Identifiers of connected nodes in order to be able to update them with simulated flows
   InputsModel$id <- id
   InputsModel$down <- node$down
-  if(length(UpstreamNodeRows) > 0) {
+  if (length(UpstreamNodeRows) > 0) {
     InputsModel$UpstreamNodes <- griwrm$id[UpstreamNodeRows]
     InputsModel$UpstreamIsModeled <- !is.na(griwrm$model[UpstreamNodeRows])
     names(InputsModel$UpstreamIsModeled) <- InputsModel$UpstreamNodes
@@ -386,7 +386,7 @@ getModelTimeStep <- function(InputsModel) {
     return(NA)
   })
   TS <- TS[!is.na(TS)]
-  if(length(unique(TS)) != 1) {
+  if (length(unique(TS)) != 1) {
     stop("Time steps of the model of all nodes should be identical")
   }
   return(unique(TS))
@@ -401,7 +401,7 @@ getModelTimeStep <- function(InputsModel) {
 #' @return the selected column or value in respect to `id`
 #' @noRd
 getInputBV <- function(x, id, unset = NULL) {
-  if(is.null(x)) {
+  if (is.null(x)) {
     return(unset)
   }
   if (is.matrix(x) || is.data.frame(x)) {
@@ -412,7 +412,7 @@ getInputBV <- function(x, id, unset = NULL) {
     # vector (for ZInputs and NLayers)
     if (length(x) == 1 && is.null(names(x))) {
       return(x)
-    } else if(!id %in% names(x)) {
+    } else if (!id %in% names(x)) {
       return(unset)
     } else {
       return(x[id])
@@ -491,7 +491,7 @@ getNodeBasinArea <- function(i, griwrm) {
   if (i %in% which(Diversions)) return(NA)
   UpstreamNodeRows <-
     which(griwrm$down == griwrm$id[i] & !is.na(griwrm$down) & !Diversions)
-  if(length(UpstreamNodeRows) > 0) {
+  if (length(UpstreamNodeRows) > 0) {
     upstreamAreas <- sapply(UpstreamNodeRows, getNodeBasinArea, griwrm = griwrm)
     return(sum(upstreamAreas, na.rm = TRUE))
   } else {
diff --git a/R/CreateRunOptions.GRiwrmInputsModel.R b/R/CreateRunOptions.GRiwrmInputsModel.R
index 7b49faeab85ac0d216c7f85612644705449ef0e8..8bc88b9b2450b222eee899286ab9ec2bcfd39438 100644
--- a/R/CreateRunOptions.GRiwrmInputsModel.R
+++ b/R/CreateRunOptions.GRiwrmInputsModel.R
@@ -6,7 +6,7 @@ CreateRunOptions.GRiwrmInputsModel <- function(x, IniStates = NULL, ...) {
   RunOptions <- list()
   class(RunOptions) <- append(class(RunOptions), "GRiwrmRunOptions")
 
-  for(id in names(x)) {
+  for (id in names(x)) {
     RunOptions[[id]] <- CreateRunOptions(x[[id]], IniStates = IniStates[[id]], ...)
     RunOptions[[id]]$id <- id
   }
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index c74cc68539bd87d676e1d2428b9cc1ddd63c239a..ae1fe25edd9ddc791fcfbc3d5df07684b38ac221 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -37,7 +37,7 @@ CreateRunOptions.InputsModel <- function(x, ...) {
 
   # Add FUN_MOD in parameters if carried by InputsModel
   if (!"FUN_MOD" %in% names(dots)) {
-    if(!is.null(x$FUN_MOD)) {
+    if (!is.null(x$FUN_MOD)) {
       dots$FUN_MOD <- x$FUN_MOD
     } else {
       stop(" The parameter `FUN_MOD` must be defined")
diff --git a/R/CreateSupervisor.R b/R/CreateSupervisor.R
index 24c9f5b0d1d546cdc80683c3fa4e45b99662c5a0..f6da59b712e99856024ccfa114e1d1dbfb1a68b6 100644
--- a/R/CreateSupervisor.R
+++ b/R/CreateSupervisor.R
@@ -13,7 +13,7 @@
 #'
 #' @example man-examples/RunModel.Supervisor.R
 CreateSupervisor <- function(InputsModel, TimeStep = 1L) {
-  if(!inherits(InputsModel, "GRiwrmInputsModel")) {
+  if (!inherits(InputsModel, "GRiwrmInputsModel")) {
     stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateInputsModel.GRiwrm)")
   }
   if (!is.integer(TimeStep)) stop("`TimeStep` parameter must be an integer")
diff --git a/R/RunModel.GRiwrmInputsModel.R b/R/RunModel.GRiwrmInputsModel.R
index bd0f9833c69ca0be68aa9ad51006aa498ea41487..6fd9b6d0fdb271ab938720009651e21d7fec0a30 100644
--- a/R/RunModel.GRiwrmInputsModel.R
+++ b/R/RunModel.GRiwrmInputsModel.R
@@ -15,11 +15,11 @@ RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
   OutputsModel <- list()
   class(OutputsModel) <- c("GRiwrmOutputsModel", class(OutputsModel))
 
-  for(id in names(x)) {
+  for (id in names(x)) {
     message("RunModel.GRiwrmInputsModel: Processing sub-basin ", x[[id]]$id, "...")
 
     # Update x[[id]]$Qupstream with simulated upstream flows
-    if(any(x[[id]]$UpstreamIsModeled)) {
+    if (any(x[[id]]$UpstreamIsModeled)) {
       x[[id]] <- UpdateQsimUpstream(x[[id]], RunOptions[[id]], OutputsModel)
     }
     # Run the model for the sub-basin
diff --git a/R/RunModel.InputsModel.R b/R/RunModel.InputsModel.R
index 8fb7a70e4c1271bb9dd5077bc7d57cd071af732a..17f290794c72909dc327cb76194193f9438e602a 100644
--- a/R/RunModel.InputsModel.R
+++ b/R/RunModel.InputsModel.R
@@ -2,7 +2,14 @@
 #'
 #' @details This function calls [airGR::RunModel] (See [airGR::RunModel] for further details).
 #'
-#' The list produced by the function (See Value section of [airGR::RunModel_GR4J]) is here completed by an item *$Qsim_m3* storing the simulated discharge series in m3/s.
+#' The list produced by the function (See Value section of [airGR::RunModel_GR4J])
+#' is here completed by:
+#'
+#' - an item `$Qsim_m3` storing the simulated discharge series in m3/s
+#' - an item `$Qover_m3` storing the volumes of over abstraction which occurs
+#' when `RunModel_Lag` warns for negative simulated flows. The latter reflects the volume
+#' that was planned to be drawn from the sub-basin but could not be drawn because
+#' of the lack of water.
 #'
 #' @inheritParams airGR::RunModel
 #' @param x \[object of class \emph{InputsModel}\] see [airGR::CreateInputsModel] for details
@@ -29,7 +36,7 @@ RunModel.InputsModel <- function(x = NULL,
     }
   }
 
-  if(is.null(FUN_MOD)) {
+  if (is.null(FUN_MOD)) {
     if (x$isReservoir) {
       FUN_MOD <- "RunModel_Reservoir"
     } else {
@@ -39,26 +46,19 @@ RunModel.InputsModel <- function(x = NULL,
 
   FUN_MOD <- match.fun(FUN_MOD)
   if (identical(FUN_MOD, RunModel_Lag)) {
-    QcontribDown <- list(
-      RunOptions = list(
-        WarmUpQsim = rep(0, length(RunOptions$IndPeriod_WarmUp))
-      ),
-      Qsim = rep(0, length(RunOptions$IndPeriod_Run))
-    )
-    class(QcontribDown) <- c("OutputsModel", class(RunOptions)[-1])
-    x$BasinAreas[length(x$BasinAreas)] <- 1
-    OutputsModel <- RunModel_Lag(x, RunOptions, Param, QcontribDown)
-    OutputsModel$DatesR <- x$DatesR[RunOptions$IndPeriod_Run]
-  } else if((inherits(x, "GR") & is.null(x$UpstreamNodes)) | identical(FUN_MOD, RunModel_Reservoir)) {
+    OutputsModel <- RunModel.SD(x, RunOptions, Param)
+  } else if ((inherits(x, "GR") & is.null(x$UpstreamNodes)) | identical(FUN_MOD, RunModel_Reservoir)) {
     # Upstream basins and Reservoir are launch directly
     OutputsModel <- FUN_MOD(x, RunOptions, Param)
   } else {
-    # Intermediate basins (other than reservoir) are launch with SD capabilities
+    # Intermediate basins (other than reservoir) are launched with SD capabilities
     if (!is.null(x$UpstreamNodes) & !inherits(x, "SD")) {
       # For calibration of node with diversion
       class(x) <- c(class(x), "SD")
     }
     OutputsModel <- airGR::RunModel(x, RunOptions, Param, FUN_MOD)
+    OutputsModel <- calcOverAbstraction(OutputsModel, FALSE)
+    OutputsModel$RunOptions <- calcOverAbstraction(OutputsModel$RunOptions, TRUE)
   }
   OutputsModel$RunOptions$TimeStep <- RunOptions$FeatFUN_MOD$TimeStep
   if (is.null(OutputsModel$Qsim_m3)) {
diff --git a/R/RunModel.SD.R b/R/RunModel.SD.R
index 5f8b7816d913311e8d7223196c19f9c7c0115b8b..cb9853146c1660a49afe63ca2db7cd6a8e259abe 100644
--- a/R/RunModel.SD.R
+++ b/R/RunModel.SD.R
@@ -1,23 +1,38 @@
-#' Run a semi-distributed model from rainfall-runoff model outputs
+#' Run a hydraulic routing model from rainfall-runoff model outputs
 #'
 #' @inheritParams airGR::RunModel_Lag
 #' @param x \[object of class `InputsModel`\] used as `InputsModel` parameter for [airGR::RunModel_Lag]
 #' @param ... further arguments passed to or from other methods
 #'
 #' @return `OutputsModel` object. See [airGR::RunModel_Lag]
-#' @export
+#' @noRd
 #'
-RunModel.SD <- function(x, RunOptions, Param, QcontribDown, ...) {
-  if (x$isReservoir) {
-    OutputsModel <- RunModel_Reservoir(x,
-                                       RunOptions = RunOptions,
-                                       Param = Param[1:2])
-  } else {
-    OutputsModel <- airGR::RunModel_Lag(x,
-                                        RunOptions = RunOptions,
-                                        Param = Param[1],
-                                        QcontribDown = QcontribDown)
+RunModel.SD <- function(x, RunOptions, Param, QcontribDown = NULL, ...) {
+  if (is.null(QcontribDown)) {
+    QcontribDown <- list(
+      Qsim = rep(0, length(RunOptions$IndPeriod_Run))
+    )
+    if (!is.null(RunOptions$IndPeriod_WarmUp) &&
+        !(length(RunOptions$IndPeriod_WarmUp) == 1 && RunOptions$IndPeriod_WarmUp == 0L)) {
+      QcontribDown$RunOptions = list(WarmUpQsim = rep(0, length(RunOptions$IndPeriod_WarmUp)))
+    }
+    class(QcontribDown) <- c("OutputsModel", class(RunOptions)[-1])
+    x$BasinAreas[length(x$BasinAreas)] <- 1E-6
   }
+  OutputsModel <- airGR::RunModel_Lag(x,
+                                      RunOptions = RunOptions,
+                                      Param = Param[1],
+                                      QcontribDown = QcontribDown)
+  if (is.null(OutputsModel$DatesR)) {
+    OutputsModel$DatesR <- x$DatesR[RunOptions$IndPeriod_Run]
+  }
+  if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$RunOptions$WarmUpQsim_m3 <-
+      OutputsModel$RunOptions$WarmUpQsim * sum(x$BasinAreas, na.rm = TRUE) * 1e3
+  }
+  OutputsModel <- calcOverAbstraction(OutputsModel, FALSE)
+  OutputsModel$RunOptions <- calcOverAbstraction(OutputsModel$RunOptions, TRUE)
+
   OutputsModel$RunOptions$TimeStep <- RunOptions$FeatFUN_MOD$TimeStep
   return(OutputsModel)
 }
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index d6b20474312aa148cae6f8e64e17c3cbf3bee9e1..1c67d692cb5ced4adf8ac6e1e08ef6348b001821 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -33,7 +33,7 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
 
   # Run runoff model for each sub-basin
   x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
-    if(inherits(IM, "GR")) {
+    if (inherits(IM, "GR")) {
       OM_GR <- RunModel.GR(IM,
                            RunOptions = RunOptions[[IM$id]],
                            Param = Param[[IM$id]])
@@ -44,19 +44,19 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   class(x$OutputsModel) <- c("GRiwrmOutputsModel", class(x$OutputsModel))
 
   # Copy simulated pure runoff flows (no SD nor Diversion nodes) to Qupstream
-  for(id in getNoSD_Ids(x$InputsModel, include_diversion = FALSE)) {
+  for (id in getNoSD_Ids(x$InputsModel, include_diversion = FALSE)) {
     updateQupstream.Supervisor(x, id, IndPeriod_Run)
   }
 
-  # Store OutputsModel for step by step simulation
-  x$storedOutputs <- initStoredOutputs(x)
-
   # Initialization of model states by running the model with no supervision on warm-up period
   RunOptionsWarmUp <- RunOptions
-  for(id in names(x$InputsModel)) {
+  for (id in names(x$InputsModel)) {
     RunOptionsWarmUp[[id]]$IndPeriod_Run <- RunOptionsWarmUp[[id]]$IndPeriod_WarmUp
     RunOptionsWarmUp[[id]]$IndPeriod_WarmUp <- 0L
     RunOptionsWarmUp[[id]]$Outputs_Sim <- c("StateEnd", "Qsim")
+    if (x$InputsModel[[id]]$isReservoir) {
+      RunOptionsWarmUp[[id]]$Outputs_Sim <- c(RunOptionsWarmUp[[id]]$Outputs_Sim, "Qsim_m3")
+    }
   }
   OM_WarmUp <- suppressMessages(
     RunModel.GRiwrmInputsModel(x$InputsModel,
@@ -67,7 +67,7 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   # Adapt RunOptions to step by step simulation and copy states
   SD_Ids <- getSD_Ids(x$InputsModel, add_diversions = TRUE)
   names(SD_Ids) <- SD_Ids
-  for(id in SD_Ids) {
+  for (id in SD_Ids) {
     RunOptions[[id]]$IndPeriod_WarmUp <- 0L
     RunOptions[[id]]$Outputs_Sim <- c("Qsim_m3", "StateEnd")
     x$OutputsModel[[id]]$StateEnd <- serializeIniStates(OM_WarmUp[[id]]$StateEnd)
@@ -75,20 +75,23 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
 
   # Set Outputs to archive for final restitution
   outputVars <- lapply(SD_Ids, function(id) {
-    ov <- "Qsim_m3"
+    ov <- c("Qsim_m3", "Qover_m3")
     if (x$InputsModel[[id]]$hasDiversion) {
       ov <- c(ov, "Qdiv_m3", "Qnat")
-    } else if (x$InputsModel[[id]]$isReservoir) {
+    }
+    if (x$InputsModel[[id]]$isReservoir) {
       ov <- c(ov, "Qinflows_m3", "Vsim")
     }
     return(ov)
   })
+  # Store OutputsModel for step by step simulation
+  x$storedOutputs <- initStoredOutputs(x, outputVars)
 
   message("Processing: 0%", appendLF = FALSE)
   iProgressSteps <- round(length(lSuperTS) * seq(0.1, 0.9, 0.1))
 
   # Loop over time steps with a step equal to the supervision time step
-  for(i in seq_along(lSuperTS)) {
+  for (i in seq_along(lSuperTS)) {
     iProgressMessage <- which(i == iProgressSteps)
     if (length(iProgressMessage) == 1) {
       message(" ", 10 * iProgressMessage, "%", appendLF = FALSE)
@@ -98,24 +101,28 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
     x$ts.index <- iTS - x$ts.index0
     x$ts.date <- x$InputsModel[[1]]$DatesR[iTS]
     # Regulation occurs from second time step
-    if(iTS[1] > ts.start) {
+    if (iTS[1] > ts.start) {
       doSupervision(x)
     }
     # Loop over sub-basin using SD model
-    for(id in SD_Ids) {
+    for (id in SD_Ids) {
       # Run model for the sub-basin and one time step
       RunOptions[[id]]$IniStates <- serializeIniStates(x$OutputsModel[[id]]$StateEnd)
       RunOptions[[id]]$IndPeriod_Run <- iTS
-      if (RunOptions[[id]]$FeatFUN_MOD$IsSD) {
-        # Route upstream flows for SD nodes
+      # Route upstream flows for SD nodes
+      if (x$InputsModel[[id]]$isReservoir) {
+        x$OutputsModel[[id]] <- RunModel_Reservoir(
+          x$InputsModel[[id]],
+          RunOptions = RunOptions[[id]],
+          Param = Param[[id]]
+        )
+      } else {
         x$OutputsModel[[id]] <- RunModel.SD(
           x$InputsModel[[id]],
           RunOptions = RunOptions[[id]],
           Param = Param[[id]],
           QcontribDown = x$storedOutputs$QcontribDown[x$ts.index, id]
         )
-      } else {
-        x$OutputsModel[[id]]$Qsim_m3 <- x$storedOutputs$Qsim_m3[x$ts.index, id]
       }
       if (x$InputsModel[[id]]$hasDiversion) {
         # Compute diverted and simulated flows on Diversion nodes
@@ -136,7 +143,7 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
 
   message(" 100%")
 
-  for(id in SD_Ids) {
+  for (id in SD_Ids) {
     x$OutputsModel[[id]]$DatesR <- x$DatesR[IndPeriod_Run]
     for (outputVar in outputVars[[id]]) {
       x$OutputsModel[[id]][[outputVar]] <- x$storedOutputs[[outputVar]][, id]
@@ -152,9 +159,10 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   return(x$OutputsModel)
 }
 
+
 updateQupstream.Supervisor <- function(x, id, iTS) {
   downId <- x$InputsModel[[id]]$down
-  if(!is.null(x$InputsModel[[downId]])) {
+  if (!is.null(x$InputsModel[[downId]])) {
     x$InputsModel[[downId]]$Qupstream[iTS, id] <-
       x$OutputsModel[[id]]$Qsim_m3
   }
diff --git a/R/RunModel_Reservoir.R b/R/RunModel_Reservoir.R
index 3e20530576207bc9c1744900d5c410282e5530c6..1d4e0e52845390a0349726880403ac5f6e1062b3 100644
--- a/R/RunModel_Reservoir.R
+++ b/R/RunModel_Reservoir.R
@@ -51,7 +51,9 @@ RunModel_Reservoir <- function(InputsModel, RunOptions, Param) {
   celerity <- Param[2]
 
   # Compute inflows with RunModel_Lag
-  OutputsModel <- RunModel(InputsModel, RunOptions, celerity, FUN_MOD = "RunModel_Lag")
+  OutputsModel <- RunModel.SD(InputsModel,
+                              RunOptions,
+                              Param = celerity)
   names(OutputsModel)[names(OutputsModel) == "Qsim_m3"] <- "Qinflows_m3"
   Qinflows_m3 <- c(OutputsModel$RunOptions$WarmUpQsim_m3,
                    OutputsModel$Qinflows_m3)
@@ -73,7 +75,7 @@ RunModel_Reservoir <- function(InputsModel, RunOptions, Param) {
   }
 
   # Time series volume and release calculation
-  for(i in iPerTot) {
+  for (i in iPerTot) {
     Vsim[i] <- V0 + Qinflows_m3[i]
     if (InputsModel$hasDiversion) {
       Qdiv_m3[i] <- min(Vsim[i] + InputsModel$Qmin[IndPerTot[i]], InputsModel$Qdiv[IndPerTot[i]])
@@ -89,7 +91,7 @@ RunModel_Reservoir <- function(InputsModel, RunOptions, Param) {
   }
 
   # Format OutputsModel
-  if(length(IndPerWarmUp) > 0) {
+  if (length(IndPerWarmUp) > 0) {
     iWarmUp <- seq(length(RunOptions$IndPeriod_WarmUp))
     OutputsModel$RunOptions$WarmUpQsim_m3 <- Qsim_m3[iWarmUp]
     OutputsModel$RunOptions$WarmUpVsim <- Vsim[iWarmUp]
diff --git a/R/UpdateQsimUpstream.R b/R/UpdateQsimUpstream.R
index 9d0a012a2b4af6697babb3ddb838b9b2eb3eabd9..80148499ec3db3f10e09845305ef69007c41c6a6 100644
--- a/R/UpdateQsimUpstream.R
+++ b/R/UpdateQsimUpstream.R
@@ -12,7 +12,7 @@
 #'
 UpdateQsimUpstream <- function(InputsModel, Runoptions, OutputsModel) {
   iQ <- which(InputsModel$UpstreamIsModeled)
-  for(i in iQ) {
+  for (i in iQ) {
     InputsModel$Qupstream[Runoptions$IndPeriod_Run, i] <- OutputsModel[[InputsModel$UpstreamNodes[i]]][[InputsModel$UpstreamVarQ[i]]]
     varWarmupQ_m3 <- paste0("WarmUp", InputsModel$UpstreamVarQ[i])
     if (!is.null(OutputsModel[[InputsModel$UpstreamNodes[i]]]$RunOptions[[varWarmupQ_m3]])) {
diff --git a/R/globals.R b/R/globals.R
deleted file mode 100644
index 87c0cf2e82f44b5c8fdb4f696f5702907593f56b..0000000000000000000000000000000000000000
--- a/R/globals.R
+++ /dev/null
@@ -1,6 +0,0 @@
-#' getNodeRanking: no visible binding for global variable 'donor'
-#' updateParameters4Ungauged: no visible binding for global variable
-#' down'
-#' updateParameters4Ungauged: no visible binding for global variable
-#' model'
-utils::globalVariables(c("donor", "down ", "model "))
diff --git a/R/plot.GRiwrm.R b/R/plot.GRiwrm.R
index cddb2f90ac08552773ada031ebf6caa0eb605709..eadf1e624273ea862483a4041374d481628bef3b 100644
--- a/R/plot.GRiwrm.R
+++ b/R/plot.GRiwrm.R
@@ -59,7 +59,7 @@ plot.GRiwrm <- function(x,
     x$id[x$nodeclass == nc]
   })
   names(node_class) <- unique(x$nodeclass)
-  node_class <- lapply(node_class, function(id) if(length(id) > 0) paste0("id_", id))
+  node_class <- lapply(node_class, function(id) if (length(id) > 0) paste0("id_", id))
   node_class <- paste("class", sapply(node_class, paste, collapse = ","), names(node_class))
   css <- c(
     paste("classDef default", defaultClassDef),
diff --git a/R/plot.Qm3s.R b/R/plot.Qm3s.R
index 0b378cacdbfde6546b236de4378ddff7c77cc1d1..870a8acc8247c2706b16d522c140fac8ea0d077e 100644
--- a/R/plot.Qm3s.R
+++ b/R/plot.Qm3s.R
@@ -22,7 +22,6 @@
 #'
 #' @example man-examples/RunModel.GRiwrmInputsModel.R
 #'
-#' @importFrom grDevices rainbow
 #' @importFrom graphics matplot
 #' @export plot.Qm3s
 #' @export
@@ -32,7 +31,7 @@ plot.Qm3s <- function(x,
                       xlab = "Date",
                       ylab = expression("Flow rate (m"^"3"*"/s)"),
                       main = "Simulated flows",
-                      col = rainbow(ncol(x) - 1),
+                      col = grDevices::hcl.colors(ncol(x) - 1, "Zissou 1"),
                       legend = colnames(x)[-1],
                       legend.cex = 0.7,
                       legend.x = "topright",
@@ -53,7 +52,7 @@ plot.Qm3s <- function(x,
     main = main,
     col = col, ...
   )
-  if(!is.null(legend)) {
+  if (!is.null(legend)) {
     legend(x = legend.x,
            y = legend.y,
            legend = legend,
diff --git a/R/utils.CreateInputsModel.R b/R/utils.CreateInputsModel.R
index c5a4ab9481327ec8df388d96361540e45ea50f34..f473014d4d1c14cc4af43dbb63f65f823b0ab077 100644
--- a/R/utils.CreateInputsModel.R
+++ b/R/utils.CreateInputsModel.R
@@ -2,7 +2,7 @@ updateQObsQrelease <- function(g, Qobs, Qrelease) {
   reservoirIds <- g$id[!is.na(g$model) & g$model == "RunModel_Reservoir"]
   # Fill Qrelease with Qobs
   warn_ids <- NULL
-  for(id in reservoirIds) {
+  for (id in reservoirIds) {
     if (!id %in% colnames(Qrelease)) {
       if (id %in% colnames(Qobs)) {
         if (!any(g$id == id & (!is.na(g$model) & g$model == "Diversion"))) {
diff --git a/R/utils.RunModel.R b/R/utils.RunModel.R
index 054fa0fa8b75d88c2742ce33562515725676a8ea..fc3e58b203ea9c1fcafa606029c5a8d29d0e318e 100644
--- a/R/utils.RunModel.R
+++ b/R/utils.RunModel.R
@@ -7,9 +7,9 @@
 #' @param Param [list] of containing model parameter values of each node of the network
 #' @noRd
 checkRunModelParameters <- function(InputsModel, RunOptions, Param) {
-  if(!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateRunOptions.GRiwrmInputsModel)")
-  if(!inherits(RunOptions, "GRiwrmRunOptions")) stop("Argument `RunOptions` parameter must of class 'GRiwrmRunOptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
-  if(!is.list(Param) || !all(names(InputsModel) %in% names(Param))) stop("Argument `Param` must be a list with names equal to nodes IDs")
+  if (!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateRunOptions.GRiwrmInputsModel)")
+  if (!inherits(RunOptions, "GRiwrmRunOptions")) stop("Argument `RunOptions` parameter must of class 'GRiwrmRunOptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
+  if (!is.list(Param) || !all(names(InputsModel) %in% names(Param))) stop("Argument `Param` must be a list with names equal to nodes IDs")
 }
 
 
@@ -64,3 +64,27 @@ serializeIniStates <- function(IniStates) {
   unlist(IniStates)
 }
 
+
+#' Cap negative `OutputsModel$Qsim_m3` to zero and fill `OutputsModel$Qover_m3`
+#' with over-abstracted volumes
+#'
+#' @param O Either `OutputsModel` or `OutputsModel$RunOptions` (for warm-up Qsim)
+#' @param WarmUp `TRUE` if `O` is `OutputsModel$RunOptions`
+#'
+#' @return Modified `OutputsModel` or `OutputsModel$RunOptions`
+#' @noRd
+#'
+calcOverAbstraction <- function(O, WarmUp) {
+  f <- list(sim = "Qsim_m3", over = "Qover_m3")
+  if (WarmUp) {
+    f <- lapply(f, function(x) paste0("WarmUp", x))
+  }
+  if (!is.null(O[[f$sim]])) {
+    O[[f$over]] <- rep(0, length(O[[f$sim]]))
+    if (any(!is.na(O[[f$sim]]) & O[[f$sim]] < 0)) {
+      O[[f$over]][O[[f$sim]] < 0] <- - O[[f$sim]][!is.na(O[[f$sim]]) & O[[f$sim]] < 0]
+      O[[f$sim]][!is.na(O[[f$sim]]) & O[[f$sim]] < 0] <- 0
+    }
+  }
+  return(O)
+}
diff --git a/R/utils.Supervisor.R b/R/utils.Supervisor.R
index 4b80681c0a7dfd50ab2ddc47e7d1a64bf68c3728..8856aba41939d26075dcb90961b25dd96b54a5d5 100644
--- a/R/utils.Supervisor.R
+++ b/R/utils.Supervisor.R
@@ -43,7 +43,7 @@ setDataToLocation <- function(ctrlr, sv) {
       # Direct injection node => update Qusptream of downstream node
       node <- sv$griwrm4U$down[sv$griwrm4U$id == locU]
       # ! Qupstream contains warm up period and run period => the index is shifted
-      if(!is.null(sv$InputsModel[[node]])) {
+      if (!is.null(sv$InputsModel[[node]])) {
         sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index, locU] <- U
       }
     } else if (sv$nodeProperties[locU, "Diversion"]){
@@ -73,7 +73,7 @@ doSupervision <- function(supervisor) {
     # Run logic
     supervisor$controllers[[id]]$U <-
       supervisor$controllers[[id]]$FUN(supervisor$controllers[[id]]$Y)
-    if(is.vector(supervisor$controllers[[id]]$U)) {
+    if (is.vector(supervisor$controllers[[id]]$U)) {
       supervisor$controllers[[id]]$U <- matrix(supervisor$controllers[[id]]$U, nrow = 1)
     }
     # Check U output
@@ -97,30 +97,21 @@ doSupervision <- function(supervisor) {
 }
 
 
-initStoredOutputs <- function(x) {
-  np <- getAllNodesProperties(x$griwrm)
-  so <- list()
-  so$QcontribDown <- do.call(
+initStoredOutputs <- function(x, outputVars) {
+  QcontribDown <- do.call(
     cbind,
     lapply(x$OutputsModel, "[[", "Qsim")
   )
-  so$Qsim_m3 <- do.call(
-    cbind,
-    lapply(x$OutputsModel, "[[", "Qsim_m3")
-  )
-  if (sum(np$Diversion) > 0) {
-    # Outputs of Diversion nodes
-    so$Qdiv_m3 <- so$Qsim_m3[, np$id[np$Diversion], drop = FALSE] * NA
-    so$Qnat <- so$Qdiv_m3
-  }
-  if (sum(np$Reservoir) > 0) {
-    # Specific Outputs of RunModel_Reservoir
-    so$Vsim <- matrix(rep(NA, sum(np$Reservoir) * nrow(so$Qsim_m3)),
-                      nrow = nrow(so$Qsim_m3))
-    colnames(so$Vsim) <- np$id[np$Reservoir]
-    so$Qinflows_m3 <- so$Vsim
-    # Add columns Qsim_m3 at reservoir (out of the scope of GR models calculated above)
-    so$Qsim_m3 <- cbind(so$Qsim_m3, so$Vsim)
-  }
+  so <- lapply(setNames(nm = unique(unlist(outputVars))), function(ov) {
+    s <- sapply(outputVars, function(y) "Qsim_m3" %in% y)
+    ids <- names(s)[s]
+    if (length(ids) > 0) {
+      m <- matrix(NA, nrow = nrow(QcontribDown), ncol = length(ids))
+      colnames(m) <- ids
+      return(m)
+    }
+    return(NULL)
+  })
+  so$QcontribDown <- QcontribDown
   return(so)
 }
diff --git a/data-raw/Severn.R b/data-raw/Severn.R
index db720be1031c9481c2b04cb91611cba091d76c42..ad1cd3384152d03dc365eaaafa5393fe3231e97c 100644
--- a/data-raw/Severn.R
+++ b/data-raw/Severn.R
@@ -13,7 +13,7 @@ files <- file.path(
 )
 
 read_CAMELS_ts <- function(file, cut) {
-  if(!file.exists(file)) {
+  if (!file.exists(file)) {
     stop("File not found: ", file, "\n",
          "Please download it from:\n",
          "https://catalogue.ceh.ac.uk/documents/8344e4f3-d2ea-44f5-8afa-86d2987543a9")
diff --git a/man/RunModel.InputsModel.Rd b/man/RunModel.InputsModel.Rd
index 7b4887d9aadddbbe3a6a3e2415c8de1f00903461..32bec4a7217d8aab658e945021367073d58bd182 100644
--- a/man/RunModel.InputsModel.Rd
+++ b/man/RunModel.InputsModel.Rd
@@ -30,5 +30,13 @@ Wrapper for \link[airGR:RunModel]{airGR::RunModel} for one sub-basin
 \details{
 This function calls \link[airGR:RunModel]{airGR::RunModel} (See \link[airGR:RunModel]{airGR::RunModel} for further details).
 
-The list produced by the function (See Value section of \link[airGR:RunModel_GR4J]{airGR::RunModel_GR4J}) is here completed by an item \emph{$Qsim_m3} storing the simulated discharge series in m3/s.
+The list produced by the function (See Value section of \link[airGR:RunModel_GR4J]{airGR::RunModel_GR4J})
+is here completed by:
+\itemize{
+\item an item \verb{$Qsim_m3} storing the simulated discharge series in m3/s
+\item an item \verb{$Qover_m3} storing the volumes of over abstraction which occurs
+when \code{RunModel_Lag} warns for negative simulated flows. The latter reflects the volume
+that was planned to be drawn from the sub-basin but could not be drawn because
+of the lack of water.
+}
 }
diff --git a/man/RunModel.SD.Rd b/man/RunModel.SD.Rd
deleted file mode 100644
index 161bb4722a7696ae3ef7f82c93b0fcd8aabe39b8..0000000000000000000000000000000000000000
--- a/man/RunModel.SD.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RunModel.SD.R
-\name{RunModel.SD}
-\alias{RunModel.SD}
-\title{Run a semi-distributed model from rainfall-runoff model outputs}
-\usage{
-\method{RunModel}{SD}(x, RunOptions, Param, QcontribDown, ...)
-}
-\arguments{
-\item{x}{[object of class \code{InputsModel}] used as \code{InputsModel} parameter for \link[airGR:RunModel_Lag]{airGR::RunModel_Lag}}
-
-\item{RunOptions}{[object of class \emph{RunOptions}] see \code{\link[airGR]{CreateRunOptions}} for details}
-
-\item{Param}{[numeric] vector of 1 parameter
-    \tabular{ll}{
-      Velocity \tab mean flow velocity [m/s]
-    }
-  }
-
-\item{QcontribDown}{[numeric] vector or [OutputsModel] containing the time series of the runoff contribution of the downstream sub-basin}
-
-\item{...}{further arguments passed to or from other methods}
-}
-\value{
-\code{OutputsModel} object. See \link[airGR:RunModel_Lag]{airGR::RunModel_Lag}
-}
-\description{
-Run a semi-distributed model from rainfall-runoff model outputs
-}
diff --git a/man/plot.Qm3s.Rd b/man/plot.Qm3s.Rd
index 74f969ca3e86e14e0d7f5894913b89df0d851914..6d2f87b8d446d87f39f7e2ba30a01b923893b0d0 100644
--- a/man/plot.Qm3s.Rd
+++ b/man/plot.Qm3s.Rd
@@ -10,7 +10,7 @@
   xlab = "Date",
   ylab = expression("Flow rate (m"^"3" * "/s)"),
   main = "Simulated flows",
-  col = rainbow(ncol(x) - 1),
+  col = grDevices::hcl.colors(ncol(x) - 1, "Zissou 1"),
   legend = colnames(x)[-1],
   legend.cex = 0.7,
   legend.x = "topright",
diff --git a/tests/testthat/helper_1_RunModel.R b/tests/testthat/helper_1_RunModel.R
index ef1d0b7cf85fa9ebf32c16f8739a224e2616a20d..07ae1a53ed0e69cff32824f7caf1a3ecb5a43574 100644
--- a/tests/testthat/helper_1_RunModel.R
+++ b/tests/testthat/helper_1_RunModel.R
@@ -8,7 +8,7 @@
 #' e <- setupRunModel()
 # variables are copied from environment 'e' to the current environment
 # https://stackoverflow.com/questions/9965577/r-copy-move-one-environment-to-another
-#' for(x in ls(e)) assign(x, get(x, e))
+#' for (x in ls(e)) assign(x, get(x, e))
 #'
 setupRunModel <-
   function(runInputsModel = TRUE,
@@ -36,7 +36,7 @@ setupRunModel <-
     }))
 
     # Set network
-    if(is.null(griwrm)) {
+    if (is.null(griwrm)) {
       nodes <- loadSevernNodes()
       griwrm <-
         CreateGRiwrm(nodes)
diff --git a/tests/testthat/helper_Calibration.R b/tests/testthat/helper_Calibration.R
index 63ed941664de80ec0929522cd682cc5182e09a42..e4dba204170b2aab5b8afb9f2a4481c0a87041fd 100644
--- a/tests/testthat/helper_Calibration.R
+++ b/tests/testthat/helper_Calibration.R
@@ -16,7 +16,7 @@ runCalibration <- function(nodes = NULL,
                      runRunModel = runRunModel,
                      Qobs2 = Qobs2,
                      IsHyst = IsHyst)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   rm(e)
   np <- getAllNodesProperties(griwrm)
 
diff --git a/tests/testthat/helper_RunModel_Reservoir.R b/tests/testthat/helper_RunModel_Reservoir.R
index 48a8167df7815fd1a9e76cefa82aad1d7c3aa5f0..45bd3c7b4e2b80e3de3ac1f3d0f01d9230ce6175 100644
--- a/tests/testthat/helper_RunModel_Reservoir.R
+++ b/tests/testthat/helper_RunModel_Reservoir.R
@@ -69,7 +69,7 @@ testDerivedUngauged <- function(donorByDerivation) {
   CalibOptions <- CreateCalibOptions(InputsModel,
                                      FixedParam = list(Dam = c(650E6, 1)))
   e <- runCalibration(g, Qobs2 = Qobs2, CalibOptions = CalibOptions)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   expect_equal(Param[["54095"]][1:3],
                Param[[ifelse(donorByDerivation, "54029", "54001")]][2:4])
 }
diff --git a/tests/testthat/test-CreateInputsModel.R b/tests/testthat/test-CreateInputsModel.R
index fa6acf9d2a2575f7ebcadebd14c8e00929b125f4..9c67fd18f2a16dcf1bde88cd36cc5c4b56f7ba22 100644
--- a/tests/testthat/test-CreateInputsModel.R
+++ b/tests/testthat/test-CreateInputsModel.R
@@ -266,7 +266,7 @@ test_that("Node with upstream nodes having area = NA should return correct Basin
     Dam = rep(0,11536)
   )
   e <- setupRunModel(griwrm = g, runInputsModel = FALSE, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   InputsModel <-
     suppressWarnings(CreateInputsModel(g, DatesR, Precip, PotEvap, Qobs = Qobs2))
   expect_equal(sum(InputsModel$`54001`$BasinAreas),
@@ -276,7 +276,7 @@ test_that("Node with upstream nodes having area = NA should return correct Basin
 test_that("Use of Qobs for Qrelease should raise a warning",  {
   g <- CreateGRiwrm(n_rsrvr)
   e <- setupRunModel(griwrm = g, runInputsModel = FALSE)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   expect_warning(CreateInputsModel(griwrm, DatesR, Precip, PotEvap,
                                    TempMean = TempMean,
                                    Qobs = Qobs_rsrvr))
diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R
index d3f39504b41aa91046e89fda5fd958e149742467..cdd1d675d49ed37dc2b1070d4f88653e07745120 100644
--- a/tests/testthat/test-RunModel.R
+++ b/tests/testthat/test-RunModel.R
@@ -236,7 +236,7 @@ test_that("A transparent upstream node with area=NA should return same result #1
   )
   g <- CreateGRiwrm(nodes)
   e <- setupRunModel(griwrm = g, runRunModel = FALSE)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   Param <- ParamMichel[c("54095", "54001")]
   Param[["P"]] <- 1
   OM <- RunModel(InputsModel,
@@ -246,3 +246,23 @@ test_that("A transparent upstream node with area=NA should return same result #1
   names(OM[["54001"]]$StateEnd$SD[[1]]) <- "54095" # For matching upstream IDs with ref
   expect_equal(OM[["54001"]], OM_GriwrmInputs[["54001"]])
 })
+
+test_that("RunModel should return water deficit (Qover_m3)", {
+  nodes <- loadSevernNodes()
+  nodes <- nodes[nodes$id %in% c("54095", "54001"), ]
+  nodes[nodes$id == "54001", c("down", "length")] <- c(NA, NA)
+  nodes <- rbind(
+    nodes,
+    data.frame(id = "P", down = "54001", length = 10, area = NA, model = NA)
+  )
+  g <- CreateGRiwrm(nodes)
+  Qobs2 <- data.frame(P = rep(-2E6, length(DatesR)))
+  expect_warning(e <- setupRunModel(griwrm = g, runRunModel = TRUE, Qobs2 = Qobs2))
+  for (x in ls(e)) assign(x, get(x, e))
+  expect_false(any(OM_GriwrmInputs$`54001`$Qsim_m3 < 0))
+  expect_true(all(OM_GriwrmInputs$`54001`$Qover_m3 >= 0))
+  sv <- CreateSupervisor(InputsModel)
+  OM_sv <- RunModel(sv, RunOptions, ParamMichel)
+  expect_equal(OM_sv$`54001`$Qsim_m3, OM_GriwrmInputs$`54001`$Qsim_m3)
+  expect_equal(OM_sv$`54001`$Qover_m3, OM_GriwrmInputs$`54001`$Qover_m3)
+})
diff --git a/tests/testthat/test-RunModel.Supervisor.R b/tests/testthat/test-RunModel.Supervisor.R
index 8f8675cb8d24d8e3aeb4fc63e3dd5a8b21b029a7..06916956a4dff5410ef601f26ac6def968e0056f 100644
--- a/tests/testthat/test-RunModel.Supervisor.R
+++ b/tests/testthat/test-RunModel.Supervisor.R
@@ -94,7 +94,7 @@ test_that("RunModel.Supervisor with diversion node should not produce NAs", {
   Qobs2 <- matrix(data = rep(0, length(DatesR)), ncol = 1)
   colnames(Qobs2) <- "54001"
   e <- setupRunModel(griwrm = g_div, runRunModel = FALSE, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   sv <- CreateSupervisor(InputsModel, TimeStep = 1L)
   logicFunFactory <- function(sv) {
     #' @param Y Flow measured at "54002" the previous time step
diff --git a/tests/testthat/test-RunModel_Ungauged.R b/tests/testthat/test-RunModel_Ungauged.R
index 7ee25b3865e7534298030854f5e59d95d60b3032..4d8d88c802f6c78d859e5611b1eb1008e6f23259 100644
--- a/tests/testthat/test-RunModel_Ungauged.R
+++ b/tests/testthat/test-RunModel_Ungauged.R
@@ -8,7 +8,7 @@ test_that("RunModel_Ungauged should act as RunModel", {
   nodes$model[nodes$id == "54095"] <- "Ungauged"
   g <- CreateGRiwrm(nodes)
   e <- setupRunModel(runRunModel = FALSE, griwrm = g)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
 
   Param <- ParamMichel["54001"]
   Param[["54095"]] <-
@@ -65,7 +65,7 @@ test_that("RunModel_Ungauged works with a diversion as donor (#110)", {
   Qobs2 <- matrix(0, ncol = 1, nrow = 11536)
   colnames(Qobs2) <- "54032"
   e <- runCalibration(nodes, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   OCdiv <- OutputsCalib
   expect_equal(OCdiv, OC)
 })
@@ -101,7 +101,7 @@ test_that("RunModel_Ungauged works with a diversion as upstream node (#113)", {
   Qobs2 <- matrix(0, ncol = 1, nrow = 11536)
   colnames(Qobs2) <- "54095"
   e <- runCalibration(nodes, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   expect_equal(OutputsCalib$`54032`$CritFinal, CritValue)
 })
 
@@ -111,7 +111,7 @@ test_that("RunModel_Ungauged works with a diversion as upstream node (#113)", {
   Qobs2 <- matrix(0, ncol = 1, nrow = 11536)
   colnames(Qobs2) <- "54095"
   e <- runCalibration(nodes, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   expect_equal(OutputsCalib$`54032`$CritFinal, CritValue)
 })
 
@@ -123,7 +123,7 @@ test_that("Ungauged node with diversion outside the sub-network should work", {
 
   # First without Diversion
   e <- runCalibration(nodes, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   OC1 <- OutputsCalib
   Param1 <- Param
   OM1 <- RunModel(
@@ -146,7 +146,7 @@ test_that("Ungauged node with diversion outside the sub-network should work", {
   Qobs2 <- matrix(0, ncol = 1, nrow = 11536)
   colnames(Qobs2) <- "54095"
   e <- runCalibration(nodes, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   OC2 <- OutputsCalib
   expect_equal(OC2$`54001`$CritFinal, OC1$`54001`$CritFinal)
   expect_equal(OC2$`54032`$CritFinal, OC1$`54032`$CritFinal)
@@ -184,7 +184,7 @@ test_that("Ungauged node with upstream node with diversion should work", {
   Qobs2 <- matrix(0, ncol = length(g$id[g$model == "Diversion"]), nrow = 11536)
   colnames(Qobs2) <- g$id[g$model == "Diversion"]
   e <- setupRunModel(griwrm = g, runRunModel = FALSE, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
 
   ParamRef[["P"]] <- 1
   OM <- RunModel(InputsModel,
@@ -230,7 +230,7 @@ test_that("Donor node with diversion should work", {
   Qobs2 <- matrix(0, ncol = length(g$id[g$model == "Diversion"]), nrow = 11536)
   colnames(Qobs2) <- g$id[g$model == "Diversion"]
   e <- runCalibration(nodes, Qobs2 = Qobs2)
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   expect_equal(OC_ref$`54032`$CritFinal, OutputsCalib$`54032`$CritFinal, tolerance = 1E-3)
 })
 
@@ -261,7 +261,7 @@ test_that("Cemaneige with hysteresis works",  {
   e <- suppressWarnings(
     setupRunModel(griwrm = griwrm, runRunModel = FALSE, IsHyst = TRUE)
   )
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
 
   expect_true(all(sapply(InputsModel, function(x) x$model$hasX4)))
 
@@ -285,7 +285,7 @@ test_that("Cemaneige with hysteresis works",  {
   e <- suppressWarnings(
     runCalibration(nodes, InputsCrit = InputsCrit, CalibOptions = CO, IsHyst = TRUE)
   )
-  for(x in ls(e)) assign(x, get(x, e))
+  for (x in ls(e)) assign(x, get(x, e))
   expect_equal(sapply(Param, length),
                c("54057" = 9, "54032" = 9, "54001" = 8))
 })
diff --git a/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd b/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
index 0d761d1bf9d70b3ad716a10646a54fbec0634ada..df821f4d62ae59a4946667ecc33258987a5adbf9 100644
--- a/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
+++ b/vignettes/V04_Closed-loop_regulated_withdrawal.Rmd
@@ -257,7 +257,7 @@ As for the previous model, we need to set up an `GRiwrmInputsModel` object conta
 ```{r}
 # Flow time series are needed for all direct injection nodes in the network
 # even if they may be overwritten after by a controller
-QobsIrrig <- data.frame(Irrigation1 = rep(0, length(DatesR)), 
+QobsIrrig <- data.frame(Irrigation1 = rep(0, length(DatesR)),
                         Irrigation2 = rep(0, length(DatesR)))
 
 # Creation of the GRiwrmInputsModel object
@@ -298,7 +298,7 @@ fIrrigationFactory <- function(supervisor,
     month <- as.numeric(format(supervisor$ts.date, "%m"))
     U <- irrigationObjective[month, c(2,3)] # m3/s
     meanU <- mean(rowSums(U))
-    if(meanU > 0) {
+    if (meanU > 0) {
       # calculate the naturalized flow from the measured flow and the abstracted flow of the previous week
       lastU <- supervisor$controllers[[supervisor$controller.id]]$U # m3/day
       Qnat <- (Y - rowSums(lastU)) / 86400 # m3/s
diff --git a/vignettes/V05_Modelling_ungauged_nodes.Rmd b/vignettes/V05_Modelling_ungauged_nodes.Rmd
index 190fedfefd03af1d93c2f9d90d88933009655c99..6b0cafc6524b92eb11bcd3d9d717003ca7b8359c 100644
--- a/vignettes/V05_Modelling_ungauged_nodes.Rmd
+++ b/vignettes/V05_Modelling_ungauged_nodes.Rmd
@@ -30,19 +30,19 @@ Ungauged nodes in the semi-distributed model can be used to reach two different
 - increase spatial resolution of the rain fall to improve streamflow simulation [@lobligeoisWhenDoesHigher2014].
 - simulate streamflows in location of interest for management purpose
 
-This vignette introduces the implementation in airGRiwrm of the method developped by @lobligeoisMieuxConnaitreDistribution2014 for calibrating ungauged nodes in a 
+This vignette introduces the implementation in airGRiwrm of the method developped by @lobligeoisMieuxConnaitreDistribution2014 for calibrating ungauged nodes in a
 semi-distributed model.
 
 ## Presentation of the study case
 
-Using the study case of the vignette #1 and #2, we considere this time that nodes `54001` and 
-`54029` are ungauged. We simulate the streamflow at these locations by sharing 
+Using the study case of the vignette #1 and #2, we considere this time that nodes `54001` and
+`54029` are ungauged. We simulate the streamflow at these locations by sharing
 hydrological parameters of the gauged node `54032`.
 
 ```{r network, echo = FALSE}
 mmd <- function(x, ...) {
   # For avoiding crash of R CMD build in console mode
-  if(Sys.getenv("RSTUDIO") == "1") {
+  if (Sys.getenv("RSTUDIO") == "1") {
     DiagrammeR::mermaid(x, ...)
   }
 }
@@ -157,7 +157,7 @@ OC_U <- suppressWarnings(
   Calibration(IM_U, RunOptions, InputsCrit, CalibOptions))
 ```
 
-Hydrological parameters for sub-basins 
+Hydrological parameters for sub-basins
 
 ## Run the model with the optimized model parameters
 
@@ -166,9 +166,9 @@ The hydrologic model uses parameters herited from the calibration of the gauged
 ```{r param}
 ParamV05 <- sapply(griwrmV05$id, function(x) {OC_U[[x]]$Param})
 dfParam <- do.call(
-  rbind, 
-  lapply(ParamV05, function(x) 
-    if(length(x)==4) {return(c(NA, x))} else return(x))
+  rbind,
+  lapply(ParamV05, function(x)
+    if (length(x)==4) {return(c(NA, x))} else return(x))
 )
 colnames(dfParam) <- c("velocity", paste0("X", 1:4))
 knitr::kable(round(dfParam, 3))
diff --git a/vignettes/airGRiwrm.bib b/vignettes/airGRiwrm.bib
index 03dc558861d32041e9eb745a3582e294c2bd67f8..914ff752d0862550a64039ed08f051c307f3f774 100644
--- a/vignettes/airGRiwrm.bib
+++ b/vignettes/airGRiwrm.bib
@@ -213,7 +213,7 @@
   author = {Lobligeois, Florent},
   year = {2014},
   month = mar,
-  url = {https://pastel.archives-ouvertes.fr/tel-01134990/document},
+  url = {https://pastel.hal.science/tel-01134990/document},
   urldate = {2017-04-06},
   school = {AgroParisTech},
 }
diff --git a/vignettes/seinebasin/V02_First_run.Rmd b/vignettes/seinebasin/V02_First_run.Rmd
index 6f3bc3435cc97f4f76ae972bf4fc429de3f1c29b..a774a184a9c0bd3085de8cbc2e438cf52f78a4a3 100644
--- a/vignettes/seinebasin/V02_First_run.Rmd
+++ b/vignettes/seinebasin/V02_First_run.Rmd
@@ -56,7 +56,7 @@ ParamClimAware <- sapply(griwrm$id, function(id) {
   Param <- unlist(nodeParam[c("X1", "X2", "X3", "X4")])
   # Add lag model parameter if upstream nodes exist
   UpstrNodes <- which(griwrm$down == id & !is.na(griwrm$down))
-  if(length(UpstrNodes) > 0) {
+  if (length(UpstrNodes) > 0) {
     maxLength <- max(griwrm$length[UpstrNodes])
     Param <- c(
       maxLength * 1000 / ((nodeParam$Tau0 + nodeParam$K0) * 3600),
diff --git a/vignettes/seinebasin/V04_Open-loop_influenced_flow.Rmd b/vignettes/seinebasin/V04_Open-loop_influenced_flow.Rmd
index a8aa3d08210333d167b29f55904f827d4db79259..baefcf70cbbaf7b37b09d05fc9fcb6383a44abaf 100644
--- a/vignettes/seinebasin/V04_Open-loop_influenced_flow.Rmd
+++ b/vignettes/seinebasin/V04_Open-loop_influenced_flow.Rmd
@@ -23,14 +23,14 @@ library(airGRiwrm)
 
 ## Integration of the reservoir connections into the model
 
-The aim of this vignette is to add reservoir connections to the Seine River model 
-previously set up. A complete description of the reservoir system of the Seine 
+The aim of this vignette is to add reservoir connections to the Seine River model
+previously set up. A complete description of the reservoir system of the Seine
 River can be found in @dehayEtudeImpactChangement2012
 
 ### Add connections to the reservoirs in the gauging station network
 
-Because of the **airGR** SD model structure, the topology of the network will 
-differ between the physical one and the modeled one. In the following, we provide 
+Because of the **airGR** SD model structure, the topology of the network will
+differ between the physical one and the modeled one. In the following, we provide
 details for each of the 4 lakes, and then we present the complete SD network.
 
 First, the physical topology of the Aube lake is represented below:
@@ -39,7 +39,7 @@ First, the physical topology of the Aube lake is represented below:
 
 mmd <- function(x, ...) {
   # For avoiding crash of R CMD build in console mode
-  if(Sys.getenv("RSTUDIO") == "1") {
+  if (Sys.getenv("RSTUDIO") == "1") {
     mermaid(x, ...)
   }
 }
@@ -56,8 +56,8 @@ style AUBE fill:#ccf
 ", width = "100%")
 ```
 
-In the SD model, we do not simulate intermediate flows as well as the reservoir 
-in the catchment. As a result, an equivalent topology compatible with **airGR** 
+In the SD model, we do not simulate intermediate flows as well as the reservoir
+in the catchment. As a result, an equivalent topology compatible with **airGR**
 will be the one below:
 
 ```{r mmd_aube_model, echo = FALSE}
@@ -71,7 +71,7 @@ style AUBE_R3 fill:#fcc
 ", width = "100%")
 ```
 
-`AUBE_P2` will propagate negative flows downstream (reservoir inflows) while 
+`AUBE_P2` will propagate negative flows downstream (reservoir inflows) while
 `AUBE_R3` will propagate positive flows downstream (reservoir releases).
 
 The configuration on the lake Seine is similar:
@@ -115,8 +115,8 @@ style P fill:#ccf
 ", width = "100%")
 ```
 
-We can keep the same structure to model it. `PANNEC_R` corresponds to the flow 
-released by the Pannecière lake, it is acting as an upstream node which means 
+We can keep the same structure to model it. `PANNEC_R` corresponds to the flow
+released by the Pannecière lake, it is acting as an upstream node which means
 that the flow simulated at `CHAUM_07` is no longer routed to downstream.
 
 ```{r mmd_pannec_model, echo = FALSE}
@@ -194,7 +194,7 @@ plot(griwrm2)
 
 ## Loading reservoir observation time series
 
-Description of the files, the columns and the type of connection (inlet / outlet) 
+Description of the files, the columns and the type of connection (inlet / outlet)
 are defined in the list below:
 
 ```{r lCfgReservoirs}
@@ -230,7 +230,7 @@ load("_cache/V03.RData")
 
 ### How to handle former upstream sub-basins with upstream flows ?
 
-A lag parameter is now mandatory for these sub-basins. As no calibration is 
+A lag parameter is now mandatory for these sub-basins. As no calibration is
 possible at that stage an arbitrary one will be used (1 m/s).
 
 ```{r ParamMichel}