diff --git a/R/CreateCalibOptions.R b/R/CreateCalibOptions.R
index fc71940be53f63ee7ff9d7b93ccf398e1b434ce3..d944d9dc7954eec2567409b17a2e4cbebf6ec3b2 100644
--- a/R/CreateCalibOptions.R
+++ b/R/CreateCalibOptions.R
@@ -23,17 +23,21 @@ CreateCalibOptions <- function(x, ...) {
 #' @export
 CreateCalibOptions.InputsModel <- function(x,
                                            ...) {
-  if (!exists("FUN_MOD") && !is.null(x$FUN_MOD)) {
-    airGR::CreateCalibOptions(
-      FUN_MOD = x$FUN_MOD,
-      IsSD = !is.null(x$Qupstream) & x$FUN_MOD != "RunModel_Lag",
-      ...
-    )
-  } else {
-    airGR::CreateCalibOptions(
-      ...
-    )
+  dots <- list(...)
+  # Add FUN_MOD in parameters if carried by InputsModel
+  if (!"FUN_MOD" %in% names(dots)) {
+    if(!is.null(x$FUN_MOD)) {
+      dots$FUN_MOD <- x$FUN_MOD
+    } else {
+      stop(" The parameter `FUN_MOD` must be defined")
+    }
   }
+  # Automatically define IsSD for intermediate basin GR models
+  dots$IsSD = !is.null(x$Qupstream) & dots$FUN_MOD != "RunModel_Lag"
+  # Add IsHyst in parameters if carried by InputsModel
+  if (!is.null(x$model$IsHyst)) dots$IsHyst <- x$model$IsHyst
+  # Call airGR function
+  do.call(airGR::CreateCalibOptions, dots)
 }
 
 #' @rdname CreateCalibOptions
diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 50edf1fbd53f26d05818c1278dba3c69eaa90001..e36d26bea6eb412ce1fb6c0e44b683e6ced6c946 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -35,6 +35,8 @@
 #' @param NLayers (optional) named [vector] of [numeric] integer giving the number
 #'        of elevation layers requested [-], required to create CemaNeige module
 #'        inputs, default=5
+#' @param IsHyst [boolean] boolean indicating if the hysteresis version of
+#'        CemaNeige is used. See details of [airGR::CreateRunOptions].
 #' @param ... used for compatibility with S3 methods
 #'
 #' @details Meteorological data are needed for the nodes of the network that
@@ -70,7 +72,8 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
                                      PrecipScale = TRUE,
                                      TempMean = NULL, TempMin = NULL,
                                      TempMax = NULL, ZInputs = NULL,
-                                     HypsoData = NULL, NLayers = 5, ...) {
+                                     HypsoData = NULL, NLayers = 5,
+                                     IsHyst = FALSE, ...) {
 
   # Check and format inputs
   varNames <- c("Precip", "PotEvap", "TempMean", "Qobs", "Qmin",
@@ -326,7 +329,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
 
   # Add the model function
   InputsModel$FUN_MOD <- FUN_MOD
-  featModel <- .GetFeatModel(InputsModel)
+  featModel <- .GetFeatModel(InputsModel, IsHyst)
   InputsModel$isUngauged <- griwrm$model[griwrm$id == id] == "Ungauged"
   InputsModel$gaugedId <- griwrm$donor[griwrm$id == id]
   InputsModel$hasUngaugedNodes <- hasUngaugedNodes(id, griwrm)
@@ -335,7 +338,8 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
       indexParamUngauged = ifelse(inherits(InputsModel, "SD"), 0, 1) +
         seq.int(featModel$NbParam),
       hasX4 = grepl("RunModel_GR[456][HJ]", FUN_MOD),
-      iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4)
+      iX4 = ifelse(inherits(InputsModel, "SD"), 5, 4),
+      IsHyst = featModel$IsHyst
     )
   InputsModel$hasDiversion <- np$Diversion
   InputsModel$isReservoir <- np$Reservoir
@@ -446,7 +450,7 @@ hasUngaugedNodes <- function(id, griwrm) {
 #' function to extract model features partially copied from airGR:::.GetFeatModel
 #' @importFrom utils tail
 #' @noRd
-.GetFeatModel <- function(InputsModel) {
+.GetFeatModel <- function(InputsModel, IsHyst) {
   path <- system.file("modelsFeatures/FeatModelsGR.csv", package = "airGR")
   FeatMod <- read.table(path, header = TRUE, sep = ";", stringsAsFactors = FALSE)
   NameFunMod <- ifelse(test = FeatMod$Pkg %in% "airGR",
@@ -461,6 +465,13 @@ hasUngaugedNodes <- function(id, griwrm) {
   if (FeatMod$IsSD) {
     FeatMod$NbParam <- FeatMod$NbParam + 1
   }
+  FeatMod$IsHyst <- FALSE
+  if (grepl("CemaNeige", FeatMod$NameMod)) {
+    FeatMod$IsHyst <- IsHyst
+    if (IsHyst) {
+      FeatMod$NbParam <- FeatMod$NbParam + 2
+    }
+  }
   return(FeatMod)
 }
 
diff --git a/R/CreateRunOptions.R b/R/CreateRunOptions.R
index e6762391d008c9c850e558160ca951d6ca8f56c7..c74cc68539bd87d676e1d2428b9cc1ddd63c239a 100644
--- a/R/CreateRunOptions.R
+++ b/R/CreateRunOptions.R
@@ -43,6 +43,8 @@ CreateRunOptions.InputsModel <- function(x, ...) {
       stop(" The parameter `FUN_MOD` must be defined")
     }
   }
+  # Add IsHyst in parameters if carried by InputsModel
+  if (!is.null(x$model$IsHyst)) dots$IsHyst <- x$model$IsHyst
 
   # Temporary fix waiting for resolution of HYCAR-Hydro/airgr#167
   if (identical(match.fun(dots$FUN_MOD), RunModel_Lag)) {
diff --git a/man/CreateInputsModel.GRiwrm.Rd b/man/CreateInputsModel.GRiwrm.Rd
index f77b5cebec78cbbdab61ae6f8839cfb17ede6077..f0aa1bcc12b9c9d3fb65b5aa679301e902c92175 100644
--- a/man/CreateInputsModel.GRiwrm.Rd
+++ b/man/CreateInputsModel.GRiwrm.Rd
@@ -18,6 +18,7 @@
   ZInputs = NULL,
   HypsoData = NULL,
   NLayers = 5,
+  IsHyst = FALSE,
   ...
 )
 }
@@ -70,6 +71,9 @@ if not defined a single elevation is used for CemaNeige}
 of elevation layers requested \link{-}, required to create CemaNeige module
 inputs, default=5}
 
+\item{IsHyst}{\link{boolean} boolean indicating if the hysteresis version of
+CemaNeige is used. See details of \link[airGR:CreateRunOptions]{airGR::CreateRunOptions}.}
+
 \item{...}{used for compatibility with S3 methods}
 }
 \value{