diff --git a/DESCRIPTION b/DESCRIPTION
index d8945af7ee56f0c1db0b68f63d4a4a6bc005a894..c5cb0158580669be56e81f7f5aff425babb2c944 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -29,6 +29,7 @@ Imports:
 Suggests:
   knitr, rmarkdown,
   caRamel, coda, DEoptim, dplyr, FME, ggmcmc, hydroPSO, imputeTS, Rmalschains,
+  ggplot2, GGally,
   testthat
 Description: Hydrological modelling tools developed at INRAE-Antony (HYCAR Research Unit, France). The package includes several conceptual rainfall-runoff models (GR4H, GR5H, GR4J, GR5J, GR6J, GR2M, GR1A), a snow accumulation and melt model (CemaNeige) and the associated functions for their calibration and evaluation. Use help(airGR) for package description and references.
 License: GPL-2
diff --git a/NEWS.md b/NEWS.md
index 13d0bb6ee6644fa4c44acb8efd26b448790015b0..242e55295f85b8464c550877bcbd272003a68286 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,10 @@
 
 ### 1.6.9.36 Release Notes (2021-01-28)
 
+#### New features
+
+- Added a section 'param_optim' vignette to explain how to manage with multiobjective optimization using the 'CaRamel' package. ([#61](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/61))
+
 #### Major user-visible changes
 
 - `Imax()` now returns an error message when `IndPeriod_Run` doesn't select 24 hours by day, instead of `numeric(0)`. ([#92](https://gitlab.irstea.fr/HYCAR-Hydro/airgr/-/issues/92))
diff --git a/inst/vignettesData/vignetteParamOptimCaramel.rda b/inst/vignettesData/vignetteParamOptimCaramel.rda
new file mode 100644
index 0000000000000000000000000000000000000000..7d2e9fcc41220cab22fa5b1e57e7c666ef74803e
Binary files /dev/null and b/inst/vignettesData/vignetteParamOptimCaramel.rda differ
diff --git a/vignettes/V02.1_param_optim.Rmd b/vignettes/V02.1_param_optim.Rmd
index 462b7d4d73b9885545835d8ca7bd74cfcaf2c04b..2dd121e92d89e6838e3c7db3af73322e74080400 100644
--- a/vignettes/V02.1_param_optim.Rmd
+++ b/vignettes/V02.1_param_optim.Rmd
@@ -1,5 +1,6 @@
 ---
 title: "Plugging in new calibration algorithms in airGR"
+author: "François Bourgin, Guillaume Thirel"
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteEngine{knitr::rmarkdown}
@@ -20,6 +21,7 @@ library(GGally)
 # source("airGR.R")
 set.seed(321)
 load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
 ```
 
 
@@ -160,7 +162,7 @@ As it can be seen in the table below, the four additional optimization strategie
 ```{r, 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"),
@@ -189,9 +191,9 @@ InputsCrit_inv <- InputsCrit
 InputsCrit_inv$transfo <- "inv"
   
 MOptimGR4J <- function(i) {
-  
-  if (algo == "caRamel") ParamOptim = x[i, ]
-  
+  if (algo == "caRamel") {
+    ParamOptim <- x[i, ]
+  }
   ## Transformation of the parameter set to real space
   RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
                                             Direction = "TR")
@@ -224,7 +226,7 @@ optMO <- caRamel::caRamel(nobj = 2,
                           popsize = 100,
                           archsize = 100,
                           maxrun = 15000,
-                          prec = rep(1.e-3,2),
+                          prec = rep(1.e-3, 2),
                           carallel = FALSE,
                           graph = FALSE)
 ```
@@ -234,34 +236,39 @@ The algorithm returns parameter sets that describe the pareto front, illustratin
 ```{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_inv") +
+  coord_equal(xlim = c(0.4, 0.9), ylim = c(0.4, 0.9)) + 
+  xlab("KGE") + ylab("KGE [1/Q]") +
   theme_bw()
 ```
 
 The pameter 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, 1, FUN = function(x) {
-  airGR::TransfoParam_GR4J(x, Direction = "TR")})
-ggpairs(data.frame(t(param_optMO)), diag = NULL) + theme_bw()
+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, 1, FUN = function(x) {
+RunOptions$Outputs_Sim <- "Qsim"
+run_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
   airGR::RunModel_GR4J(InputsModel = InputsModel,
                        RunOptions = RunOptions,
-                       x)}$Qsim)
-run_optMO = data.frame(run_optMO)
+                       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") +
+  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_classic()
+  theme_bw()
 ```