diff --git a/R/CreateErrorCrit_GAPX.R b/R/CreateErrorCrit_GAPX.R
index 33b424fa2c0634e6008488c8738b130a1f273619..ac064ae9eb3ce7af433442bd616bb8b2f8cdce45 100644
--- a/R/CreateErrorCrit_GAPX.R
+++ b/R/CreateErrorCrit_GAPX.R
@@ -6,7 +6,7 @@ CreateErrorCrit_GAPX <- function(FUN_TRANSFO) {
       stop("'OutputsModel' must be of class 'OutputsModel'")
     }
 
-    OutputsModel$ParamT <- FUN_TRANSFO(OutputsModel$Param, "RT")
+    OutputsModel$RunOptions$ParamT <- FUN_TRANSFO(OutputsModel$RunOptions$Param, "RT")
 
     EC <- .ErrorCrit(InputsCrit = InputsCrit, Crit = "GAPX", OutputsModel = OutputsModel, warnings = warnings)
 
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 63c99bcbde3e487a1caad606c40dee665b59a068..34ea35b623e42e79aa61d91f299672762d10a943 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -48,10 +48,10 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
 
   if (inherits(QcontribDown, "OutputsModel")) {
     OutputsModel <- QcontribDown
-    if (is.null(OutputsModel$WarmUpQsim)) {
-      OutputsModel$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
+    if (is.null(OutputsModel$RunOptions$WarmUpQsim)) {
+      OutputsModel$RunOptions$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
     }
-    QsimDown <- c(OutputsModel$WarmUpQsim, OutputsModel$Qsim)
+    QsimDown <- c(OutputsModel$RunOptions$WarmUpQsim, OutputsModel$Qsim)
   } else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
     OutputsModel <- list()
     class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
@@ -156,11 +156,11 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
     # message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
   }
   if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
-    OutputsModel$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+    OutputsModel$RunOptions$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))] / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
   }
 
   if ("Param" %in% RunOptions$Outputs_Sim) {
-    OutputsModel$Param <- c(Param, OutputsModel$Param)
+    OutputsModel$RunOptions$Param <- c(Param, OutputsModel$RunOptions$Param)
   }
 
   class(OutputsModel) <- c(class(OutputsModel), "SD")
diff --git a/R/Utils.R b/R/Utils.R
index 19892f3ebc110cf4e5a590361ca795e263fb82aa..505f47f3208666db0c54d652dc54ea37ff2dae31 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -213,6 +213,9 @@
     }
     return(res0)
   })
+  if (!is.null(x$RunOptions)) {
+    res$RunOptions <- x$RunOptions
+  }
   if (!is.null(x$StateEnd)) {
     res$StateEnd <- x$StateEnd
   }
@@ -221,7 +224,7 @@
 }
 
 .IndexOutputsModel <- function(x, i) {
-# '[.OutputsModel' <- function(x, i) {
+  # '[.OutputsModel' <- function(x, i) {
   if (!inherits(x, "OutputsModel")) {
     stop("'x' must be of class 'OutputsModel'")
   }
diff --git a/R/UtilsErrorCrit.R b/R/UtilsErrorCrit.R
index 8c8609f003934bdfa7f798fff18f30cdde989e30..7854146b87391ca8d24650108a7c30408c8c1064 100644
--- a/R/UtilsErrorCrit.R
+++ b/R/UtilsErrorCrit.R
@@ -49,7 +49,7 @@
     Q = OutputsModel$Qsim,
     SCA = rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "Gratio")),
     SWE = rowMeans(sapply(OutputsModel$CemaNeigeLayers[InputsCrit$idLayer], FUN = "[[", "SnowPack")),
-    ParamT = OutputsModel$ParamT
+    ParamT = OutputsModel$RunOptions$ParamT
   )
 
   VarSim[!InputsCrit$BoolCrit] <- NA
diff --git a/R/UtilsRunModel.R b/R/UtilsRunModel.R
index 846f018214f2b4ac21805d727b3974e6be1896ac..17476edeb1060f36c19d33ccb48b1bfea2c14f13 100644
--- a/R/UtilsRunModel.R
+++ b/R/UtilsRunModel.R
@@ -39,18 +39,19 @@
   }
 
   if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
-    OutputsModel$WarmUpQsim <- RESULTS$Outputs[seq_len(length(RunOptions$IndPeriod_WarmUp)),
+    OutputsModel$RunOptions$WarmUpQsim <- RESULTS$Outputs[seq_len(length(RunOptions$IndPeriod_WarmUp)),
                                                which(FortranOutputs == "Qsim")]
-    class(OutputsModel$WarmUpQsim) <- c("WarmUpOutputsModelItem", class(OutputsModel$WarmUpQsim))
+    # class(OutputsModel$RunOptions$WarmUpQsim) <- c("WarmUpOutputsModelItem", class(OutputsModel$RunOptions$WarmUpQsim))
+  }
+
+  if ("Param" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$RunOptions$Param <- Param
   }
 
   if ("StateEnd" %in% RunOptions$Outputs_Sim) {
     OutputsModel$StateEnd <- RESULTS$StateEnd
   }
 
-  if ("Param" %in% RunOptions$Outputs_Sim) {
-    OutputsModel$Param <- Param
-  }
 
   class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
 
diff --git a/vignettes/V05_sd_model.Rmd b/vignettes/V05_sd_model.Rmd
index 6c17770e9c503ee166a9ecc34c723106e65986e5..693dffb38e9f6a38d4abfe19b3bee15b2e5eff23 100644
--- a/vignettes/V05_sd_model.Rmd
+++ b/vignettes/V05_sd_model.Rmd
@@ -122,7 +122,7 @@ For using upstream simulated flows, we should concatenate a vector with the simu
 ```{r}
 Qsim_upstream <- rep(NA, length(BasinObs$DatesR))
 # Simulated flow during warm-up period (365 days before run period)
-Qsim_upstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$WarmUpQsim
+Qsim_upstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$RunOptions$WarmUpQsim
 # Simulated flow during run period
 Qsim_upstream[Ind_Run] <- OutputsModelUp$Qsim