diff --git a/vignettes/v00_teaching_hydrology.rmd b/vignettes/v00_teaching_hydrology.rmd
index 70d13262488b9452eff96abdea6e5ca8af5f0152..5137a5a2df9ef66b7fc11d24e5b660f10b77da75 100644
--- a/vignettes/v00_teaching_hydrology.rmd
+++ b/vignettes/v00_teaching_hydrology.rmd
@@ -83,7 +83,7 @@ Rainfall-runoff models are composed of different elements, e.g. conceptual reser
 
 To illustrate the production and the routing parts of hydrological modelling which are present in any model, it is possible to use the different GR models available using `r airGRteaching` and to produce rainfall-runoff transformations considering different model parameter values.
 
-The GR4J model [@perrin_improvement_2003] comprises a production store (X1 parameter, the maximal capacity of the production store [mm]), which determines the actual evapotranspiration and the net rainfall. The routing of net rainfall is determined through two unit hydrographs (X4 parameter, the unit hydrograph time constant [days]) and a routing store (X3 parameter, the one-day ahead maximal capacity [mm]). A final component represents the intercatchment groundwater exchange (X2 parameter, the exchange coefficient [mm/day]).
+The GR4J model [@perrin_improvement_2003] comprises a production store (X1 parameter, the maximal capacity of the production store [mm]), which determines the actual evapotranspiration and the net rainfall. The routing of net rainfall is determined through two unit hydrographs (X4 parameter, the unit hydrograph time constant [d]) and a routing store (X3 parameter, the one-day ahead maximal capacity [mm]). A final component represents the intercatchment groundwater exchange (X2 parameter, the exchange coefficient [mm/d]).
 
 The following command lines allow to prepare the hydro-meteorological data for use with GR4J.
 
@@ -111,7 +111,7 @@ per_sim <- range(prep_no_q$InputsModel$DatesR)
 As an example, the following commands and figure illustrate the role of the X2 parameter of GR4J on the production part of the rainfall-runoff transformation, by testing different X2 values. We observe higher streamflow values simulated with higher X2 values. 
 
 ```{r, v00_fig_gr4j_x2, echo=TRUE, eval=TRUE, warning=FALSE, fig.width=3*3, fig.height=3*1.7, out.width='98%'}
-# Different X2 values around its median values (0 [mm/day])
+# Different X2 values around its median values (0 [mm/d])
 param_x2 <- seq(from = -2, to = 2, by = 1)
 
 # Combination of parameter values (X1, X3 and X4 are fixed; X2 changes)
@@ -134,19 +134,19 @@ ind_zoom <- 400:430
 col_param_x2 <- colorRampPalette(c("green1", "green4"))(ncol(sim_x2))
 matplot(x = as.POSIXct(prep_no_q$InputsModel$DatesR[ind_zoom]),
         y = sim_x2[ind_zoom, ],
-        xlab = "time [days]", ylab = "flow [mm/day]",
+        xlab = "time [d]", ylab = "flow [mm/d]",
         type = "l", lty = 1, lwd = 2, col = col_param_x2)
 legend("topright",
        legend = sprintf("% .1f", param_x2),
        lwd = 2, col = col_param_x2,
-       title = "X2 values [mm/day]",
+       title = "X2 values [mm/d]",
        bg = "white")
 ```
 
 Moreover, the following commands and figure illustrate the role of the X4 parameter on the routing part of the rainfall-runoff transformation, with delayed flood peak values when considering higher X4 values.
 
 ```{r, v00_fig_gr4j_x4, echo=TRUE, eval=TRUE, warning=FALSE, fig.width=3*3, fig.height=3*1.7, out.width='98%'}
-# Different X4 values around its median values (1.4 [days])
+# Different X4 values around its median values (1.4 [d])
 param_x4 <- seq(from = 1.0, to = 3.0, by = 0.5)
 
 # Combination of parameter values (X1, X2 and X3 are fixed; X4 changes)
@@ -169,12 +169,12 @@ ind_zoom <- 400:430
 col_param_x4 <- colorRampPalette(c("steelblue1", "steelblue4"))(ncol(sim_x4))
 matplot(x = as.POSIXct(prep_no_q$InputsModel$DatesR[ind_zoom]),
         y = sim_x4[ind_zoom, ],
-        xlab = "time [days]", ylab = "flow [mm/day]",
+        xlab = "time [d]", ylab = "flow [mm/d]",
         type = "l", lty = 1, lwd = 2, col = col_param_x4)
 legend("topright",
        legend = sprintf("% .1f", param_x4),
        lwd = 2,col = col_param_x4, 
-       title = "X4 values [days]",
+       title = "X4 values [d]",
        bg = "white")
 ```
 
@@ -196,7 +196,7 @@ sim_x4_y <- SeriesAggreg(x = sim_x4_y,
 # Graphical comparison
 matplot(x = sim_x2_y$DatesR, y = sim_x2_y[, -1], 
         type = "l", lty = 1, lwd = 2, col = col_param_x2, 
-        xlab = "time [years]", ylab = "flow [mm/year]")
+        xlab = "time [yr]", ylab = "flow [mm/yr]")
 matlines(x = sim_x4_y$DatesR, y = sim_x4_y[, -1], 
          type = "l", lty = 1, lwd = 2, col = col_param_x4)
 legend("topright", 
@@ -393,7 +393,7 @@ tab_cal_reg <- SeriesAggreg(tab_cal,
 col_snow <- c("black", rep("orangered", 2))
 lty_snow <- c(1, 1:2)
 matplot(y = tab_cal_reg[, grep("^Q", colnames(tab_cal))],
-        xlab = "time [months]", ylab = "flow [mm/day]", 
+        xlab = "time [months]", ylab = "flow [mm/d]", 
         type = "l", lty = lty_snow, lwd = 2, col = col_snow)
 legend("topright", 
        legend = c("Qobs", "Qsim without snow mod.", "Qsim with snow mod."),
@@ -459,7 +459,7 @@ tab_sim_reg <- SeriesAggreg(tab_sim_trsf,
 col_trsf <- c("black", rep("orangered", 3))
 lty_trsf <- c(1, 1:3)
 matplot(y = tab_sim_reg[, -1],
-        xlab = "time [months]", ylab = "flow [mm/day]", 
+        xlab = "time [months]", ylab = "flow [mm/d]", 
         type = "l", lty = lty_trsf, lwd = 2, col = col_trsf)
 legend("bottomleft", 
        legend = c("Qobs", "Qsim", "sqrt(Qsim)", "log(Qsim)"),
@@ -497,7 +497,7 @@ tab_crit <- data.frame(Date     = as.POSIXct(cal_nse$OutputsModel$DatesR),
 col_crit <- c("black", rep("orangered", 2))
 lty_crit <- c(1, 1:2)
 matplot(x = tab_crit$Date, y = tab_crit[, -1],
-        xlab = "time [days]", ylab = "flow [mm/day]", 
+        xlab = "time [d]", ylab = "flow [mm/d]", 
         type = "l", lty = lty_crit, lwd = 2, col = col_crit,
         xlim = as.POSIXct(x = c("2004-01-01", "2004-03-01"), tz = "UTC"))
 legend("topleft", 
@@ -587,7 +587,7 @@ ts_obs_y$Arid <- ts_obs_y$Evap / ts_obs_y$Ptot
 # Identification of wetter and dryer hydrological years
 barplot(height = ts_obs_y$Arid, 
         names.arg = format(ts_obs_y$Date, format = "%Y"), 
-        xlab = "time [years]", ylab = "aridity index [-]",
+        xlab = "time [yr]", ylab = "aridity index [-]",
         col = "royalblue")
 
 # Wet and dry periods
diff --git a/vignettes/v01_EN_flow_reconstruction.rmd b/vignettes/v01_EN_flow_reconstruction.rmd
index fcb35f72b29fd43a5f8bc6e98edd8517ead25cb2..5bd796595d6326d9252b0d2be8bff2c2674b1ddb 100644
--- a/vignettes/v01_EN_flow_reconstruction.rmd
+++ b/vignettes/v01_EN_flow_reconstruction.rmd
@@ -206,9 +206,9 @@ The final aims of this exercise is to reconstruct the streamflows of the catchme
 
 The data set available to the rainfall-runoff modelling consists of:
 
-* a time series of daily total precipitation (liquid + solid) [mm/day] (`Ptot`);
-* a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/day] (`Evap`);
-* a time series of daily streamflows expressed as a specific discharge [mm/day] (`Qmmd`).
+* a time series of daily total precipitation (liquid + solid) [mm/d] (`Ptot`);
+* a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/d] (`Evap`);
+* a time series of daily streamflows expressed as a specific discharge [mm/d] (`Qmmd`).
 
 The daily time series can be aggregated to a monthly time step using the `SeriesAggreg()` function.
 
diff --git a/vignettes/v02_EN_flow_forecasting.rmd b/vignettes/v02_EN_flow_forecasting.rmd
index eedc457c7654fd0fd1303a1d6c11c2f2fe2f31ec..255492e26c57769422ff5cf953cbc051f9aa450e 100644
--- a/vignettes/v02_EN_flow_forecasting.rmd
+++ b/vignettes/v02_EN_flow_forecasting.rmd
@@ -169,7 +169,7 @@ ind_xlim <- c(nrow(ts_plot)-300, nrow(ts_plot))
 
 # Plotting hyrdograph
 plot(x = ts_plot$Date, y = ts_plot$Qmmd,
-     xlab = "time [days]", ylab = "flow [mm/day]",
+     xlab = "time [d]", ylab = "flow [mm/d]",
      main = paste("Year", year_max),
      type = "l", lwd = 1,
      log = "y",
@@ -259,9 +259,9 @@ This last step aims at performing several rainfall-runoff simulations considerin
 
 The data set available to the rainfall-runoff modelling consists of:
 
-* a time series of daily total precipitation (liquid + solid) [mm/day] (`Ptot`);
-* a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/day] (`Evap`);
-* a time series of daily streamflows expressed as a specific discharge [mm/day] (`Qmmd`).
+* a time series of daily total precipitation (liquid + solid) [mm/d] (`Ptot`);
+* a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/d] (`Evap`);
+* a time series of daily streamflows expressed as a specific discharge [mm/d] (`Qmmd`).
 
 
 
@@ -379,9 +379,9 @@ plot(cal_hist)
 The six parameters and the value of the calibration criterion ($NSE_{ln}$) obtained after the automatic calibration procedure are:
 
 * X1 = `r round(param_cal_hist[1], digits = 0)` [mm]
-* X2 = `r round(param_cal_hist[2], digits = 3)` [mm/day]
+* X2 = `r round(param_cal_hist[2], digits = 3)` [mm/d]
 * X3 = `r round(param_cal_hist[3], digits = 0)` [mm]
-* X4 = `r round(param_cal_hist[4], digits = 3)` [day]
+* X4 = `r round(param_cal_hist[4], digits = 3)` [d]
 * X5 = `r round(param_cal_hist[5], digits = 3)` [-]
 * X6 = `r round(param_cal_hist[6], digits = 2)` [mm]
 * $NSE_{ln}$ = `r round(crit_cal_hist, digits = 3)` [-]
diff --git a/vignettes/v02_FR_flow_forecasting.rmd b/vignettes/v02_FR_flow_forecasting.rmd
index 798f9e66d38aa9622c9f11dbd09d388431e6ec6e..a56587731c01a1ee80f1cbbb850b378c7df7dd83 100644
--- a/vignettes/v02_FR_flow_forecasting.rmd
+++ b/vignettes/v02_FR_flow_forecasting.rmd
@@ -169,7 +169,7 @@ ind_xlim <- c(nrow(ts_plot)-300, nrow(ts_plot))
 
 # Plotting hyrdograph
 plot(x = ts_plot$Date, y = ts_plot$Qmmd,
-     xlab = "time [days]", ylab = "flow [mm/day]",
+     xlab = "time [d]", ylab = "flow [mm/d]",
      main = paste("Year", year_max),
      type = "l", lwd = 1,
      log = "y",
diff --git a/vignettes/v03_EN_impact_CC_flow.rmd b/vignettes/v03_EN_impact_CC_flow.rmd
index 1004d632e30b98980eb3af011b8bdeb84577c39b..7cf314a6498be03c39b1cff561db5edfb4e76891 100644
--- a/vignettes/v03_EN_impact_CC_flow.rmd
+++ b/vignettes/v03_EN_impact_CC_flow.rmd
@@ -161,13 +161,13 @@ par(mfrow = c(1, 2))
 
 # Present time regime
 plot(x = reg_d$Date, y = reg_d$Qmmd, 
-     xlab = "time (days)", ylab = "flow [mm/day]", 
+     xlab = "time [d]", ylab = "flow [mm/d]", 
      main = "Present time climate",
      type = "l")
 
 # Future time regime
 plot(x = reg_d$Date, y = reg_d$Qmmd, 
-     xlab = "time (days)", ylab = "", 
+     xlab = "time [d]", ylab = "", 
      main = "Future climate",
      type = "n", yaxt = "n")
 axis(side = 2, labels = FALSE)
@@ -239,10 +239,10 @@ The parameters obtained after calibration will then be used to perform rainfall-
 
 The data set available to the rainfall-runoff modelling consists of:
 
-* a time series of daily total precipitation (liquid + solid) [mm/day] (`Ptot`);
+* a time series of daily total precipitation (liquid + solid) [mm/d] (`Ptot`);
 * a time series of daily air temperatures [°C] (`Temp`);
-* a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/day] (`Evap`);
-* a time series of daily streamflows expressed as a specific discharge [mm/day] (`Qmmd`);
+* a time series of daily potential evapotranspiration calculated with the @oudin_which_2005 formula [mm/d] (`Evap`);
+* a time series of daily streamflows expressed as a specific discharge [mm/d] (`Qmmd`);
 * a distribution of the altitude on the catchment area (hypsometric curve) [m] (`Hypso`).
 
 You also have the tables of average monthly changes estimated by your fellow climatologists presented in the introduction.
@@ -328,7 +328,7 @@ The regimes calculated in this way can be represented graphically (see following
 ```{r, v03_fig_regime, echo=FALSE, eval=TRUE, fig.width=3*3, fig.height=3*1.7, out.width='98%'}
 # Graphical parameters
 plot(x = reg_hist_m$Dates, y = reg_hist_m$Qobs,
-     xlab = "time (months)", ylab = "flow [mm/month]",
+     xlab = "time [months]", ylab = "flow [mm/month]",
      type = "l")
 lines(x = reg_hist_m$Dates, y = reg_hist_m$Qsim,
       type = "l", col = "orangered")
@@ -480,7 +480,7 @@ The evolutions are more marked at the daily time step (see following figure).
 ```{r, v03_fig_regime_day_cc, echo=FALSE, eval=TRUE, fig.width=3*3, fig.height=3*1.7, out.width='98%'}
 # Plot framework
 matplot(x = reg_cc_d$Dates, y = reg_cc_d[, -1], 
-        xlab = "time (days)", ylab = "flow [mm/day]", 
+        xlab = "time [d]", ylab = "flow [mm/d]", 
         type = "l", lty = 1, 
         col = c("black", "orangered", col_qscen))
 
diff --git a/vignettes/v03_FR_impact_CC_flow.rmd b/vignettes/v03_FR_impact_CC_flow.rmd
index 25bb7e06493cd8c0c853fba803007da3109763f4..ed4410c1e9ff0fc79bbc6b8ef707a9bd14f50dc7 100644
--- a/vignettes/v03_FR_impact_CC_flow.rmd
+++ b/vignettes/v03_FR_impact_CC_flow.rmd
@@ -160,13 +160,13 @@ par(mfrow = c(1, 2))
 
 # Present time regime
 plot(x = reg_d$Date, y = reg_d$Qmmd, 
-     xlab = "time (days)", ylab = "flow [mm/day]", 
+     xlab = "time [d]", ylab = "flow [mm/d]", 
      main = "Present time climate",
      type = "l")
 
 # Future time regime
 plot(x = reg_d$Date, y = reg_d$Qmmd, 
-     xlab = "time (days)", ylab = "", 
+     xlab = "time [d]", ylab = "", 
      main = "Future climate",
      type = "n", yaxt = "n")
 axis(side = 2, labels = FALSE)
@@ -328,7 +328,7 @@ Les régimes ainsi calculés peuvent être représentés graphiquement (cf. figu
 ```{r, v03_fig_regime, echo=FALSE, eval=TRUE, fig.width=3*3, fig.height=3*1.7, out.width='98%'}
 # Graphical parameters
 plot(x = reg_hist_m$Dates, y = reg_hist_m$Qobs,
-     xlab = "time (months)", ylab = "flow [mm/month]",
+     xlab = "time [months]", ylab = "flow [mm/month]",
      type = "l")
 lines(x = reg_hist_m$Dates, y = reg_hist_m$Qsim,
       type = "l", col = "orangered")
@@ -480,7 +480,7 @@ Les évolutions sont plus marquées au pas de temps journalier (cf. figure suiva
 ```{r, v03_fig_regime_day_cc, echo=FALSE, eval=TRUE, fig.width=3*3, fig.height=3*1.7, out.width='98%'}
 # Plot framework
 matplot(x = reg_cc_d$Dates, y = reg_cc_d[, -1], 
-        xlab = "time (days)", ylab = "flow [mm/day]", 
+        xlab = "time [d]", ylab = "flow [mm/d]", 
         type = "l", lty = 1, 
         col = c("black", "orangered", col_qscen))