tests/testthat/test-SeriesAggreg.R

+context("SeriesAggreg")
+
+## load catchment data
+data(L0123002)
+
+# Test removed because of #171 to reintegrated later...
+# test_that("No warning with InputsModel Cemaneige'", {
+#   ## preparation of the InputsModel object
+#   InputsModel <- CreateInputsModel(
+#     FUN_MOD = RunModel_CemaNeige,
+#     DatesR = BasinObs$DatesR,
+#     Precip = BasinObs$P,
+#     TempMean = BasinObs$T,
+#     ZInputs = BasinInfo$HypsoData[51],
+#     HypsoData = BasinInfo$HypsoData,
+#     NLayers = 5
+#   )
+#   # Expect no warning: https://stackoverflow.com/a/33638939/5300212
+#   expect_warning(SeriesAggreg(InputsModel, Format = "%m"),
+#                  regexp = NA)
+# })
+
+
+test_that("Warning: deprecated 'TimeFormat' argument", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E
+  )
+  expect_warning(SeriesAggreg(InputsModel, Format = "%Y%m", TimeFormat = "daily"),
+                 regexp = "deprecated 'TimeFormat' argument")
+})
+
+
+test_that("Warning: deprecated 'NewTimeFormat' argument: please use 'Format' instead",
+          {
+            InputsModel <- CreateInputsModel(
+              FUN_MOD = RunModel_GR4J,
+              DatesR = BasinObs$DatesR,
+              Precip = BasinObs$P,
+              PotEvap = BasinObs$E
+            )
+            expect_warning(SeriesAggreg(InputsModel, NewTimeFormat = "monthly"),
+                           regexp = "deprecated 'NewTimeFormat' argument: please use 'Format' instead")
+          })
+
+
+test_that("Warning: deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
+          {
+            InputsModel <- CreateInputsModel(
+              FUN_MOD = RunModel_GR4J,
+              DatesR = BasinObs$DatesR,
+              Precip = BasinObs$P,
+              PotEvap = BasinObs$E
+            )
+            expect_warning(SeriesAggreg(InputsModel, Format = "%Y%m", NewTimeFormat = "monthly"),
+                           regexp = "deprecated 'NewTimeFormat' argument: 'Format' argument is used instead")
+          })
+
+
+test_that("Check SeriesAggreg output values on yearly aggregation", {
+  TabSeries <- data.frame(
+    DatesR = BasinObs$DatesR,
+    P = BasinObs$P,
+    E = BasinObs$E,
+    Qmm = BasinObs$Qmm
+  )
+  GoodValues <- apply(BasinObs[BasinObs$DatesR >= as.POSIXct("1984-09-01", tz = "UTC") &
+                                 BasinObs$DatesR < as.POSIXct("1985-09-01", tz = "UTC"),
+                               c("P", "E", "Qmm")],
+                      MARGIN = 2, FUN = sum)
+  TestedValues <- unlist(SeriesAggreg(TabSeries,
+                                      Format = "%Y",
+                                      YearFirstMonth = 9,
+                                      ConvertFun = rep("sum", 3))[1, c("P", "E", "Qmm")])
+  expect_equal(GoodValues, TestedValues)
+})
+
+
+test_that("Regime calculation should switch ConvertFun to 'mean' for InputsModel", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E
+  )
+  expect_equal(SeriesAggreg(InputsModel, Format = "%m")$Precip,
+               SeriesAggreg(BasinObs[, c("DatesR", "P")], Format = "%m", ConvertFun = "mean")$P)
+})
+
+
+test_that("No DatesR should warning", {
+  TabSeries <- list(
+    Dates = BasinObs$DatesR,
+    P = BasinObs$P,
+    E = BasinObs$E,
+    Qmm = BasinObs$Qmm
+  )
+  expect_warning(
+    SeriesAggreg(TabSeries, Format = "%Y%m", ConvertFun = "sum"),
+    regexp = "has been automatically chosen"
+  )
+})
+
+
+test_that("Check SeriesAggreg.list 'DatesR' argument", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E
+  )
+  DatesR <- InputsModel$DatesR
+  # No InputsModel$DatesR
+  InputsModel$DatesR <- NULL
+  expect_error(SeriesAggreg(InputsModel, Format = "%Y%m"), regexp = "'POSIXt'")
+  # Other list item chosen
+  InputsModel$SuperDates <- DatesR
+  expect_warning(SeriesAggreg(InputsModel, Format = "%Y%m"), regexp = "SuperDates")
+  # Wrong InputsModel$DatesR
+  InputsModel$DatesR <- BasinObs$P
+  expect_error(SeriesAggreg(InputsModel, Format = "%Y%m"), regexp = "'POSIXt'")
+
+})
+
+
+test_that("Check SeriesAggreg.list with embedded lists", {
+  InputsModel <-
+    CreateInputsModel(
+      FUN_MOD = RunModel_CemaNeige,
+      DatesR = BasinObs$DatesR,
+      Precip = BasinObs$P,
+      TempMean = BasinObs$T,
+      ZInputs = BasinInfo$HypsoData[51],
+      HypsoData = BasinInfo$HypsoData,
+      NLayers = 5
+    )
+  I2 <- SeriesAggreg(InputsModel, Format = "%Y%m")
+  expect_equal(length(I2$ZLayers), 5)
+  expect_null(I2$LayerPrecip$DatesR)
+  expect_equal(length(I2$DatesR), length(I2$LayerPrecip$L1))
+})
+
+
+test_that("Check SeriesAggreg.outputsModel", {
+  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
+  )
+
+  ## run period selection
+  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"))
+
+  ## preparation of the RunOptions object
+  suppressWarnings(
+    RunOptions <- CreateRunOptions(
+      FUN_MOD = RunModel_CemaNeigeGR4J,
+      InputsModel = InputsModel,
+      IndPeriod_Run = Ind_Run
+    )
+  )
+
+  ## simulation
+  Param <- c(
+    X1 = 408.774,
+    X2 = 2.646,
+    X3 = 131.264,
+    X4 = 1.174,
+    CNX1 = 0.962,
+    CNX2 = 2.249
+  )
+  OutputsModel <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
+                                         RunOptions = RunOptions,
+                                         Param = Param)
+
+  O2 <- SeriesAggreg(OutputsModel, Format = "%Y%m")
+  expect_equal(length(O2$StateEnd), 3)
+  expect_equal(length(O2$DatesR),
+               length(O2$CemaNeigeLayers$Layer01$Pliq))
+})
+
+
+test_that("Check data.frame handling in SeriesAggreg.list", {
+  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 and the downstream one in km
+    BasinAreas = c(180, 180) # Upstream and downstream areas in km²
+  )
+  # Test removed because of #171 to reintegrated later...
+  # expect_warning(SeriesAggreg(InputsModelDown1, Format = "%Y%m"),
+  #                regexp = NA)
+  I2 <- SeriesAggreg(InputsModelDown1, Format = "%Y%m")
+  expect_equal(length(I2$DatesR), nrow(I2$Qupstream))
+  InputsModelDown1$Qupstream <- InputsModelDown1$Qupstream[-1, , drop = FALSE]
+  expect_warning(SeriesAggreg(InputsModelDown1, Format = "%Y%m"),
+                 regexp = "it will be ignored in the aggregation")
+})
+
+
+test_that("SeriesAggreg from and to the same time step should return initial time series", {
+  InputsModel <- CreateInputsModel(
+    FUN_MOD = RunModel_GR4J,
+    DatesR = BasinObs$DatesR,
+    Precip = BasinObs$P,
+    PotEvap = BasinObs$E
+  )
+  I2 <- SeriesAggreg(InputsModel, Format = "%Y%m")
+  expect_warning(SeriesAggreg(I2, Format = "%Y%m"), regexp = "No time-step conversion was performed")
+  expect_equal(I2, suppressWarnings(SeriesAggreg(I2, Format = "%Y%m")))
+})
+
+
+test_that("SeriesAggreg.data.frame with first column not named DatesR should work", {
+            expect_warning(SeriesAggreg(
+              data.frame(BasinObs$DatesR, BasinObs$Qmm),
+              Format = "%Y%m",
+              ConvertFun = "sum"
+            ),
+            regexp = NA)
+})
+
+
+test_that("SeriesAggreg should work with ConvertFun 'min', 'max' and 'median'", {
+  Qls <- BasinObs[, c("DatesR", "Qls")]
+  test_ConvertFunRegime <- function(x, ConvertFun, Format) {
+    expect_equal(nrow(SeriesAggreg(x, Format, ConvertFun = ConvertFun)),
+                 length(unique(format(BasinObs$DatesR, format = "%Y"))))
+  }
+  lapply(c("max", "min", "median"), function(x) {test_ConvertFunRegime(Qls, x, Format = "%Y")})
+})
+
+
+test_that("Error on convertFun Q without 0-100", {
+  Qls <- BasinObs[, c("DatesR", "Qls")]
+  expect_error(SeriesAggreg(Qls, Format = "%Y", "q101"))
+  expect_error(SeriesAggreg(Qls, Format = "%Y", "q-2"))
+  expect_error(SeriesAggreg(Qls, Format = "%Y", "q12.5"))
+})
+
+
+test_that("ConvertFun q50 should be equal to median", {
+  Qls <- BasinObs[, c("DatesR", "Qls")]
+  expect_equal(SeriesAggreg(Qls, Format = "%Y", "q50"),
+               SeriesAggreg(Qls, Format = "%Y", "median"))
+  expect_equal(SeriesAggreg(Qls, Format = "%Y", "q50"),
+               SeriesAggreg(Qls, Format = "%Y", "q050"))
+})
+

tests/testthat/test-evap.R

+context("Test evaporation")
+
+
+rm(list = ls())
+data(L0123001); BasinObs_L0123001 <- BasinObs
+data(L0123002); BasinObs_L0123002 <- BasinObs
+data(L0123003); BasinObs_L0123003 <- BasinObs
+
+comp_evap <- function(BasinObs,
+                      Lat, LatUnit,
+                      TimeStepIn = "daily",
+                      TimeStepOut = "daily") {
+  PotEvap <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+                      Temp = BasinObs$T,
+                      Lat = Lat, LatUnit = LatUnit,
+                      TimeStepIn = TimeStepIn, TimeStepOut = TimeStepOut)
+  PotEvapFor <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+                         Temp = BasinObs$T,
+                         Lat = Lat, LatUnit = LatUnit,
+                         TimeStepIn = TimeStepIn, TimeStepOut = TimeStepOut,
+                         RunFortran = TRUE)
+  all(range(PotEvap - PotEvapFor) < 0.000001)
+}
+
+
+test_that("PE_Oudin works", {
+  skip_on_cran()
+
+  expect_true(comp_evap(BasinObs = BasinObs_L0123001,
+                        Lat = 0.8, LatUnit = "rad",
+                        TimeStepIn = "daily", TimeStepOut = "daily"))
+  expect_true(comp_evap(BasinObs = BasinObs_L0123001,
+                        Lat = 0.8, LatUnit = "rad",
+                        TimeStepIn = "daily", TimeStepOut = "hourly"))
+  expect_true(comp_evap(BasinObs = BasinObs_L0123002,
+                        Lat = 0.9, LatUnit = "rad",
+                        TimeStepIn = "daily", TimeStepOut = "daily"))
+  expect_true(comp_evap(BasinObs = BasinObs_L0123002,
+                        Lat = 0.9, LatUnit = "rad",
+                        TimeStepIn = "daily", TimeStepOut = "hourly"))
+
+  ## check with several catchments using different values for Lat
+
+  ## one by one
+  PotEvapFor1 <- PE_Oudin(JD = as.POSIXlt(BasinObs_L0123001$DatesR)$yday + 1,
+                          Temp = BasinObs_L0123001$T,
+                          Lat = 0.8, LatUnit = "rad",
+                          RunFortran = TRUE)
+  PotEvapFor2 <- PE_Oudin(JD = as.POSIXlt(BasinObs_L0123002$DatesR)$yday + 1,
+                          Temp = BasinObs_L0123002$T,
+                          Lat = 0.9, LatUnit = "rad",
+                          RunFortran = TRUE)
+
+  ## all in one
+  BasinObs_L0123001$Lat <- 0.8
+  BasinObs_L0123002$Lat <- 0.9
+  BasinObs <- rbind(BasinObs_L0123001, BasinObs_L0123002)
+  PotEvapFor <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+                         Temp = BasinObs$T,
+                         Lat = BasinObs$Lat, LatUnit = "rad",
+                         RunFortran = TRUE)
+
+  expect_equal(PotEvapFor, c(PotEvapFor1, PotEvapFor2))
+})
+
+test_that("Inconsitent time series", {
+  skip_on_cran()
+
+  msgDaily <- "each day should have only one identical value of Julian days. The time series is not sorted, or contains duplicate or missing dates"
+  msgHoury <- "each day must have 24 identical values of Julian days (one for each hour). The time series is not sorted, or contains duplicate or missing dates"
+
+  # duplicated dates
+  DatesFor1Dupl <- BasinObs_L0123001$DatesR
+  DatesFor1Dupl[5L] <- DatesFor1Dupl[4L]
+  expect_warning(object = PE_Oudin(JD = as.POSIXlt(DatesFor1Dupl)$yday + 1,
+                                   Temp = BasinObs_L0123001$T,
+                                   Lat = 0.8, LatUnit = "rad"),
+                 regexp = msgDaily)
+
+  # not ordered daily dates
+  DatesFor1Messy <- sample(BasinObs_L0123001$DatesR)
+  expect_warning(object = PE_Oudin(JD = as.POSIXlt(DatesFor1Messy)$yday + 1,
+                                   Temp = BasinObs_L0123001$T,
+                                   Lat = 0.8, LatUnit = "rad"),
+                 regexp = msgDaily)
+
+  # not ordered hourly dates
+  DatesFor3Messy <- sample(BasinObs_L0123003$DatesR)
+  expect_error(object = PE_Oudin(JD = as.POSIXlt(DatesFor3Messy)$yday + 1,
+                                   Temp = seq_along(BasinObs_L0123003$T),
+                                   Lat = 0.8, LatUnit = "rad", TimeStepIn = "hourly"),
+                 regexp = msgHoury, fixed = TRUE)
+
+
+})
+
+

tests/testthat/test-vignettes.R

+context("Test vignette chunks")
+
+test_that("V01_get_started works", {
+  skip_on_cran()
+  rm(list = ls())
+  expect_true(RunVignetteChunks("V01_get_started"))
+  TestQmmQlsConversion(BasinObs, BasinInfo$BasinArea)
+})
+
+test_that("V02.1_param_optim works", {
+  skip_on_cran()
+  rm(list = ls())
+  load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+  load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
+  rda_resGLOB <- resGLOB
+  rda_resPORT <- resPORT
+  rda_optMO <- optMO
+  expect_true(RunVignetteChunks("V02.1_param_optim"))
+  expect_equal(summary(resGLOB), summary(rda_resGLOB), tolerance = 1e-7)
+  expect_equal(resGLOB[, -1], rda_resGLOB[, -1], tolerance = 1e-2) # High tolerance due to randomisation in optimisations
+  expect_equal(summary(optMO$parameters), summary(rda_optMO$parameters), tolerance = 1e-7)
+})
+
+test_that("V02.2_param_mcmc works", {
+  skip_on_cran()
+  rm(list = ls())
+  load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))
+  rda_gelRub <- gelRub
+  rda_multDRAM <- multDRAM
+  expect_true(RunVignetteChunks("V02.2_param_mcmc"))
+  expect_equal(gelRub, rda_gelRub, tolerance = 1e-7)
+  expect_equal(multDRAM, rda_multDRAM, tolerance = 1e-7)
+})
+
+test_that("V03_param_sets_GR4J works", {
+  skip_on_cran()
+  rm(list = ls())
+  expect_true(RunVignetteChunks("V03_param_sets_GR4J"))
+})
+
+test_that("V04_cemaneige_hysteresis works", {
+  skip_on_cran()
+  rm(list = ls())
+  load(system.file("vignettesData/vignetteCNHysteresis.rda", package = "airGR"))
+  rda_OutputsCrit_Cal <- OutputsCrit_Cal
+  rda_OutputsCrit_Cal_NoHyst <- OutputsCrit_Cal_NoHyst
+  rda_OutputsCrit_Val <- OutputsCrit_Val
+  rda_OutputsCrit_Val_NoHyst <- OutputsCrit_Val_NoHyst
+  expect_true(RunVignetteChunks("V04_cemaneige_hysteresis"))
+  TestQmmQlsConversion(BasinObs, BasinInfo$BasinArea)
+  expect_equal(OutputsCrit_Cal, rda_OutputsCrit_Cal, tolerance = 1e-7)
+  expect_equal(OutputsCrit_Cal_NoHyst, rda_OutputsCrit_Cal_NoHyst, tolerance = 1e-7)
+  expect_equal(OutputsCrit_Val, rda_OutputsCrit_Val, tolerance = 1e-7)
+  expect_equal(OutputsCrit_Val_NoHyst, rda_OutputsCrit_Val_NoHyst, tolerance = 1e-7)
+})

vignettes/V00_airgr_ref.bib

 
+@phdthesis{ficchi_adaptive_2017,
+	title = {An adaptive hydrological model for multiple time-steps: diagnostics and improvements based on fluxes consistency},
+	shorttitle = {An adaptive hydrological model for multiple time-steps},
+	url = {https://www.theses.fr/2017PA066097},
+	abstract = {Cette thèse vise à explorer la question du changement d'échelle temporelle en modélisation hydrologique conceptuelle. Les principaux objectifs sont : (i) étudier les effets du changement du pas de temps sur les performances, les paramètres et la structure des modèles hydrologiques ; (ii) mettre au point un modèle pluie-débit applicable à différents pas de temps. Notre point de départ est le modèle global journalier GR4J, développé à Irstea. Ce modèle a été choisi comme le modèle de référence à adapter à d'autres résolutions plus fines, jusqu'à des pas de temps infra-horaires, en suivant une approche descendante. Pour nos tests, nous avons construit une base de données de 240 bassins versants non influencés en France, à différents pas de temps allant de 6 minutes à 1 jour, en utilisant: (i) les données pluviométriques à 6 minutes et la réanalyse des lames d'eau journalières à plus haute résolution spatiale ; (ii) les données de température journalière pour le calcul de l'évapotranspiration potentielle ; (iii) les données hydrométriques à pas de temps variable. Nous avons étudié l'impact de la distribution temporelle des entrées sur les performances du modèle en se focalisant sur la simulation de crue, sur la base de 2400 événements. Ensuite, notre évaluation du modèle a porté sur l'analyse de la cohérence des flux internes du modèle à différents pas de temps, afin d'assurer une performance satisfaisante à travers un fonctionnement du modèle cohérent. Notre diagnostic du modèle nous a permis d'identifier une amélioration de la structure du modèle à différents pas de temps infra-journaliers basée sur la complexification de la composante d'interception du modèle.},
+	urldate = {2017-11-24},
+	school = {Université Pierre et Marie Curie, Paris 6},
+	author = {Ficchi, Andrea},
+	month = feb,
+	year = {2017},
+	keywords = {Rainfall-runoff modelling, airGRcite, Cohérence structurelle, Événements de crue, Interception, Modèles GR, Modélisation pluie-Débit, Pas de temps fin, Short time step}
+}
+
+@article{ficchi_hydrological_2019,
+	title = {Hydrological modelling at multiple sub-daily time steps: model improvement via flux-matching},
+	issn = {00221694},
+	shorttitle = {Hydrological modelling at multiple sub-daily time steps},
+	url = {https://linkinghub.elsevier.com/retrieve/pii/S0022169419305281},
+	doi = {10.1016/j.jhydrol.2019.05.084},
+	language = {en},
+	urldate = {2019-06-11},
+	journal = {Journal of Hydrology},
+	author = {Ficchì, Andrea and Perrin, Charles and Andréassian, Vazken},
+	month = jun,
+	year = {2019},
+	keywords = {airGR, airGRcite}
+}
+
 @phdthesis{mathevet_quels_2005,
 	address = {Paris},
 	title = {Quels modèles pluie-débit globaux au pas de temps horaire ? {Développements} empiriques et comparaison de modèles sur un large échantillon de bassins versants},
@@ -150,4 +178,25 @@
 	year = {2019},
 	keywords = {airGRcite},
 	pages = {70--81}
-}
\ No newline at end of file
+}
+
+@article{delavenne_regularization_2019,
+  ids = {lavenneRegularizationApproachImprove2019a},
+  title = {A {{Regularization Approach}} to {{Improve}} the {{Sequential Calibration}} of a {{Semidistributed Hydrological Model}}},
+  author = {de Lavenne, Alban and Andréassian, Vazken and Thirel, Guillaume and Ramos, Maria-Helena and Perrin, Charles},
+  date = {2019},
+  journaltitle = {Water Resources Research},
+  volume = {55},
+  pages = {8821--8839},
+  issn = {1944-7973},
+  doi = {10.1029/2018WR024266},
+  url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018WR024266},
+  urldate = {2021-06-07},
+  abstract = {In semidistributed hydrological modeling, sequential calibration usually refers to the calibration of a model by considering not only the flows observed at the outlet of a catchment but also the different gauging points inside the catchment from upstream to downstream. While sequential calibration aims to optimize the performance at these interior gauged points, we show that it generally fails to improve performance at ungauged points. In this paper, we propose a regularization approach for the sequential calibration of semidistributed hydrological models. It consists in adding a priori information on optimal parameter sets for each modeling unit of the semidistributed model. Calibration iterations are then performed by jointly maximizing simulation performance and minimizing drifts from the a priori parameter sets. The combination of these two sources of information is handled by a parameter k to which the method is quite sensitive. The method is applied to 1,305 catchments in France over 30 years. The leave-one-out validation shows that, at locations considered as ungauged, model simulations are significantly improved (over all the catchments, the median KGE criterion is increased from 0.75 to 0.83 and the first quartile from 0.35 to 0.66), while model performance at gauged points is not significantly impacted by the use of the regularization approach. Small catchments benefit most from this calibration strategy. These performances are, however, very similar to the performances obtained with a lumped model based on similar conceptualization.},
+  annotation = {\_eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2018WR024266},
+  file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\QNCGSJXB\\Lavenne et al. - 2019 - A Regularization Approach to Improve the Sequentia.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\5DMRIV8Y\\2018WR024266.html},
+  keywords = {regularization,semidistributed model,stepwise calibration,ungauged basins},
+  langid = {english},
+  number = {11},
+  options = {useprefix=true}
+}

vignettes/V01_get_started.Rmd

 ---
 title: "Get Started with airGR"
+author: "Guillaume Thirel, Olivier Delaigue, Laurent Coron"
 bibliography: V00_airgr_ref.bib
 output: rmarkdown::html_vignette
 vignette: >
@@ -10,29 +11,32 @@ vignette: >
 
 # Introduction
 
-**airGR** is a package that brings into the [**R software**](https://cran.r-project.org/) the hydrological modelling tools used and developed at the [Catchment Hydrology Research Group](https://webgr.irstea.fr/en/) at [Irstea (France)](http://www.irstea.fr/en/), including the [**GR rainfall-runoff models**](https://webgr.irstea.fr/en/modeles/) and a snowmelt and accumulation model, [**CemaNeige**](https://webgr.irstea.fr/en/modeles/modele-de-neige/). Each model core is coded in **Fortran** to ensure low computational time. The other package functions (i.e. mainly the calibration algorithm and the efficiency criteria calculation) are coded in **R**. 
+**airGR** is a package that brings into the [**R software**](https://cran.r-project.org/) the hydrological modelling tools used and developed at the [Catchment Hydrology Research Group](https://webgr.inrae.fr/home/) at [INRAE (France)](https://www.inrae.fr/en), including the [**GR rainfall-runoff models**](https://webgr.inrae.fr/models/) that can be applied either on a **lumped** or **semi-distributed** way. A snow accumulation and melt model ([**CemaNeige**](https://webgr.inrae.fr/models/snow-model/)) and the associated functions for the calibration and evaluation of models are also included. Each model core is coded in **Fortran** to ensure low computational time. The other package functions (i.e. mainly the calibration algorithm and the efficiency criteria calculation) are coded in **R**. 
 
 
 The **airGR** package has been designed to fulfill two major requirements: to facilitate the use by non-expert users and to allow flexibility regarding the addition of external criteria, models or calibration algorithms. The names of the functions and their arguments were chosen to this end. **airGR** also contains basics plotting facilities. 
 
 
 
-Six hydrological models and one snowmelt and accumulation model are implemented in **airGR**. The snow model can be used alone or together with the daily hydrological models. 
+Seven hydrological models and one snow melt and accumulation model are implemented in **airGR**. The hydrological models can be applied either on a lumped way or on a semi-distributed way (on sub-catchments). The snow model can either be used alone or with the daily or hourly hydrological models. Naturally each hydrological model can also be used alone.
 
 The models can be called within **airGR** using the following functions: 
 
   * `RunModel_GR4H()`: four-parameter hourly lumped hydrological model [@mathevet_quels_2005] 
+  * `RunModel_GR5H()`: five-parameter hourly lumped hydrological model [@ficchi_adaptive_2017; @ficchi_hydrological_2019] 
   * `RunModel_GR4J()`: four-parameter daily lumped hydrological model [@perrin_improvement_2003] 
   * `RunModel_GR5J()`: five-parameter daily lumped hydrological model [@le_moine_bassin_2008] 
   * `RunModel_GR6J()`: six-parameter daily lumped hydrological model [@pushpalatha_downward_2011] 
   * `RunModel_GR2M()`: two-parameter monthly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_stepwise_2006] 
   * `RunModel_GR1A()`: one-parameter yearly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_linking_2006] 
-  * `RunModel_CemaNeige()`: two-parameter degree-day snowmelt and accumulation model [@valery_as_2014] 
+  * `RunModel_CemaNeige()`: two-parameter degree-day snowmelt and accumulation model [@valery_as_2014; @riboust_revisiting_2019] 
+  * `RunModel_CemaNeigeGR4H()`: combined use of **GR4H** and **CemaNeige** 
+  * `RunModel_CemaNeigeGR5H()`: combined use of **GR5H** and **CemaNeige** 
   * `RunModel_CemaNeigeGR4J()`: combined use of **GR4J** and **CemaNeige** 
   * `RunModel_CemaNeigeGR5J()`: combined use of **GR5J** and **CemaNeige** 
   * `RunModel_CemaNeigeGR6J()`: combined use of **GR6J** and **CemaNeige** 
 
-The [**GRP**](https://webgr.irstea.fr/en/modeles/modele-de-prevision-grp/) forecasting model and the [**Otamin**](https://webgr.irstea.fr/en/modeles/otamin/) predictive uncertainty tool are not available in **airGR**.
+The [**GRP**](https://webgr.inrae.fr/models/hydrological-forecasting-model-grp/) forecasting model and the [**Otamin**](https://webgr.inrae.fr/software/otamin/) predictive uncertainty tool are not available in **airGR**.
 
 In this vignette, we show how to prepare and run a calibration and a simulation with airGR hydrological models.
 
@@ -44,7 +48,7 @@ In the following example, we use a data sample contained in the package. For rea
 
 First, it is necessary to load the **airGR** package:
 
-```{r}
+```{r, warning=FALSE}
 library(airGR)
 ```
 
@@ -59,7 +63,7 @@ Below is presented an example of a `data.frame` of daily hydrometeorological obs
   
 ```{r}
 data(L0123001)
-summary(BasinObs)
+summary(BasinObs, digits = 2)
 ```
 The usual functions (e.g. `read.table()`) can be used to load real-case data sets.
 
@@ -82,7 +86,7 @@ To facilitate the use of the package, there are several functions dedicated to t
 To run a GR hydrological model or CemaNeige, the user has to prepare the input data with the `CreateInputsModel()` function.
 As arguments, this function needs the function name corresponding to the model the user wants to run, a vector of dates, a vector of precipitation and a vector of potential evapotranspiration.
 
-In the example below, we already have the potential evapotranspiration. If the user does not have these data, it is possible to compute it with the [Oudin's formula](http://dx.doi.org/10.1016/j.jhydrol.2004.08.026) with the `PEdaily_Oudin()` function (this function only needs Julian days, daily average air temperature and latitude).
+In the example below, we already have the potential evapotranspiration. If the user does not have these data, it is possible to compute it with the [Oudin's formula](http://dx.doi.org/10.1016/j.jhydrol.2004.08.026) with the `PE_Oudin()` function (this function only needs Julian days, daily average air temperature and latitude).
 
 Missing values (`NA`) of precipitation (or potential evapotranspiration) are **not allowed**.
 
@@ -94,7 +98,6 @@ str(InputsModel)
 ```
 
 
-
 ## RunOptions object
 
 The `CreateRunOptions()` function allows to prepare the options required to the `RunModel*()` functions, which are the actual models functions.
@@ -102,14 +105,14 @@ The `CreateRunOptions()` function allows to prepare the options required to the
 The user must at least define the following arguments:
 
   * `FUN_MOD`: the name of the model function to run
-  * `InputsModel`: the associated inputs data
+  * `InputsModel`: the associated input data
   * `IndPeriod_Run`: the period on which the model is run
   
 To select a period for which the user wants to run the model, select the corresponding indexes for different time periods (not the POSIXt dates), as follows:
 
 ```{r}
-Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y")=="01/01/1990"), 
-               which(format(BasinObs$DatesR, format = "%d/%m/%Y")=="31/12/1999"))
+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"))
 str(Ind_Run)
 ```
 
@@ -136,22 +139,25 @@ The `CreateRunOptions()` function returns warnings if the default initialization
 
 ## InputsCrit object
 
-
 The `CreateInputsCrit()` function allows to prepare the input in order to calculate a criterion. It is possible to define the following arguments:
 
   * `FUN_CRIT`: the name of the error criterion function (the available functions are introduced later on)
   * `InputsModel`: the inputs of the hydrological model previously prepared by the `CreateInputsModel()` function
   * `RunOptions`: the options of the hydrological model previously prepared by the `CreateRunOptions()` function 
-  * `Obs`: the observed discharge expressed in *mm/time step*
+  * `VarObs`: the name of the considered variable (by default `"Q"` for the discharge)
+  * `Obs`: the observed variable time serie (e.g. the discharge expressed in *mm/time step*)
 
 Missing values (`NA`) are **allowed** for observed discharge.
 
+It is possible to compute a composite criterion (e.g. the average between NSE computed on discharge and NSE computed on log of discharge). In this case, users have to provide lists to the following arguments (some of the are optional): `FUN_CRIT`, `Obs`, `VarObs`, `BoolCrit`, `transfo`, `Weights.` 
+
 ```{r}
 InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, 
-                               RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])
+                               RunOptions = RunOptions, VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run])
 str(InputsCrit)
 ```
 
+
 ## CalibOptions object
 
 Before using the automatic calibration tool, the user needs to prepare the calibration options with the `CreateCalibOptions()` function. For that, it is necessary to define the following arguments:
@@ -217,7 +223,6 @@ Performing a model control with **airGR** is similar to running a simulation (se
 
 # Simulation
 
-
 ## Simulation run
 
 To run a model, the user has to use the `RunModel*()` functions (`InputsModel`, `RunOptions` and parameters).
@@ -229,10 +234,9 @@ str(OutputsModel)
 ```
 
 
-
 ## Results preview
 
-Although it is possible for the user to design its own graphics from the outputs of the `RunModel*()` functions, the **airGR** package offers the possibility to make use of the `plot.OutputsModel()` function (or `plot()` with an `OutputsModel` object). This function returns a dashboard of results including various graphs (depending on the model used): 
+Although it is possible for the user to design its own graphics from the outputs of the `RunModel*()` functions, the **airGR** package offers the possibility to make use of the `plot()` function. This function returns a dashboard of results including various graphs (depending on the model used): 
 
   * time series of total precipitation and simulated discharge (and observed discharge if provided)
   * interannual average daily simulated discharge (and daily observed discharge if provided) and interannual average monthly precipitation 
@@ -252,11 +256,19 @@ To evaluate the efficiency of the model, it is possible to use the same criterio
 
 ```{r}
 OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
+str(OutputsCrit)
 ```
 ```{r}
 OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
+str(OutputsCrit)
 ```
 
 
+## Features diagram
+
+The figure below presents the implementation flow of the **airGR** functions.
+<br><br>
+<img src="fig/airGR_features_diagram.svg" alt="airGR features diagram" width=75% />
+
 
 # References

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,30 +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}
+
+```{r local1, warning=FALSE, results='hide', eval=FALSE}
 startGR4J <- c(4.1, 3.9, -0.9, -8.7)
-optPORT <- stats::nlminb(start = startGR4J, 
+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))
@@ -106,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)
 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}
+
+```{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.
@@ -129,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))
@@ -138,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}
-resGLOB <- 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}
+```{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: >
@@ -69,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)
@@ -82,7 +82,7 @@ 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, 
@@ -101,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, Obs = BasinObs$Qmm[Ind_Val])
-
+                                    RunOptions = RunOptions_Val, Obs = BasinObs$Qmm[Ind_Val])
 ```
 
 # Calibration of the GR4J model with the generalist parameter sets
@@ -133,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)
 ```
 
@@ -141,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.
 
@@ -155,17 +155,16 @@ 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)
 
 
@@ -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: >
@@ -11,6 +12,7 @@ vignette: >
 
 ```{r, warning=FALSE, include=FALSE}
 library(airGR)
+load(system.file("vignettesData/vignetteCNHysteresis.rda", package = "airGR"))
 ```
 
 
@@ -36,7 +38,7 @@ We load an example data set from the package. Please note that this data set inc
   
   
 ## loading catchment data
-```{r, warning=FALSE}
+```{r, warning=FALSE, eval=FALSE}
 data(X0310010)
 summary(BasinObs)
 ```
@@ -49,7 +51,7 @@ We assume that the R global environment contains data and functions from the [Ge
 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}
+```{r, warning=FALSE, eval=FALSE}
 ## preparation of the InputsModel object
 InputsModel <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR4J,
                                  DatesR = BasinObs$DatesR, Precip = BasinObs$P,
@@ -77,7 +79,7 @@ Ind_Val <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2005-09-01
 
 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}
+```{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,
@@ -97,7 +99,7 @@ CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
 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}
+```{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,
@@ -124,7 +126,7 @@ InputsCrit_Val  <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
 
 We can now calibrate the model.
 
-```{r, warning=FALSE}
+```{r, warning=FALSE, eval=FALSE}
 ## calibration
 OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
                             InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
@@ -135,7 +137,7 @@ OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_C
 
 Now we can run it on the calibration period and assess it. 
 
-```{r, warning=FALSE, message=FALSE}
+```{r, warning=FALSE, message=FALSE, eval=FALSE}
 ## run on the calibration period
 OutputsModel_Cal <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
                                            RunOptions = RunOptions_Cal, 
@@ -154,7 +156,7 @@ str(OutputsCrit_Cal, max.level = 2)
 
 Now we can run the model on the validation period and assess it. 
 
-```{r, warning=FALSE, message=FALSE}
+```{r, warning=FALSE, message=FALSE, eval=FALSE}
 ## run on the validation period
 OutputsModel_Val <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
                                            RunOptions = RunOptions_Val, 
@@ -176,7 +178,7 @@ str(OutputsCrit_Val, max.level = 2)
 
 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}
+```{r, warning=FALSE, eval=FALSE}
 ## preparation of RunOptions object
 RunOptions_Cal_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
                                           InputsModel = InputsModel,
@@ -201,16 +203,19 @@ CalibOptions_NoHyst <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
 
 We can now calibrate the model.
 
-```{r, warning=FALSE}
+```{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)
+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}
+```{r, warning=FALSE, eval=FALSE}
 OutputsModel_Cal_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
                                                   RunOptions = RunOptions_Cal_NoHyst,  
                                                   Param = OutputsCalib_NoHyst$ParamFinalR)
@@ -223,7 +228,7 @@ OutputsModel_Val_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
 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}
+```{r, warning=FALSE, message=FALSE, eval=FALSE}
 OutputsCrit_Cal_NoHyst <- ErrorCrit(InputsCrit = InputsCrit_Cal,
                                     OutputsModel = OutputsModel_Cal_NoHyst)
 

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