diff --git a/R/CreateInputsModel.R b/R/CreateInputsModel.R
index 2a2a6488620aaacdc33a124b7e4756064e2ce79f..3bd0f14752105a78f3010a3d8cd744d646f6fa2e 100644
--- a/R/CreateInputsModel.R
+++ b/R/CreateInputsModel.R
@@ -210,7 +210,7 @@ CreateInputsModel <- function(FUN_MOD,
       stop("'Qupstream' must have same number of rows as 'DatesR' length")
     }
     if(any(is.na(Qupstream))) {
-      stop("'Qupstream' cannot contain any NA value")
+      warning("'Qupstream' contains NA values: model outputs will contain NAs")
     }
     if(any(LengthHydro > 1000)) {
       warning("The unit of 'LengthHydro' has changed from m to km in v1.7 of airGR: values superior to 1000 km seem unrealistic")
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 9344a78bee11963b6da67520e7186e97ef94ce03..e98e1b867340deacf0597364dedc981bcb46249e 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -83,10 +83,16 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
       Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
       Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
   }
-  # Warning for negative flows
-  if (any(OutputsModel$Qsim < 0)) {
-    warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
-    OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+  # Warning for negative flows or NAs only in extended outputs
+  if(length(RunOptions$Outputs_Sim <3)) {
+    if (any(OutputsModel$Qsim[!is.na(OutputsModel$Qsim)] < 0)) {
+      warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
+      OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+    }
+    # Warning for NAs
+    if (any(is.na(OutputsModel$Qsim))) {
+      warning(length(which(is.na(OutputsModel$Qsim))), " time steps with NA values")
+    }
   }
   # Convert back Qsim to mm
   OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
diff --git a/tests/testthat/test-RunModel_Lag.R b/tests/testthat/test-RunModel_Lag.R
index 179bd87e7f45cd93705ef4a2f6902c4195bfeca7..b3eca0417b91edfb8ba46d30801ff449dc3a3509 100644
--- a/tests/testthat/test-RunModel_Lag.R
+++ b/tests/testthat/test-RunModel_Lag.R
@@ -19,21 +19,6 @@ test_that("'BasinAreas' must have one more element than 'LengthHydro'", {
 
 BasinAreas <- c(BasinInfo$BasinArea, BasinInfo$BasinArea)
 
-test_that("'Qupstream' cannot contain any NA value", {
-  expect_error(
-    InputsModel <- CreateInputsModel(
-      FUN_MOD = RunModel_GR4J,
-      DatesR = BasinObs$DatesR,
-      Precip = BasinObs$P,
-      PotEvap = BasinObs$E,
-      Qupstream = matrix(BasinObs$Qmm, ncol = 1),
-      LengthHydro = 1,
-      BasinAreas = BasinAreas
-    ),
-    regexp = "'Qupstream' cannot contain any NA value"
-  )
-})
-
 # Qupstream = sinusoid synchronised on hydrological year from 0 mm to mean value of Qobs
 Qupstream <- floor((sin((seq_along(BasinObs$Qmm)/365*2*3.14))+1) * mean(BasinObs$Qmm, na.rm = TRUE))
 
@@ -85,7 +70,7 @@ test_that("'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOpt
   )
 })
 
-test_that("'InputsModel$OutputsModel$Qim' should contain no NA'", {
+test_that("'InputsModel$OutputsModel$Qsim' should contain no NA'", {
   InputsModel$OutputsModel <- OutputsGR4JOnly
   InputsModel$OutputsModel$Qsim[10L] <- NA
   expect_error(
@@ -94,6 +79,31 @@ test_that("'InputsModel$OutputsModel$Qim' should contain no NA'", {
   )
 })
 
+test_that("'Qupstream' contain NA values", {
+  expect_warning(
+    InputsModel <- CreateInputsModel(
+      FUN_MOD = RunModel_GR4J,
+      DatesR = BasinObs$DatesR,
+      Precip = BasinObs$P,
+      PotEvap = BasinObs$E,
+      Qupstream = matrix(BasinObs$Qmm, ncol = 1),
+      LengthHydro = 1,
+      BasinAreas = BasinAreas
+    ),
+    regexp = "'Qupstream' contains NA values: model outputs will contain NAs"
+  )
+
+  RunOptions <- suppressWarnings(CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                                  InputsModel = InputsModel,
+                                                  IndPeriod_Run = Ind_Run))
+  InputsModel$OutputsModel <- OutputsGR4JOnly
+
+  expect_warning(
+    RunModel_Lag(InputsModel = InputsModel, RunOptions = RunOptions, Param = 1),
+    regexp = "time steps with NA values"
+  )
+})
+
 test_that("Upstream basin with nil area should return same Qdown as GR4J alone", {
   UpstBasinArea <- InputsModel$BasinAreas[1L]
   InputsModel$BasinAreas[1L] <- 0
diff --git a/vignettes/V05_sd_model.Rmd b/vignettes/V05_sd_model.Rmd
index 430c8e3b62eeb619533a1e1c8480a4318946618e..632819c13cdd8c32d751f23f4933467a030de2d8 100644
--- a/vignettes/V05_sd_model.Rmd
+++ b/vignettes/V05_sd_model.Rmd
@@ -12,7 +12,6 @@ vignette: >
 ```{r, include=FALSE, fig.keep='none', results='hide'}
 library(airGR)
 options(digits = 3)
-library(imputeTS)
 ```
 
 # Introduction
@@ -56,7 +55,9 @@ For the observed flow at the downstream outlet, we generate it with the assumpti
 
 ```{r}
 QObsDown <- (BasinObs$Qmm + c(0, 0, BasinObs$Qmm[1:(length(BasinObs$Qmm)-2)])) / 2
-summary(cbind(QObsUp = BasinObs$Qmm, QObsDown), digits = 3)
+options(digits = 5)
+summary(cbind(QObsUp = BasinObs$Qmm, QObsDown))
+options(digits = 3)
 ```
 
 # Calibration of the upstream subcatchment
@@ -91,19 +92,13 @@ OutputsModelUp <- RunModel_GR4J(InputsModel = InputsModelUp, RunOptions = RunOpt
 
 # Calibration of the downstream subcatchment with upstream flow observations
 
-Observed flow data contain `NA` values and a complete time series is mandatory for running the Lag model. We propose to complete the observed upstream flow with linear interpolation:
-
-```{r}
-QObsUp <- imputeTS::na_interpolation(BasinObs$Qmm)
-```
-
 we need to create the `InputsModel` object completed with upstream information:
 
 ```{r}
 InputsModelDown1 <- CreateInputsModel(
   FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
   Precip = BasinObs$P, PotEvap = BasinObs$E,
-  Qupstream = matrix(QObsUp, ncol = 1), # upstream observed flow
+  Qupstream = matrix(BasinObs$Qmm, ncol = 1), # upstream observed flow
   LengthHydro = 100, # distance between upstream catchment outlet & the downstream one [km]
   BasinAreas = c(180, 180) # upstream and downstream areas [km²]
 )