vignettes/V02.1_param_optim.Rmd

 ---
 title: "Plugging in new calibration algorithms in airGR"
-author: "François Bourgin"
+author: "François Bourgin, Guillaume Thirel"
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteEngine{knitr::rmarkdown}
@@ -10,14 +10,18 @@ vignette: >
 
 
 
-```{r, warning=FALSE, include=FALSE, fig.keep='none', results='hide'}
+```{r setup, warning=FALSE, include=FALSE, fig.keep='none', results='hide'}
 library(airGR)
 library(DEoptim)
-library(hydroPSO)
+# library(hydroPSO) # Needs R version >= 3.6 or latticeExtra <= 0.6-28 on R 3.5. Archived on 2023-10-16 as requires archived packages 'hydroTSM' and 'hydroGOF'.
 library(Rmalschains)
+library(caRamel)
+library(ggplot2)
+library(GGally)
 # source("airGR.R")
 set.seed(321)
 load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
 ```
 
 
@@ -37,15 +41,23 @@ Please note that the calibration period is defined in the `CreateRunOptions()` f
 <!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
 <!-- ``` -->
 
-```{r, echo=TRUE, eval=FALSE}
+```{r Calibration_Michel, echo=TRUE, eval=FALSE}
 example("Calibration_Michel")
 ```
 
+In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
+
+```{r RunOptions, results='hide', eval=FALSE}
+RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+                                      IndPeriod_Run = Ind_Run,
+                                      Outputs_Sim = "Qsim")
+```
 
 Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used.
 
 Please note that this vignette is only for illustration purposes and does not provide any guidance about which optimization strategies is recommended for the family of the **GR** models.
 
+
 ## Definition of the necessary function
 
 Parameter estimation can be performed by defining a function that takes a parameter set as input and returns the value of the performance criterion.
@@ -54,8 +66,7 @@ Here we choose to minimize the root mean square error.
 
 The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model.
 
-
-```{r, warning=FALSE, results='hide', eval=FALSE}
+```{r OptimGR4J, warning=FALSE, results='hide', eval=FALSE}
 OptimGR4J <- function(ParamOptim) {
   ## Transformation of the parameter set to real space
   RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
@@ -74,29 +85,38 @@ OptimGR4J <- function(ParamOptim) {
 
 
 In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space:
-```{r, warning=FALSE, results='hide', eval=FALSE}
+
+```{r boundsGR4J, warning=FALSE, results='hide', eval=FALSE}
 lowerGR4J <- rep(-9.99, times = 4)
 upperGR4J <- rep(+9.99, times = 4)
 ```
 
+
+
 # Local optimization
 
 We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space:
-```{r, warning=FALSE, results='hide', eval=FALSE}
-optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7), 
+
+```{r local1, warning=FALSE, results='hide', eval=FALSE}
+startGR4J <- c(4.1, 3.9, -0.9, -8.7)
+optPORT <- stats::nlminb(start = startGR4J,
                          objective = OptimGR4J,
                          lower = lowerGR4J, upper = upperGR4J,
                          control = list(trace = 1))
 ```
+
 The RMSE value reaches a local minimum value after 35 iterations.
 
 We can also try a multi-start approach to test the consistency of the local optimization.
 Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function).
 For each starting point, a local optimization is performed.
-```{r, warning=FALSE, results='hide', eval=FALSE}
-startGR4J <- expand.grid(data.frame(CalibOptions$StartParamDistrib))
+
+```{r local2, warning=FALSE, results='hide', eval=FALSE}
+startGR4JDistrib <- TransfoParam_GR4J(ParamIn = CalibOptions$StartParamDistrib,
+                                      Direction = "RT")
+startGR4J <- expand.grid(data.frame(startGR4JDistrib))
 optPORT_ <- function(x) {
-  opt <- stats::nlminb(start = x, 
+  opt <- stats::nlminb(start = x,
                        objective = OptimGR4J,
                        lower = lowerGR4J, upper = upperGR4J,
                        control = list(trace = 1))
@@ -105,20 +125,23 @@ listOptPORT <- apply(startGR4J, MARGIN = 1, FUN = optPORT_)
 ```
 
 We can then extract the best parameter sets and the value of the performance criteria:
-```{r, warning=FALSE, results='hide', eval=FALSE}
+
+```{r local3, warning=FALSE, results='hide', eval=FALSE}
 parPORT <- t(sapply(listOptPORT, function(x) x$par))
 objPORT <- sapply(listOptPORT, function(x) x$objective)
-parPORT  <- data.frame(parPORT, RMSE = objPORT)
+resPORT <- data.frame(parPORT, RMSE = objPORT)
 ```
 
 As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets.
-```{r, warning=FALSE}
-summary(parPORT)
+
+```{r local4, warning=FALSE}
+summary(resPORT)
 ```
 
 The existence of several local minima illustrates the importance of defining an appropriate starting point or of using a multi-start strategy or a global optimization strategy.
 
 
+
 # Global optimization
 
 Global optimization is most often used when facing a complex response surface, with multiple local mimina.
@@ -128,8 +151,10 @@ Here we use the following R implementation of some popular strategies:
 * [hydroPSO: particle swarm](https://cran.r-project.org/package=hydroPSO)
 * [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains)
 
+
 ## Differential Evolution
-```{r, warning=FALSE, results='hide', eval=FALSE}
+
+```{r optDE, warning=FALSE, results='hide', eval=FALSE}
 optDE <- DEoptim::DEoptim(fn = OptimGR4J,
                           lower = lowerGR4J, upper = upperGR4J,
                           control = DEoptim::DEoptim.control(NP = 40, trace = 10))
@@ -137,38 +162,146 @@ optDE <- DEoptim::DEoptim(fn = OptimGR4J,
 
 
 ## Particle Swarm
-```{r, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
+
+```{r hydroPSO1, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
+# to install the package temporary removed from CRAN
+# Rtools needed (windows : https://cran.r-project.org/bin/windows/Rtools/)
+# install.packages("https://cran.r-project.org/src/contrib/Archive/hydroPSO/hydroPSO_0.5-1.tar.gz",
+#                  repos = NULL, type = "source", dependencies = TRUE)
+```
+
+```{r hydroPSO2, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
 optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
                              lower = lowerGR4J, upper = upperGR4J,
                              control = list(write2disk = FALSE, verbose = FALSE))
 ```
 
+
 ## MA-LS-Chains
-```{r, warning=FALSE, results='hide', eval=FALSE}
+
+```{r optMALS, warning=FALSE, results='hide', eval=FALSE}
 optMALS <- Rmalschains::malschains(fn = OptimGR4J,
-                                   lower = lowerGR4J, upper = upperGR4J, 
+                                   lower = lowerGR4J, upper = upperGR4J,
                                    maxEvals = 2000)
 ```
 
+
+
 # Results
 
 As it can be seen in the table below, the four additional optimization strategies tested lead to very close optima.
 
-```{r, warning=FALSE, echo=FALSE, eval=FALSE}
-parGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"), 
+```{r resGLOB1, warning=FALSE, echo=FALSE, eval=FALSE}
+resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
                       round(rbind(
-                        OutputsCalib$ParamFinalR                          ,
+                        OutputsCalib$ParamFinalR,
                         airGR::TransfoParam_GR4J(ParamIn = optPORT$par                    , Direction = "TR"),
                         airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
                         airGR::TransfoParam_GR4J(ParamIn = as.numeric(optPSO$par)         , Direction = "TR"),
                         airGR::TransfoParam_GR4J(ParamIn = optMALS$sol                    , Direction = "TR")),
                         digits = 3))
 ```
-```{r, warning=FALSE, echo=FALSE}
-parGLOB
+```{r resGLOB2, warning=FALSE, echo=FALSE}
+resGLOB
 ```
 
 <!-- This is an expected result because the response surface for quadratic performance criteria of the **GR4J** model is generally sufficiently well defined in the transformed parameter space to allow using a local optimization strategy instead of a more time consuming global optimization strategy. -->
 
 
+# Multiobjective optimization
+
+Multiobjective optimization is used to explore possible trade-offs between different performances criteria.
+Here we use the following R implementation of an efficient strategy:
+
+* [caRamel: Automatic Calibration by Evolutionary Multi Objective Algorithm](https://cran.r-project.org/package=caRamel)
+
+Motivated by using the rainfall-runoff model for low flow simulation, we explore the trade-offs between the KGE values obtained without any data transformation and with the inverse transformation.
 
+First, the `OptimGR4J()` function previously used is modified to return two values.
+
+```{r, warning=FALSE, results='hide', eval=FALSE}
+InputsCrit_inv <- InputsCrit
+InputsCrit_inv$transfo <- "inv"
+
+MOptimGR4J <- function(i) {
+  if (algo == "caRamel") {
+    ParamOptim <- x[i, ]
+  }
+  ## Transformation of the parameter set to real space
+  RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
+                                            Direction = "TR")
+  ## Simulation given a parameter set
+  OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
+                                       RunOptions = RunOptions,
+                                       Param = RawParamOptim)
+  ## Computation of the value of the performance criteria
+  OutputsCrit1 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit,
+                                       OutputsModel = OutputsModel,
+                                       verbose = FALSE)
+  ## Computation of the value of the performance criteria
+  OutputsCrit2 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit_inv,
+                                       OutputsModel = OutputsModel,
+                                       verbose = FALSE)
+  return(c(OutputsCrit1$CritValue, OutputsCrit2$CritValue))
+}
+```
+
+
+## caRamel
+
+caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm.
+
+```{r, warning=FALSE, results='hide', eval=FALSE}
+algo <- "caRamel"
+optMO <- caRamel::caRamel(nobj = 2,
+                          nvar = 4,
+                          minmax = rep(TRUE, 2),
+                          bounds = matrix(c(lowerGR4J, upperGR4J), ncol = 2),
+                          func = MOptimGR4J,
+                          popsize = 100,
+                          archsize = 100,
+                          maxrun = 15000,
+                          prec = rep(1.e-3, 2),
+                          carallel = FALSE,
+                          graph = FALSE)
+```
+
+The algorithm returns parameter sets that describe the pareto front, illustrating the trade-off between overall good performance and good performance for low flow.
+
+```{r, fig.width=6, fig.height=6, warning=FALSE}
+ggplot() +
+  geom_point(aes(optMO$objectives[, 1], optMO$objectives[, 2])) +
+  coord_equal(xlim = c(0.4, 0.9), ylim = c(0.4, 0.9)) +
+  xlab("KGE") + ylab("KGE [1/Q]") +
+  theme_bw()
+```
+
+The parameter sets can be viewed in the parameter space, illustrating different populations.
+
+```{r fig.height=6, fig.width=6, message=FALSE, warning=FALSE}
+param_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
+  airGR::TransfoParam_GR4J(x, Direction = "TR")
+  })
+GGally::ggpairs(data.frame(t(param_optMO)), diag = NULL) + theme_bw()
+```
+
+```{r fig.height=6, fig.width=12, message=FALSE, warning=FALSE}
+RunOptions$Outputs_Sim <- "Qsim"
+run_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
+  airGR::RunModel_GR4J(InputsModel = InputsModel,
+                       RunOptions = RunOptions,
+                       Param = x)
+  }$Qsim)
+run_optMO <- data.frame(run_optMO)
+
+ggplot() +
+  geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
+                y = run_optMO$X1)) +
+  geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
+                y = run_optMO$X54),
+            colour = "darkred") +
+  scale_x_datetime(limits = c(as.POSIXct("1998-01-01"), NA)) +
+  ylab("Discharge [mm/d]") + xlab("Date") +
+  scale_y_sqrt() +
+  theme_bw()
+```

vignettes/V02.2_param_mcmc.Rmd

@@ -15,8 +15,6 @@ library(airGR)
 library(coda)
 library(FME)
 library(ggmcmc)
-library(dplyr)
-# source("airGR.R")
 set.seed(123)
 load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))
 ```
@@ -38,6 +36,14 @@ We use the **GR4J** model and we assume that the R global environment contains d
 example("Calibration_Michel")
 ```
 
+In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
+
+```{r, results='hide', eval=FALSE}
+RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
+                                      IndPeriod_Run = Ind_Run,
+                                      Outputs_Sim = "Qsim")
+```
+
 Please refer to the [2.1 Plugging in new calibration](V02.1_param_optim.html) vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**.
 
 
@@ -48,8 +54,9 @@ Please note that this vignette is only for illustration purposes and does not pr
 We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modMCMC()` function of the [FME](https://cran.r-project.org/package=FME) package.
 First, we need to define a function that returns twice the opposite of the log-likelihood for a given parameter set.
 
-Nota: in the `LogLikeGR4J()` function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation  of the two following lines. 
-```{r, echo=TRUE, eval=FALSE}
+Nota: in the `LogLikeGR4J()` function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation  of the two following lines.
+
+```{r, echo=TRUE, eval=FALSE, purl=FALSE}
 Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2)
 LogLike <- -2 * log(Likelihood)
 ```
@@ -62,13 +69,13 @@ Note that we do not use here any discharge transformation, although applying Box
 LogLikeGR4J <- function(ParamOptim) {
   ## Transformation to real space
   RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
-                                              Direction = "TR")
+                                            Direction = "TR")
   ## Simulation given a parameter set
   OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
                                        RunOptions = RunOptions,
                                        Param = RawParamOptim)
   ## Computation of the log-likelihood: N * log(SS)
-  ObsY <- InputsCrit$Qobs
+  ObsY <- InputsCrit$Obs
   ModY <- OutputsModel$Qsim
   LogLike <- sum(!is.na(ObsY)) * log(sum((ObsY - ModY)^2, na.rm = TRUE))
 }
@@ -79,7 +86,7 @@ LogLikeGR4J <- function(ParamOptim) {
 ## Estimation of the best-fit parameters as a starting point
 We start by using the PORT optimization routine to estimate the best-fit parameters.
 ```{r, results='hide', eval=FALSE}
-optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7), 
+optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7),
                          objective = LogLikeGR4J,
                          lower = rep(-9.9, times = 4), upper = rep(9.9, times = 4),
                          control = list(trace = 1))
@@ -96,8 +103,10 @@ Nota: in this example, there are relatively few iterations (2000), in order to l
 With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied.
 
 ```{r, results='hide', eval=FALSE}
-iniParPORT <- data.frame(Chain1 = iniParPORT, Chain2 = iniParPORT, Chain3 = iniParPORT,
-                           row.names = paste0("X", 1:4))
+iniParPORT <- data.frame(Chain1 = iniParPORT,
+                         Chain2 = iniParPORT,
+                         Chain3 = iniParPORT,
+                         row.names = paste0("X", 1:4))
 iniParPORT <- sweep(iniParPORT, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*")
 iniParPORT[iniParPORT < -9.9] <- -9.9
 iniParPORT[iniParPORT > +9.9] <- +9.9
@@ -111,8 +120,8 @@ mcmcDRAM <- apply(iniParPORT, MARGIN = 2, FUN = function(iIniParPORT) {
                jump         = 0.01,
                outputlength = 2000,
                burninlength = 0,
-               updatecov    = 100, ## Adaptative Metropolis
-               ntrydr       = 2)   ## Delayed Rejection
+               updatecov    = 100, ## adaptative Metropolis (AM)
+               ntrydr       = 2)   ## delayed rejection (RD)
 })
 ```
 
@@ -143,13 +152,13 @@ In addition, graphical tools can be used, with for example the [ggmcmc](https://
 First, the evolution of the Markov chains can be seen with a traceplot:
 
 ```{r, fig.width=6, fig.height=9, warning=FALSE}
-parDRAM <- ggmcmc::ggs(multDRAM) ## to convert objet for using by all ggs_* graphical functions
+parDRAM <- ggmcmc::ggs(multDRAM) ## to convert object for using by all ggs_* graphical functions
 ggmcmc::ggs_traceplot(parDRAM)
 ```
 
 The posterior density for each parameter can then be visualised:
 ```{r, fig.width=6, fig.height=9, warning=FALSE}
-burnParDRAM <- dplyr::filter(parDRAM, Iteration > 1000) # to keep only the second half of the series
+burnParDRAM <- parDRAM[parDRAM$Iteration > 1000, ] # to keep only the second half of the series
 ggmcmc::ggs_density(burnParDRAM)
 ```
 

vignettes/V03_param_sets_GR4J.Rmd

 ---
 title: "Generalist parameter sets for the GR4J model"
+author: "Olivier Delaigue, Guillaume Thirel"
 bibliography: V00_airgr_ref.bib
 output: rmarkdown::html_vignette
 vignette: >
@@ -36,7 +37,6 @@ In this vignette, we show a method using a known parameter set that can be used
 We load an example data set from the package and the GR4J parameter sets.
 
 ```{r, warning=FALSE}
-
 ## loading catchment data
 data(L0123001)
 
@@ -70,7 +70,6 @@ The calibration period has been defined from **1990-01-01** to **1990-02-28**, a
 As a consequence, in this example the calibration period is very short, less than 6 months.
 
 ```{r, warning=FALSE, include=TRUE}
-
 ## preparation of the InputsModel object
 InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR, 
                                  Precip = BasinObs$P, PotEvap = BasinObs$E)
@@ -83,11 +82,11 @@ Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/
 
 ## preparation of the RunOptions object for the calibration period
 RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
-                               InputsModel = InputsModel, IndPeriod_Run = Ind_Cal)
+                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Cal)
 
 ## efficiency criterion: Nash-Sutcliffe Efficiency
 InputsCrit_Cal  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, 
-                                    RunOptions = RunOptions_Cal, Qobs = BasinObs$Qmm[Ind_Cal])
+                                    RunOptions = RunOptions_Cal, Obs = BasinObs$Qmm[Ind_Cal])
 
 
 ## ---- validation step
@@ -102,8 +101,7 @@ RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
 
 ## efficiency criterion (Nash-Sutcliffe Efficiency) on the validation period
 InputsCrit_Val  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, 
-                                RunOptions = RunOptions_Val, Qobs = BasinObs$Qmm[Ind_Val])
-
+                                    RunOptions = RunOptions_Val, Obs = BasinObs$Qmm[Ind_Val])
 ```
 
 # Calibration of the GR4J model with the generalist parameter sets
@@ -134,7 +132,8 @@ Param_Best <- unlist(Param_Sets_GR4J[which.max(OutputsCrit_Loop), ])
 Param_Best
 
 ## validation
-OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val, Param = Param_Best)
+OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val, 
+                                  Param = Param_Best)
 OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
 ```
 
@@ -142,7 +141,7 @@ Now we can compute the Nash-Sutcliffe Efficiency criterion on the validation per
 
 
 
-# Calibration of the GR4J model with the built-in `Calibration_Michel()` function
+# Calibration of the GR4J model with the built-in Calibration_Michel function
 
 As seen above, the Michel's calibration algorithm is based on a local search procedure.
 
@@ -156,22 +155,22 @@ By default, the predefined grid used by the `Calibration_Michel()` function cont
 
 ```{r, warning=FALSE, message=FALSE}
 CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
-
 ```
 
 ```{r, warning=FALSE, message=FALSE, include=FALSE}
 ## calibration
 OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
                             InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
-                            FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE, 
+                            FUN_MOD = RunModel_GR4J,
                             FUN_CALIB = Calibration_Michel)
 OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal, 
-                             Param = OutputsCalib$ParamFinalR, FUN = RunModel_GR4J)
+                             Param = OutputsCalib$ParamFinalR, FUN_MOD = RunModel_GR4J)
 OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
 
 
 ## validation 
-OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val, Param = OutputsCalib$ParamFinalR)
+OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val,
+                                  Param = OutputsCalib$ParamFinalR)
 OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
 ```
 
@@ -199,10 +198,10 @@ CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibrat
 ## calibration
 OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
                             InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
-                            FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE, 
+                            FUN_MOD = RunModel_GR4J,
                             FUN_CALIB = Calibration_Michel)
 OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal, 
-                         Param = OutputsCalib$ParamFinalR, FUN = RunModel_GR4J)
+                             Param = OutputsCalib$ParamFinalR, FUN_MOD = RunModel_GR4J)
 OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
 
 

vignettes/V04_cemaneige_hysteresis.Rmd

+---
+title: "Using satellite snow cover area data for calibrating and improving CemaNeige"
+author: "Guillaume Thirel, Olivier Delaigue"
+bibliography: V00_airgr_ref.bib
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteIndexEntry{Using CemaNeige with Hysteresis}
+  %\VignetteEncoding{UTF-8}
+---
+
+
+```{r, warning=FALSE, include=FALSE}
+library(airGR)
+load(system.file("vignettesData/vignetteCNHysteresis.rda", package = "airGR"))
+```
+
+
+# Introduction
+
+## Scope
+
+Rainfall-runoff models that include a snow accumulation and melt module are still often calibrated using only discharge observations. This can result in poor snow modeling as the swnow module parameters can distorted in order to allow skilful discharge simulation.
+  
+After the work of @riboust_revisiting_2019, we propose now in **airGR** an improved version of the degree-day CemaNeige snow and accumulation module. This new version is based on a more accurate representation of the relationship that exists at the basin scale between the Snow Water Equivalent (SWE) and the Snow Cover Area (SCA). To do so, a linear SWE-SCA hysteresis, which represents the fact that snow accumulation is rather homogeneous on the basin and snow melt is more heterogeneous, was implemented.
+  
+This new CemaNeige version has two more parameters to calibrate. It also presents the advantage of allowing using satellite snow data to constrain the calibration in addition to discharge. 
+@riboust_revisiting_2019 show that while the simulated discharge is not significantly improved, the snow simulation is much improved. In addition, they show that the model is more robust (i.e. transferable in time) in terms of discharge, which has many implications for climate change impact studies.
+  
+The configuration that was identified as optimal by @riboust_revisiting_2019 includes a CemaNeige module run on 5 elevation bands and an objective function determine by a composite function of KGE' calculated on discharge (75 % weight) and KGE' calculated on each elevation band (5 % for each).
+  
+In this page, we show how to use and calibrate this new CemaNeige version. 
+  
+  
+## Data preparation
+  
+We load an example data set from the package. Please note that this data set includes MODIS data that was pre-calculated for 5 elevation bands and for which days with few data (more than 40 % cloud coverage) were assigned as missing values. 
+  
+  
+## loading catchment data
+```{r, warning=FALSE, eval=FALSE}
+data(X0310010)
+summary(BasinObs)
+```
+
+
+## Object model preparation
+
+We assume that the R global environment contains data and functions from the [Get Started](V01_get_started.html) vignette.
+
+The calibration period has been defined from 2000-09-01 to 2005-08-31, and the validation period from 2005-09-01 to 2010-07-31. CemaNeige will be used in coupling with GR4J in this vignette.
+
+
+```{r, warning=FALSE, eval=FALSE}
+## preparation of the InputsModel object
+InputsModel <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR4J,
+                                 DatesR = BasinObs$DatesR, Precip = BasinObs$P,
+                                 PotEvap = BasinObs$E, TempMean = BasinObs$T,
+                                 ZInputs = median(BasinInfo$HypsoData),
+                                 HypsoData = BasinInfo$HypsoData, NLayers = 5)
+
+## ---- calibration step
+
+## calibration period selection
+Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2000-09-01"), 
+               which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2005-08-31"))
+
+
+## ---- validation step
+
+## validation period selection
+Ind_Val <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2005-09-01"), 
+               which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2010-07-31"))
+```
+
+
+
+# Calibration and evaluation of the new CemaNeige module
+
+In order to use the Linear Hysteresis, a new argument (`IsHyst`) is added in the `CreateRunOptions()` and `CreateCalibOptions()` functions and has to be set to `TRUE`. 
+
+```{r, warning=FALSE, eval=FALSE}
+## preparation of the RunOptions object for the calibration period
+RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
+                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Cal,
+                                   IsHyst = TRUE)
+
+## preparation of the RunOptions object for the validation period
+RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
+                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Val,
+                                   IsHyst = TRUE)
+
+## preparation of the CalibOptions object
+CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
+                                   FUN_CALIB = Calibration_Michel,
+                                   IsHyst = TRUE)
+```
+
+In order to calibrate and assess the model performance, we will follow the recommendations of @riboust_revisiting_2019. This is now possible in **airGR** with the added functionality that permits calculating composite criteria by combining different metrics. 
+
+
+```{r, warning=FALSE, eval=FALSE}
+## efficiency criterion: 75 % KGE'(Q) + 5 % KGE'(SCA) on each of the 5 layers
+InputsCrit_Cal  <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
+                                    InputsModel = InputsModel, RunOptions = RunOptions_Cal,
+                                    Obs = list(BasinObs$Qmm[Ind_Cal],
+                                               BasinObs$SCA1[Ind_Cal],
+                                               BasinObs$SCA2[Ind_Cal],
+                                               BasinObs$SCA3[Ind_Cal],
+                                               BasinObs$SCA4[Ind_Cal],
+                                               BasinObs$SCA5[Ind_Cal]),
+                                    VarObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
+                                    Weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))
+
+InputsCrit_Val  <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
+                                    InputsModel = InputsModel, RunOptions = RunOptions_Val,
+                                    Obs = list(BasinObs$Qmm[Ind_Val],
+                                               BasinObs$SCA1[Ind_Val],
+                                               BasinObs$SCA2[Ind_Val],
+                                               BasinObs$SCA3[Ind_Val],
+                                               BasinObs$SCA4[Ind_Val],
+                                               BasinObs$SCA5[Ind_Val]),
+                                    VarObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
+                                    Weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))
+``` 
+
+We can now calibrate the model.
+
+```{r, warning=FALSE, eval=FALSE}
+## calibration
+OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
+                            InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
+                            FUN_MOD = RunModel_CemaNeigeGR4J,
+                            FUN_CALIB = Calibration_Michel)
+```
+
+
+Now we can run it on the calibration period and assess it. 
+
+```{r, warning=FALSE, message=FALSE, eval=FALSE}
+## run on the calibration period
+OutputsModel_Cal <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
+                                           RunOptions = RunOptions_Cal, 
+                                           Param = OutputsCalib$ParamFinalR)
+
+## evaluation 
+OutputsCrit_Cal <- ErrorCrit(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
+```
+
+
+Find below the performance of the model over the calibration period. 
+
+```{r, warning=FALSE}
+str(OutputsCrit_Cal, max.level = 2)
+```
+
+Now we can run the model on the validation period and assess it. 
+
+```{r, warning=FALSE, message=FALSE, eval=FALSE}
+## run on the validation period
+OutputsModel_Val <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
+                                           RunOptions = RunOptions_Val, 
+                                           Param = OutputsCalib$ParamFinalR)
+
+## evaluation 
+OutputsCrit_Val <- ErrorCrit(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
+```
+
+Find below the performance of the model over the validation period. 
+
+```{r, warning=FALSE}
+str(OutputsCrit_Val, max.level = 2)
+```
+
+
+
+# Comparison with the performance of the initial CemaNeige version 
+
+Here we use the same InputsModel object and calibration and validation periods. However, we have to redefine the way we run the model (`RunOptions` argument), calibrate and assess it (`InputsCrit` argument). The objective function is only based on KGE'(Q). Note how we set the `IsHyst` argument to `FALSE` in the `CreateRunOptions()`  and the  `CreateCalibOptions()` functions.
+
+```{r, warning=FALSE, eval=FALSE}
+## preparation of RunOptions object
+RunOptions_Cal_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
+                                          InputsModel = InputsModel,
+                                          IndPeriod_Run = Ind_Cal,
+                                          IsHyst = FALSE)
+
+RunOptions_Val_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
+                                          InputsModel = InputsModel, 
+                                          IndPeriod_Run = Ind_Val,
+                                          IsHyst = FALSE)
+
+InputsCrit_Cal_NoHyst <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
+                                          InputsModel = InputsModel,
+                                          RunOptions = RunOptions_Cal_NoHyst,
+                                          Obs = BasinObs$Qmm[Ind_Cal], VarObs = "Q")
+
+## preparation of CalibOptions object
+CalibOptions_NoHyst <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
+                                          FUN_CALIB = Calibration_Michel,
+                                          IsHyst = FALSE)
+```
+
+We can now calibrate the model.
+
+```{r, warning=FALSE, eval=FALSE}
+## calibration
+OutputsCalib_NoHyst <- Calibration(InputsModel = InputsModel,
+                                   InputsCrit = InputsCrit_Cal_NoHyst,
+                                   RunOptions = RunOptions_Cal_NoHyst,
+                                   CalibOptions = CalibOptions_NoHyst,
+                                   FUN_MOD = RunModel_CemaNeigeGR4J,
+                                   FUN_CALIB = Calibration_Michel)
+```
+
+And run it over the calibration and validation periods. 
+
+```{r, warning=FALSE, eval=FALSE}
+OutputsModel_Cal_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
+                                                  RunOptions = RunOptions_Cal_NoHyst,  
+                                                  Param = OutputsCalib_NoHyst$ParamFinalR)
+
+OutputsModel_Val_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
+                                                  RunOptions = RunOptions_Val_NoHyst, 
+                                                  Param = OutputsCalib_NoHyst$ParamFinalR)
+```
+
+In order to assess the model performance over the two periods, we will use the InputsCrit objects prepared before, which allow assessing also the performance in terms of snow simulation.
+
+
+```{r, warning=FALSE, message=FALSE, eval=FALSE}
+OutputsCrit_Cal_NoHyst <- ErrorCrit(InputsCrit = InputsCrit_Cal,
+                                    OutputsModel = OutputsModel_Cal_NoHyst)
+
+OutputsCrit_Val_NoHyst <- ErrorCrit(InputsCrit = InputsCrit_Val,
+                                    OutputsModel = OutputsModel_Val_NoHyst)
+```
+
+We can check the performance over the calibration and the validation period.
+
+```{r, warning=FALSE}
+str(OutputsCrit_Cal_NoHyst, max.level = 2)
+str(OutputsCrit_Val_NoHyst, max.level = 2)
+```
+
+We can see above that the performance of the initial model is slightly better than the new one over the calibration period in terms of discharge, but also that the new version calibrated using SCA provides better performance in terms of snow. 
+However, over the validation period, we see that the discharge simulated by the new version brings better performance (in addition to improved SCA also). This shows the interests of the combined use of a linear hysteresis and of SCA data for calibration in CemaNeige. 
+
+
+
+# References

vignettes/V05_sd_model.Rmd

+---
+title: "Simulated vs observed upstream flows in calibration of semi-distributed GR4J model"
+author: "David Dorchies"
+bibliography: V00_airgr_ref.bib
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteIndexEntry{Simulated vs observed upstream flows in calibration of semi-distributed GR4J model}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include=FALSE, fig.keep='none', results='hide'}
+library(airGR)
+options(digits = 3)
+```
+
+# Introduction
+
+## Scope
+
+The **airGR** package implements semi-distributed model capabilities using a lag model between subcatchments. It allows to chain together several lumped models as well as integrating anthropogenic influence such as reservoirs or withdrawals.
+
+Here we explain how to implement the semi-distribution with **airGR**. For everyday use, however, it is easier to use the [**airGRiwrm**](https://cran.r-project.org/package=airGRiwrm) package.
+
+`RunModel_Lag` documentation gives an example of simulating the influence of a reservoir in a lumped model. Try `example(RunModel_Lag)` to get it.
+
+In this vignette, we show how to calibrate 2 sub-catchments in series with a semi-distributed model consisting of 2 GR4J models.
+For doing this we compare 3 strategies for calibrating the downstream subcatchment:
+
+- using upstream observed flows
+- using upstream simulated flows
+- using upstream simulated flows and parameter regularisation [@delavenne_regularization_2019]
+
+We finally compare these calibrations with a theoretical set of parameters.
+This comparison is based on the Kling-Gupta Efficiency computed on the root-squared discharges as performance criterion.
+
+## Model description
+
+```{r, warning=FALSE, include=FALSE}
+library(airGR)
+options(digits = 3)
+```
+
+We use an example data set from the package that unfortunately contains data for only one catchment.
+
+```{r, warning=FALSE}
+## loading catchment data
+data(L0123001)
+```
+
+Let's imagine that this catchment of 360 km² is divided into 2 subcatchments:
+
+- An upstream subcatchment of 180 km²
+- 100 km downstream another subcatchment of 180 km²
+
+We consider that meteorological data are homogeneous on the whole catchment, so we use the same pluviometry `BasinObs$P` and the same evapotranspiration `BasinObs$E` for the 2 subcatchments.
+
+For the observed flow at the downstream outlet, we generate it with the assumption that the upstream flow arrives at downstream with a constant delay of 2 days.
+
+```{r}
+QObsDown <- (BasinObs$Qmm + c(0, 0, BasinObs$Qmm[1:(length(BasinObs$Qmm)-2)])) / 2
+options(digits = 5)
+summary(cbind(QObsUp = BasinObs$Qmm, QObsDown))
+options(digits = 3)
+```
+
+With a delay of 2 days between the 2 gauging stations, the theoretical Velocity parameter should be equal to:
+
+```{r}
+Velocity <- 100 * 1e3 / (2 * 86400)
+paste("Velocity: ", format(Velocity), "m/s")
+```
+
+# Calibration of the upstream subcatchment
+
+The operations are exactly the same as the ones for a GR4J lumped model. So we do exactly the same operations as in the [Get Started](V01_get_started.html) vignette.
+
+```{r}
+InputsModelUp <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+                                   Precip = BasinObs$P, PotEvap = BasinObs$E)
+Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
+               which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
+RunOptionsUp <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                 InputsModel = InputsModelUp,
+                                 IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
+                                 IniStates = NULL, IniResLevels = NULL)
+# Error criterion is KGE computed on the root-squared discharges
+InputsCritUp <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelUp,
+                                 RunOptions = RunOptionsUp,
+                                 VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run],
+                                 transfo = "sqrt")
+CalibOptionsUp <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
+OutputsCalibUp <- Calibration_Michel(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
+                                     InputsCrit = InputsCritUp, CalibOptions = CalibOptionsUp,
+                                     FUN_MOD = RunModel_GR4J)
+```
+
+And see the result of the simulation:
+
+```{r}
+OutputsModelUp <- RunModel_GR4J(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
+                                Param = OutputsCalibUp$ParamFinalR)
+```
+
+
+# Calibration of the downstream subcatchment
+
+## Creation of the InputsModel objects
+
+We need to create `InputsModel` objects completed with upstream information with upstream observed flow for the calibration of first case and upstream simulated flows for the other cases:
+
+```{r}
+InputsModelDown1 <- CreateInputsModel(
+  FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+  Precip = BasinObs$P, PotEvap = BasinObs$E,
+  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²]
+)
+```
+
+For using upstream simulated flows, we should concatenate a vector with the simulated flows for the entire period of simulation (warm-up + run):
+
+```{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$RunOptions$WarmUpQsim
+# Simulated flow during run period
+Qsim_upstream[Ind_Run] <- OutputsModelUp$Qsim
+
+InputsModelDown2 <- CreateInputsModel(
+  FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
+  Precip = BasinObs$P, PotEvap = BasinObs$E,
+  Qupstream = matrix(Qsim_upstream, 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²]
+)
+```
+
+
+## Calibration with upstream flow observations
+
+We calibrate the combination of Lag model for upstream flow transfer and GR4J model for the runoff of the downstream subcatchment:
+
+```{r}
+RunOptionsDown <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+                                   InputsModel = InputsModelDown1,
+                                   IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
+                                   IniStates = NULL, IniResLevels = NULL)
+InputsCritDown <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelDown1,
+                                   RunOptions = RunOptionsDown,
+                                   VarObs = "Q", Obs = QObsDown[Ind_Run],
+                                   transfo = "sqrt")
+CalibOptionsDown <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
+                                       FUN_CALIB = Calibration_Michel,
+                                       IsSD = TRUE) # specify that it's a SD model
+OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
+                                        RunOptions = RunOptionsDown,
+                                        InputsCrit = InputsCritDown,
+                                        CalibOptions = CalibOptionsDown,
+                                        FUN_MOD = RunModel_GR4J)
+```
+
+`RunModel` is run in order to automatically combine GR4J and Lag models.
+
+```{r}
+OutputsModelDown1 <- RunModel(InputsModel = InputsModelDown2,
+                              RunOptions = RunOptionsDown,
+                              Param = OutputsCalibDown1$ParamFinalR,
+                              FUN_MOD = RunModel_GR4J)
+```
+
+Performance of the model validation is then:
+
+```{r}
+KGE_down1 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown1)
+```
+
+
+## Calibration with upstream simulated flow
+
+We calibrate the model with the `InputsModel` object previously created for substituting the observed upstream flow with the simulated one:
+
+```{r}
+OutputsCalibDown2 <- Calibration_Michel(InputsModel = InputsModelDown2,
+                                        RunOptions = RunOptionsDown,
+                                        InputsCrit = InputsCritDown,
+                                        CalibOptions = CalibOptionsDown,
+                                        FUN_MOD = RunModel_GR4J)
+ParamDown2 <- OutputsCalibDown2$ParamFinalR
+```
+
+## Calibration with upstream simulated flow and parameter regularisation
+
+The regularisation follow the method proposed by @delavenne_regularization_2019.
+
+As a priori parameter set, we use the calibrated parameter set of the upstream catchment and the theoretical velocity:
+
+```{r}
+ParamDownTheo <- c(Velocity, OutputsCalibUp$ParamFinalR)
+```
+
+The Lavenne criterion is initialised with the a priori parameter set and the value of the KGE of the upstream basin.
+
+```{r}
+IC_Lavenne <- CreateInputsCrit_Lavenne(InputsModel = InputsModelDown2,
+                                       RunOptions = RunOptionsDown,
+                                       Obs = QObsDown[Ind_Run],
+                                       AprParamR = ParamDownTheo,
+                                       AprCrit = OutputsCalibUp$CritFinal)
+```
+
+The Lavenne criterion is used instead of the KGE for calibration with regularisation
+
+```{r}
+OutputsCalibDown3 <- Calibration_Michel(InputsModel = InputsModelDown2,
+                                        RunOptions = RunOptionsDown,
+                                        InputsCrit = IC_Lavenne,
+                                        CalibOptions = CalibOptionsDown,
+                                        FUN_MOD = RunModel_GR4J)
+```
+
+The KGE is then calculated for performance comparisons:
+
+```{r}
+OutputsModelDown3 <- RunModel(InputsModel = InputsModelDown2,
+                              RunOptions = RunOptionsDown,
+                              Param = OutputsCalibDown3$ParamFinalR,
+                              FUN_MOD = RunModel_GR4J)
+KGE_down3 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown3)
+```
+
+
+# Discussion
+
+## Identification of Velocity parameter
+
+Both calibrations overestimate this parameter:
+
+```{r}
+mVelocity <- matrix(c(Velocity,
+                      OutputsCalibDown1$ParamFinalR[1],
+                      OutputsCalibDown2$ParamFinalR[1],
+                      OutputsCalibDown3$ParamFinalR[1]),
+                    ncol = 1,
+                    dimnames = list(c("theoretical",
+                                      "calibrated with observed upstream flow",
+                                      "calibrated with simulated  upstream flow",
+                                      "calibrated with sim upstream flow and regularisation"),
+                                    c("Velocity parameter")))
+knitr::kable(mVelocity)
+```
+
+## Value of the performance criteria with theoretical calibration
+
+Theoretically, the parameters of the downstream GR4J model should be the same as the upstream one with the velocity as extra parameter:
+
+```{r}
+OutputsModelDownTheo <- RunModel(InputsModel = InputsModelDown2,
+                                 RunOptions = RunOptionsDown,
+                                 Param = ParamDownTheo,
+                                 FUN_MOD = RunModel_GR4J)
+KGE_downTheo <- ErrorCrit_KGE(InputsCritDown, OutputsModelDownTheo)
+```
+
+
+## Parameters and performance of each subcatchment for all calibrations
+
+```{r}
+comp <- matrix(c(0, OutputsCalibUp$ParamFinalR,
+                 rep(OutputsCalibDown1$ParamFinalR, 2),
+                 OutputsCalibDown2$ParamFinalR,
+                 OutputsCalibDown3$ParamFinalR,
+                 ParamDownTheo),
+               ncol = 5, byrow = TRUE)
+comp <- cbind(comp, c(OutputsCalibUp$CritFinal,
+                      OutputsCalibDown1$CritFinal,
+                      KGE_down1$CritValue,
+                      OutputsCalibDown2$CritFinal,
+                      KGE_down3$CritValue,
+                      KGE_downTheo$CritValue))
+colnames(comp) <- c("Velocity", paste0("X", 1:4), "KGE(√Q)")
+rownames(comp) <- c("Calibration of the upstream subcatchment",
+                    "Calibration 1 with observed upstream flow",
+                    "Validation 1 with simulated upstream flow",
+                    "Calibration 2 with simulated upstream flow",
+                    "Calibration 3 with simulated upstream flow and regularisation",
+                    "Validation theoretical set of parameters")
+knitr::kable(comp)
+```
+
+Even if calibration with observed upstream flows gives an improved performance criteria, in validation using simulated upstream flows the result is quite similar as the performance obtained with the calibration with upstream simulated flows. The theoretical set of parameters give also an equivalent performance but still underperforming the calibration 2 one. Regularisation allows to get similar performance as the one for calibration with simulated flows but with the big advantage of having parameters closer to the theoretical ones (Especially for the velocity parameter).
+
+# References

vignettes/fig/airGR_features_diagram.svg

+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+<svg xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:svg="http://www.w3.org/2000/svg" xmlns="http://www.w3.org/2000/svg" xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd" xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape" id="svg2" version="1.1" inkscape:version="0.91 r13725" xml:space="preserve" width="900" height="675" viewBox="0 0 900 675" sodipodi:docname="airGR_features.svg"><metadata id="metadata8"><rdf:RDF><cc:Work rdf:about=""><dc:format>image/svg+xml</dc:format><dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/><dc:title/></cc:Work></rdf:RDF></metadata><defs id="defs6"><clipPath clipPathUnits="userSpaceOnUse" id="clipPath16"><path d="m 0,1.2207e-4 720,0 0,539.99999793 -720,0 L 0,1.2207e-4 Z" id="path18" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath26"><path d="M 3.08e-7,1.2207e-4 720,1.2207e-4 720,540.00012 3.08e-7,540.00012 l 0,-539.99999793 z" id="path28" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath102"><path d="M -2.645e-6,540 720,540 720,1.2207e-4 -6.582e-6,1.2207e-4" id="path104" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath116"><path d="M -1.161e-6,540 720,540 720,1.2207e-4 l -720.000001161,0" id="path118" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath132"><path d="M -2.728e-6,540 720,540 720,1.2207e-4 l -720.000006665,0" id="path134" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath178"><path d="M -1.52e-6,540 720,540 720,6.104e-6 -1.52e-6,6.104e-6" id="path180" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath></defs><sodipodi:namedview pagecolor="#ffffff" bordercolor="#666666" borderopacity="1" objecttolerance="10" gridtolerance="10" guidetolerance="10" inkscape:pageopacity="0" inkscape:pageshadow="2" inkscape:window-width="1920" inkscape:window-height="1132" id="namedview4" showgrid="false" inkscape:zoom="0.98890193" inkscape:cx="333.64663" inkscape:cy="425.84952" inkscape:window-x="-8" inkscape:window-y="-8" inkscape:window-maximized="1" inkscape:current-layer="g10"/><g id="g10" inkscape:groupmode="layer" inkscape:label="airGR_features" transform="matrix(1.25,0,0,-1.25,0,675)"><g id="g12"><g id="g14" clip-path="url(#clipPath16)"><path d="m 0,6.104e-6 720,0 L 720,540.00001 l -720,0 L 0,6.104e-6 Z" style="fill:#ffffff;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path20" inkscape:connector-curvature="0"/></g></g><g id="g22"><g id="g24" clip-path="url(#clipPath26)"><path d="M 391.5,540 390.85,-2.1562" style="fill:none;stroke:#7f7f7f;stroke-width:2;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:8, 6;stroke-dashoffset:0;stroke-opacity:1" id="path30" inkscape:connector-curvature="0"/></g></g><path d="m 308.97,304.02 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.65 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.97 -6.62,6.62 l 0,26.46 z" style="fill:#00b0f0;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path32" inkscape:connector-curvature="0"/><g id="g34"><path d="m 308.97,304.02 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.65 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.97 -6.62,6.62 l 0,26.46 z" style="fill:none;stroke:#00b0f0;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path36" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,320.57,284.74)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text38"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 52.416 61.866001 71.351997 80.802002 86.832001 93.870003 109.26 118.746 128.196 137.214" y="0" sodipodi:role="line" id="tspan40">CreateInputsModel</tspan></text>
+<path d="m 481.47,473.17 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.61,-6.61 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.61 l 0,26.46 z" style="fill:#cc00cc;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path42" inkscape:connector-curvature="0"/><g id="g44"><path d="m 481.47,473.17 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.61,-6.61 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.61 l 0,26.46 z" style="fill:none;stroke:#cc00cc;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path46" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,493.2,453.94)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text48"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 57.438 66.096001 70.181999 74.267998 83.718002 95.634003 104.994 111.024 115.074 124.56 134.00999" y="0" sodipodi:role="line" id="tspan50">CreateCalibOptions</tspan></text>
+<path d="m 481.5,400.38 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.62 -6.62,-6.62 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.62 l 0,26.45 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path52" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,494.81,381.12)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text54"><tspan x="0 9.5939999 18.216 22.356001 26.406 35.855999 41.759998 50.273998 56.304001 60.354 69.839996 79.290001 88.290001 103.68 107.73 115.29 124.74 133.758" y="0" sodipodi:role="line" id="tspan56">Calibration_Michel</tspan></text>
+<path d="m 308.97,131.92 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.618 -6.62,-6.618 l -151.2,0 c -3.65,0 -6.61,2.958 -6.61,6.618 l 0,26.45 z" style="fill:#92d050;fill-opacity:0.30196001;fill-rule:evenodd;stroke:none" id="path58" inkscape:connector-curvature="0"/><g id="g60"><path d="m 308.97,131.92 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.618 -6.62,-6.618 l -151.2,0 c -3.65,0 -6.61,2.958 -6.61,6.618 l 0,26.45 z" style="fill:none;stroke:#92d050;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path62" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,331.25,112.61)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text64"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 52.416 61.866001 71.351997 80.802002 86.832001 93.870003 103.464 109.746 113.832" y="0" sodipodi:role="line" id="tspan66">CreateInputsCrit</tspan></text>
+<path d="m 31.145,132.03 c 0,3.66 2.961,6.62 6.615,6.62 l 49.139,0 c 3.653,0 6.615,-2.96 6.615,-6.62 l 0,-26.46 c 0,-3.65 -2.962,-6.61 -6.615,-6.61 l -49.139,0 c -3.654,0 -6.615,2.96 -6.615,6.61 l 0,26.46 z" style="fill:#92d050;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path68" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,47.784,112.73)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text70"><tspan x="0 9.4499998 13.59 23.076" y="0" sodipodi:role="line" id="tspan72">plot</tspan></text>
+<path d="m 481.5,385.53 -23.9,0 c -0.89,0 -1.62,0.73 -1.62,1.62 l 0,0.03 1.62,-1.63 -15.77,0 0,3.25 15.77,0 c 0.9,0 1.63,-0.72 1.63,-1.62 l 0,-0.03 -1.63,1.63 23.9,0 0,-3.25 z m -38.04,-3.23 -9.75,4.88 9.75,4.87 0,-9.75 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path74" inkscape:connector-curvature="0"/><path d="m 308.97,120.32 -14.17,0 c -0.9,0 -1.63,-0.73 -1.63,-1.63 l 0,-0.09 1.63,1.62 -6.05,0 0,-3.25 6.05,0 c 0.89,0 1.62,0.73 1.62,1.63 l 0,0.09 -1.62,-1.62 14.17,0 0,3.25 z m -18.6,3.15 -9.75,-4.87 9.75,-4.88 0,9.75 z" style="fill:#92d050;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path76" inkscape:connector-curvature="0"/><path d="m 120.22,388.8 -57.891,0 c -0.897,0 -1.625,-0.72 -1.625,-1.62 l 0,-240.41 3.25,0 0,240.41 -1.625,-1.63 57.891,0 0,3.25 z m -62.766,-240.4 4.875,-9.75 4.875,9.75 -9.75,0 z" style="fill:#7030a0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path78" inkscape:connector-curvature="0"/><path d="m 308.97,289.16 -77.96,0 c -0.9,0 -1.62,0.73 -1.62,1.63 l 0,68.39 3.25,0 0,-68.39 -1.63,1.62 77.96,0 0,-3.25 z m -82.83,68.4 4.9,9.74 4.85,-9.76 -9.75,0.02 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path80" inkscape:connector-curvature="0"/><path d="m 308.97,220.52 -78.19,0 c -0.89,0 -1.62,0.72 -1.62,1.62 l 0,137.07 3.25,0 0,-137.07 -1.63,1.63 78.19,0 0,-3.25 z m -83.06,137.06 4.87,9.75 4.88,-9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path82" inkscape:connector-curvature="0"/><path d="m 473.4,220.52 90.31,0 c 0.9,0 1.63,0.72 1.63,1.62 l 0,137.05 -3.25,0 0,-137.05 1.62,1.63 -90.31,0 0,-3.25 z m 95.19,137.04 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path84" inkscape:connector-curvature="0"/><path d="m 365.67,51.027 -303.341,0 c -0.897,0 -1.625,0.728 -1.625,1.625 l 0,38.183 3.25,0 0,-38.183 -1.625,1.625 303.341,0 0,-3.25 z M 57.454,89.21 l 4.875,9.75 4.875,-9.75 -9.75,0 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path86" inkscape:connector-curvature="0"/><path d="m 473.4,289.16 90.31,0 c 0.9,0 1.63,0.73 1.63,1.63 l 0,68.4 -3.25,0 0,-68.4 1.62,1.62 -90.31,0 0,-3.25 z m 95.19,68.4 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path88" inkscape:connector-curvature="0"/><path d="m 342.99,387.18 9.92,19.84 70.88,0 9.92,-19.84 -9.92,-19.85 -70.88,0 -9.92,19.85 z" style="fill:#00508c;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path90" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,365.11,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text92"><tspan x="0 8.8739996 17.496 23.4 32.021999" y="0" sodipodi:role="line" id="tspan94">Param</tspan></text>
+<path d="m 564.12,173.36 -172.94,0 c -0.89,0 -1.62,-0.73 -1.62,-1.63 l 0,-25.07 3.25,0 0,25.07 -1.63,-1.62 172.94,0 0,3.25 z m -177.81,-25.07 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path96" inkscape:connector-curvature="0"/><g id="g98"><g id="g100" clip-path="url(#clipPath102)"><path d="m 562.06,440.1 0,-16.55 c 0,-0.9 0.72,-1.63 1.62,-1.63 l 0.03,0 -1.62,1.63 0,-8.43 3.25,0 0,8.43 c 0,0.9 -0.73,1.62 -1.63,1.62 l -0.03,0 1.63,-1.62 0,16.55 -3.25,0 z m -3.22,-23.35 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path106" inkscape:connector-curvature="0"/></g></g><path d="m 473.4,117.07 124.74,0 c 0.89,0 1.62,0.73 1.62,1.62 l 0,240.49 -3.25,0 0,-240.49 1.63,1.63 -124.74,0 0,-3.25 z m 129.61,240.48 -4.87,9.75 -4.88,-9.75 9.75,0 z" style="fill:#92d050;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path108" inkscape:connector-curvature="0"/><path d="m 354.33,388.8 -34.84,0 c -0.9,0 -1.63,-0.72 -1.63,-1.62 l 0,0 1.63,1.62 -26.72,0 0,-3.25 26.72,0 c 0.9,0 1.62,0.73 1.62,1.63 l 0,0 -1.62,-1.63 34.84,0 0,3.25 z m -59.93,3.25 -9.75,-4.87 9.75,-4.88 0,9.75 z" style="fill:#00508c;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path110" inkscape:connector-curvature="0"/><g id="g112"><g id="g114" clip-path="url(#clipPath116)"><path d="m 200.03,367.33 0,-114.44 c 0,-0.9 -0.73,-1.63 -1.62,-1.63 l 0,0 1.62,1.63 0,-106.32 -3.25,0 0,106.32 c 0,0.89 0.73,1.62 1.63,1.62 l 0,0 -1.63,-1.62 0,114.44 3.25,0 z m 3.25,-219.14 -4.87,-9.75 -4.88,9.75 9.75,0 z" style="fill:#7030a0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path120" inkscape:connector-curvature="0"/></g></g><path d="m 120.22,400.41 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.62 l 0,26.46 z" style="fill:#7030a0;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path122" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,164.42,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text124"><tspan x="0 9.7740002 19.224001 28.674 44.063999 53.549999 63 71.963997" y="0" sodipodi:role="line" id="tspan126">RunModel</tspan></text>
+<g id="g128"><g id="g130" clip-path="url(#clipPath132)"><path d="m 389.56,270.94 0,-14.47 c 0,-0.9 0.73,-1.63 1.62,-1.63 l 0,0 -1.62,1.63 0,-6.36 3.25,0 0,6.36 c 0,0.89 -0.73,1.62 -1.63,1.62 l 0,0 1.63,-1.62 0,14.47 -3.25,0 z m -3.25,-19.2 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path136" inkscape:connector-curvature="0"/></g></g><text transform="matrix(1,0,0,-1,509.9,509.71)" style="font-variant:normal;font-weight:bold;font-size:24px;font-family:Calibri;-inkscape-font-specification:'Calibri Bold';writing-mode:lr-tb;fill:#5f5f5f;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text138"><tspan x="0 12.696 24.552 30.455999 36.360001 49.248001 57.240002 68.879997 77.208 83.112 96.024002" y="0" sodipodi:role="line" id="tspan140">Calibration</tspan></text>
+<text transform="matrix(1,0,0,-1,142.75,509.38)" style="font-variant:normal;font-weight:bold;font-size:24px;font-family:Calibri;-inkscape-font-specification:'Calibri Bold';writing-mode:lr-tb;fill:#5f5f5f;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text142"><tspan x="0 11.352 17.256001 36.84 49.728001 55.632 67.199997 75.528 81.431999 94.344002" y="0" sodipodi:role="line" id="tspan144">Simulation</tspan></text>
+<path d="m 308.97,235.37 c 0,3.65 2.96,6.62 6.62,6.62 l 151.19,0 c 3.66,0 6.62,-2.97 6.62,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.62,-6.61 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.61 l 0,26.46 z" style="fill:#00b0f0;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path146" inkscape:connector-curvature="0"/><g id="g148"><path d="m 308.97,235.37 c 0,3.65 2.96,6.62 6.62,6.62 l 151.19,0 c 3.66,0 6.62,-2.97 6.62,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.62,-6.61 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.61 l 0,26.46 z" style="fill:none;stroke:#00b0f0;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path150" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,324.29,216.07)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text152"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 57.618 67.068001 76.517998 88.433998 97.793999 103.824 107.874 117.36 126.81" y="0" sodipodi:role="line" id="tspan154">CreateRunOptions</tspan></text>
+<path d="m 116.19,400.41 c 0,3.65 2.96,6.61 6.62,6.61 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.62 -6.61,-6.62 l -151.2,0 c -3.66,0 -6.62,2.96 -6.62,6.62 l 0,26.46 z" style="fill:#7030a0;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path156" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,160.39,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text158"><tspan x="0 9.7740002 19.224001 28.674 44.063999 53.549999 63 71.963997" y="0" sodipodi:role="line" id="tspan160">RunModel</tspan></text>
+<path d="m 116.19,131.83 c 0,3.65 2.96,6.61 6.62,6.61 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.618 -6.61,-6.618 l -151.2,0 c -3.66,0 -6.62,2.958 -6.62,6.618 l 0,26.46 z" style="fill:#92d050;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path162" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,167.11,112.51)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text164"><tspan x="0 8.7840004 14.994 20.988001 30.474001 36.756001 46.349998 52.542 56.627998" y="0" sodipodi:role="line" id="tspan166">ErrorCrit</tspan></text>
+<path d="m 351.85,52.652 9.92,19.845 59.53,0 9.93,-19.845 -9.93,-19.844 -59.53,0 -9.92,19.844 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path168" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,372.5,46.56)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text170"><tspan x="0 12.114 21.6 30.959999" y="0" sodipodi:role="line" id="tspan172">Qobs</tspan></text>
+<g id="g174"><g id="g176" clip-path="url(#clipPath178)"><path d="m 392.81,71.77 0,13.541 c 0,0.897 -0.73,1.625 -1.63,1.625 l 0,0 1.63,-1.625 0,5.416 -3.25,0 0,-5.416 c 0,-0.897 0.73,-1.625 1.62,-1.625 l 0,0 -1.62,1.625 0,-13.541 3.25,0 z m 3.25,17.332 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path182" inkscape:connector-curvature="0"/></g></g><path d="m 563.71,222.75 0.62,-50.29" style="fill:none;stroke:#00b0f0;stroke-width:3.25;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:none;stroke-opacity:1" id="path184" inkscape:connector-curvature="0"/></g></svg>
\ No newline at end of file