diff --git a/R/RunModel.InputsModel.R b/R/RunModel.InputsModel.R
index 8fb7a70e4c1271bb9dd5077bc7d57cd071af732a..3f687095fadb22301a5262471a5e57087e251609 100644
--- a/R/RunModel.InputsModel.R
+++ b/R/RunModel.InputsModel.R
@@ -2,7 +2,12 @@
 #'
 #' @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.
 #'
 #' @inheritParams airGR::RunModel
 #' @param x \[object of class \emph{InputsModel}\] see [airGR::CreateInputsModel] for details
@@ -73,6 +78,8 @@ RunModel.InputsModel <- function(x = NULL,
   if (x$hasDiversion && !x$isReservoir) {
     OutputsModel <- RunModel_Diversion(x, RunOptions, OutputsModel)
   }
+  OutputsModel <- calcOverAbstraction(OutputsModel, FALSE)
+  OutputsModel$RunOptions <- calcOverAbstraction(OutputsModel$RunOptions, TRUE)
   return(OutputsModel)
 }
 
@@ -135,3 +142,27 @@ calc_Qdiv<- function(Qnat, Qdiv, Qmin) {
   }
   return(list(Qsim = Qsim, Qdiv = Qnat - Qsim))
 }
+
+
+#' 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 (any(O[[f$sim]] < 0)) {
+    O[[f$over]] <- rep(0, length(O[[f$sim]]))
+    O[[f$over]][O[[f$sim]] < 0] <- - O[[f$sim]][O[[f$sim]] < 0]
+    O[[f$sim]][O[[f$sim]] < 0] <- 0
+  }
+  return(O)
+}
+
diff --git a/man/RunModel.InputsModel.Rd b/man/RunModel.InputsModel.Rd
index 7b4887d9aadddbbe3a6a3e2415c8de1f00903461..f4afb51d13e983abdda889f92d4c8cc84e922118 100644
--- a/man/RunModel.InputsModel.Rd
+++ b/man/RunModel.InputsModel.Rd
@@ -30,5 +30,11 @@ 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.
+}
 }
diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R
index d3f39504b41aa91046e89fda5fd958e149742467..ba5055ede2aefc5dc6104139611bd3adbd5e2840 100644
--- a/tests/testthat/test-RunModel.R
+++ b/tests/testthat/test-RunModel.R
@@ -246,3 +246,19 @@ 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)))
+  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))
+})