diff --git a/NAMESPACE b/NAMESPACE
index ac8b38d5499127f4e777a5e8a7e08825c992d425..8c53d31c6736ba5e9110cd0d1799a12a870d4e0d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -44,6 +44,7 @@ export(max_na_rm)
 export(mean_na_rm)
 export(min_na_rm)
 export(not_na_count)
+export(pmedian)
 export(readDrias2020)
 export(readQsim)
 export(saveFlowDB)
diff --git a/R/calcQA.R b/R/calcQA.R
index b65ab29b705b99be7bead647c714efc944accd22..64eb6cc7a5b5a2a0a26c00b6db9e1c8237241e52 100644
--- a/R/calcQA.R
+++ b/R/calcQA.R
@@ -15,7 +15,6 @@
 #' @source Chazot, Sébastien, Perrin, Charles, Jean-Philippe Vidal, Eric Sauquet, Mathilde Chauveau, et Nathalie Rouchy. 2012. « Explore 2070 - Lot Hydrologie de surface - A2 - Résultats : Fiches, Cartes et Graphes Â». Ministère de l’écologie, du développement durable, des transports et du logement.
 
 #'
-#' @examples
 calcQA <- function(x, threshold = 0.8, ...) {
     UseMethod("calcQA", x)
 }
diff --git a/R/pmedian.R b/R/pmedian.R
new file mode 100644
index 0000000000000000000000000000000000000000..930dc667b9fce8754a895a2375dfe945bb4f8b88
--- /dev/null
+++ b/R/pmedian.R
@@ -0,0 +1,26 @@
+#' Median value
+#'
+#' @description
+#' Returns the (regular or parallel) median of the input values.
+#'
+#' `pmedian`  takes one or more vectors as arguments which should have same dimensions.
+#'
+#' @param ...  Several [matrix] of same dimensions
+#'
+#' @return a [matrix] of same dimension as the ones in `...` with median value calculated for each cell.
+#' @export
+#'
+pmedian <- function(...) {
+    l <- list(...)
+    m <- sapply(seq.int(length(l[[1]])), function(i) {
+        cells <- sapply(l, function(x) x[i])
+        median(cells)
+    })
+    if (is.matrix(l[[1]])) {
+        m <- matrix(m, ncol = ncol(l[[1]]))
+        colnames(m) <- colnames(l[[1]])
+    } else {
+        names(m) <- names(l[[1]])
+    }
+    return(m)
+}
diff --git a/bookdown/06-evolution_indicateurs_hydrologiques.Rmd b/bookdown/06-evolution_indicateurs_hydrologiques.Rmd
index c50585ee61ddc945423ed078190d6f98abd6a8d1..12b075457f294a53bdf4d0f35a8a1118df648da2 100644
--- a/bookdown/06-evolution_indicateurs_hydrologiques.Rmd
+++ b/bookdown/06-evolution_indicateurs_hydrologiques.Rmd
@@ -23,7 +23,7 @@ BasinsObs <- loadBasinsObs(file, cfg = cfg)
 On peut charger toutes les données climatiques de tous les scénarios GCM/RCM dans une liste avec le script suivant&nbsp;:
 
 ```{r}
-rcps <- cfg$hydroclim$drias$rcp[-1]
+rcps <- "rcp4.5"
 scenarios <- gsub("/", "_", cfg$hydroclim$drias$scenarios)
 scenarios <- rep(scenarios, length(rcps))
 
@@ -34,6 +34,19 @@ loadAllBasinsObs <- function(rcp, scenario) {
 }
 AllBasinsObs <- mapply(loadAllBasinsObs, rcp = rcps, scenario = scenarios, SIMPLIFY = FALSE)
 names(AllBasinsObs) <- paste(rcps, scenarios, sep = " - ")
+
+```
+
+On réarrange les données pour avoir une liste de data.frame pour chaque variable étudiée:
+
+```{r}
+formatObs <- function(BasinObs, item) {
+    cbind(DatesR = BasinObs$DatesR, as.data.frame(BasinObs[[item]]))
+}
+
+P_drias <- lapply(AllBasinsObs, formatObs, item = "P")
+E_drias <- lapply(AllBasinsObs, formatObs, item = "E")
+T_drias <- lapply(AllBasinsObs, formatObs, item = "Temp")
 ```
 Chaque élément de la liste `AllBasinsObs` est une liste structurée de classe *BasinsObs* (Voir `?createBasinsObs`).
 
@@ -74,7 +87,25 @@ Pour les débits&nbsp;:
 - Débit moyen de période de retour 5 ans
 
 ```{r}
+monthlyMean <- function(df) {
+    dfOut <- SeriesAggreg(df, "%m", rep("mean", ncol(df) - 1))
+    cbind(month = lubridate::month(dfOut$DatesR), as.matrix(dfOut[, -1]))
+}
 
+selectPeriod <- function(df, period) {
+    mPeriods <- data.frame(period = c("ref", "proche", "milieu", "fin"),
+                           start = as.POSIXct(c("1976-01-01", "2021-01-01", "2041-01-01", "2071-01-01"), tz = "UTC"),
+                           end = as.POSIXct(c("2005-12-31", "2050-12-31", "2070-12-31", "2100-12-31"), tz = "UTC"))
+    period <- period
+    if (!period %in% mPeriods$period) {
+        stop("`period`should be equal to ", paste(sprintf("\"%s\"", mPeriods$period), collapse = ", "))
+    }
+    selectDates <- df$DatesR >= mPeriods$start[mPeriods$period == period] & 
+                   df$DatesR <= mPeriods$end[mPeriods$period == period]
+    return(df[selectDates, ])
+}
+
+monthlyMean(selectPeriod(P_drias[[1]][, 1:3], "ref"))
 ```
 
 
@@ -117,11 +148,66 @@ Les données utilisées en entrées sont&nbsp;:
 Les tableaux de synthèse fournissent la valeur minimale, médiane et maximale des évolutions parmi les couples scénario/modèle climatiques.
 
 ```{r}
-calcComparison <- function(obs, sim, scenarios, FUN, ...) {
-    
+calcDelta <- function(sim, FUN, period, delta = "*", notFirstCol = FALSE, ...) {
+    if (!delta %in% c("*", "+")) stop("`delta`should be equal to \"*\" or \"+\"")
+    simRef <- lapply(sim, selectPeriod, period = "ref")
+    vRef <- lapply(simRef, FUN, ...)
+    simFut <- lapply(sim, selectPeriod, period = period)
+    vFut <- lapply(simFut, FUN, ...)
+    deltas <- lapply(names(sim), function(scenario) {
+        if (delta == "*") {
+            m <- (vFut[[scenario]] - vRef[[scenario]]) / vRef[[scenario]] * 100
+        } else {
+            m <- vFut[[scenario]] - vRef[[scenario]]
+        }
+        if (notFirstCol) {
+            m[, 1] <- vRef[[scenario]][, 1]
+        }
+        return(m)
+    })
+    list(
+        min = do.call(pmin, deltas),
+        med = do.call(pmedian, deltas),
+        max = do.call(pmax, deltas)
+    )
 }
 ```
 
+Exemple pour les pluies moyennes mensuelles:
+
+```{r}
+deltaP <- calcDelta(P_drias, monthlyMean, period = "fin", notFirstCol = TRUE)
+lapply(deltaP, function(x) x[, 1:3])
+```
+
+Exemple pour la température moyenne mensuelle:
+
+```{r}
+deltaT <- calcDelta(T_drias, monthlyMean, period = "fin", delta = "+", notFirstCol = TRUE)
+lapply(deltaT, function(x) x[, 1:3])
+```
+
+**N.B.:** les données climatiques correspondent ici au données moyennes du sous-bassin versant et pas le bassin versant entier. Il faudrait agréger les données des bassins amont pour avoir une moyenne du bassin pour chaque station.
+
+Exemple sur le QMNA5:
+
+```{r}
+deltaQMNA5 <- calcDelta(QsimDrias, calcQMNAn, period = "fin", return_period = 5)
+lapply(deltaQMNA5, function(x) x[1:2])
+```
+
+Exemple sur le VCN10-2:
+
+```{r}
+deltaVCN10_2 <- calcDelta(QsimDrias, calcVCNn, period = "fin", k = 10, return_period = 2)
+lapply(deltaVCN10_2, function(x) x[1:2])
+```
+
+Exemple sur le QJXA10:
+
+```{r}
+deltaQJXA10 <- calcDelta(QsimDrias, calcQJXAn, period = "fin", return_period = 2)
+lapply(deltaQJXA10, function(x) x[1:2])
+```
 
-## Lecture des données
 
diff --git a/man/pmedian.Rd b/man/pmedian.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..0839fe62c232aad4d247286e801d7efbfd24f1f0
--- /dev/null
+++ b/man/pmedian.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pmedian.R
+\name{pmedian}
+\alias{pmedian}
+\title{Median value}
+\usage{
+pmedian(...)
+}
+\arguments{
+\item{...}{Several \link{matrix} of same dimensions}
+}
+\value{
+a \link{matrix} of same dimension as the ones in \code{...} with median value calculated for each cell.
+}
+\description{
+Returns the (regular or parallel) median of the input values.
+
+\code{pmedian}  takes one or more vectors as arguments which should have same dimensions.
+}