diff --git a/R/SeriesAggreg.InputsModel.R b/R/SeriesAggreg.InputsModel.R
index a6414fb823235d3fb19904602adcd8a177b82466..37513f46d1a15322f5e9955c682a41b33904f8ba 100644
--- a/R/SeriesAggreg.InputsModel.R
+++ b/R/SeriesAggreg.InputsModel.R
@@ -1,6 +1,7 @@
-SeriesAggreg.InputsModel <- function(x, ...) {
+SeriesAggreg.InputsModel <- function(x, Format, ...) {
   SeriesAggreg.list(x,
-                    ConvertFun = .GetAggregConvertFun(names(x)),
+                    Format,
+                    ConvertFun = NA,
                     except = c("ZLayers", "LengthHydro", "BasinAreas"),
                     ...)
 }
diff --git a/R/SeriesAggreg.OutputsModel.R b/R/SeriesAggreg.OutputsModel.R
index 264bd5b7fb7f41e226004ce244d545246eefe881..3bf0d8aefd24f930c4d6c0ba36fd9786d11e640f 100644
--- a/R/SeriesAggreg.OutputsModel.R
+++ b/R/SeriesAggreg.OutputsModel.R
@@ -1,6 +1,7 @@
-SeriesAggreg.OutputsModel <- function(x, ...) {
+SeriesAggreg.OutputsModel <- function(x, Format, ...) {
   SeriesAggreg.list(x,
-                    ConvertFun = .GetAggregConvertFun(names(x)),
+                    Format,
+                    ConvertFun = NA,
                     except = "StateEnd",
                     ...)
 }
diff --git a/R/SeriesAggreg.data.frame.R b/R/SeriesAggreg.data.frame.R
index 1761ba9ac234980e3fd5f58b16860174c53aefe9..f963a0540325f30772829a6ad1c67a0aa9fbcd9d 100644
--- a/R/SeriesAggreg.data.frame.R
+++ b/R/SeriesAggreg.data.frame.R
@@ -48,14 +48,27 @@ SeriesAggreg.data.frame <- function(x,
   Format <- match.arg(Format, choices = listFormat)
 
   ## check ConvertFun
-  listConvertFun <- list(sum = sum, mean = mean)
-  ConvertFun <- names(listConvertFun)[match(ConvertFun, names(listConvertFun))]
-  if (anyNA(ConvertFun)) {
-    stop("'ConvertFun' should be a one of 'sum' or 'mean'")
-  }
   if (length(ConvertFun) != (ncol(x) - 1)) {
     stop(sprintf("'ConvertFun' must be of length %i (ncol(x)-1)", ncol(x) - 1))
   }
+  listConvertFun <- lapply(unique(ConvertFun), function(y) {
+    if (!grepl("^q\\d+$", y, ignore.case = TRUE)) {
+      return(match.fun(y))
+    }
+  })
+  names(listConvertFun) <- unique(ConvertFun)
+  lapply(ConvertFun, function(y) {
+    if (!grepl("^q\\d+$", y, ignore.case = TRUE)) {
+      TestOutput <- listConvertFun[[y]](1:10)
+      if(!is.numeric(TestOutput)) {
+        stop(sprintf("Returned value of '%s' function should be numeric", y))
+      }
+      if(length(TestOutput) != 1) {
+        stop(sprintf("Returned value of '%s' function should be of length 1", y))
+      }
+    }
+  })
+
   ## check YearFirstMonth
   msgYearFirstMonth <- "'YearFirstMonth' should be a single vector of numeric value between 1 and 12"
   YearFirstMonth <- match(YearFirstMonth, 1:12)
@@ -135,16 +148,28 @@ SeriesAggreg.data.frame <- function(x,
     }
     TabSeries2$Fac2 <- TabSeries2$Selec2
     TabSeries2$Selec <- !duplicated(TabSeries2$Selec2)
-    ConvertFun <- rep("mean", ncol(x) - 1)
   }
-
   listTsAggreg <- lapply(names(listConvertFun), function(y) {
     if (any(ConvertFun == y)) {
       colTsAggreg <- c("Fac2", colnames(x)[-1L][ConvertFun == y])
-      aggregate(. ~ Fac2,
-                data = TabSeries2[, colTsAggreg],
-                FUN = listConvertFun[[y]],
-                na.action = na.pass)
+      if (grepl("^q\\d+$", y, ignore.case = TRUE)) {
+        probs <- as.numeric(gsub("^q", "", y, ignore.case = TRUE)) / 100
+        if (probs < 0 || probs > 1) {
+          stop("'Q...' format of argument 'ConvertFun' must be an integer between 0 and 100")
+        }
+        aggregate(. ~ Fac2,
+                  data = TabSeries2[, colTsAggreg],
+                  FUN = quantile,
+                  na.action = na.pass,
+                  probs = probs,
+                  type = 8,
+                  na.rm = TRUE)
+      } else {
+        aggregate(. ~ Fac2,
+                  data = TabSeries2[, colTsAggreg],
+                  FUN = listConvertFun[[y]],
+                  na.action = na.pass)
+      }
     } else {
       NULL
     }
diff --git a/R/SeriesAggreg.list.R b/R/SeriesAggreg.list.R
index ab7df82d2cfe03ed205d8a47c31231077bb11be8..2298a9f136ab8d6117db83d100f5e17d295b9ab6 100644
--- a/R/SeriesAggreg.list.R
+++ b/R/SeriesAggreg.list.R
@@ -16,6 +16,16 @@ SeriesAggreg.list <- function(x,
     warning("deprecated 'NewTimeFormat' argument: 'Format' argument is used instead",
             call. = FALSE)
   }
+  # Check ConvertFun
+  if (any(class(x) %in% c("InputsModel", "OutputsModel"))) {
+    if (!all(is.na(ConvertFun))) {
+      warning("Argument 'ConvertFun' is ignored on 'InputsModel' and 'OutputsModel' objects")
+    }
+  } else if (length(ConvertFun)!=1) {
+    stop("Argument 'ConvertFun' must be of length 1 with 'list' object")
+  } else if (!is.character(ConvertFun)) {
+    stop("Argument 'ConvertFun' must be a character")
+  }
 
   # Determination of DatesR
   if (!is.null(x$DatesR)) {
@@ -63,13 +73,11 @@ SeriesAggreg.list <- function(x,
   }
   dfOut <- NULL
   if (length(cols)) {
-    ConvertFun2 <- .GetAggregConvertFun(names(cols))
-    if (is.null(recursive)) {
-      if (missing(ConvertFun)) {
-        stop("'ConvertFun' argument should provided if 'recursive = NULL'")
-      } else if (!is.na(ConvertFun)) {
-        ConvertFun2 <- rep(ConvertFun, length(cols))
-      }
+    # Treating aggregation at root level
+    if (is.na(ConvertFun)) {
+      ConvertFun2 <- .GetAggregConvertFun(names(cols), Format)
+    } else {
+      ConvertFun2 <- rep(ConvertFun, length(cols))
     }
     dfOut <- SeriesAggreg(cbind(DatesR, as.data.frame(cols)),
                           Format,
@@ -97,17 +105,25 @@ SeriesAggreg.list <- function(x,
       dfCols <- c(dfCols, x[sapply(x, inherits, "matrix")])
       listCols <- listCols[setdiff(names(listCols), names(dfCols))]
       if (length(listCols) > 0) {
-        # Check for predefined ConvertFun for all sub-elements
-        ConvertFun <- .GetAggregConvertFun(names(listCols))
+        if (is.na(ConvertFun)) {
+          # Check for predefined ConvertFun for all sub-elements
+          listConvertFun <- .GetAggregConvertFun(names(listCols), Format)
+        }
         # Run SeriesAggreg for each embedded list
-        listRes <- lapply(names(listCols), function(x) {
-          listCols[[x]]$DatesR <- DatesR
-          SeriesAggreg(listCols[[x]],
+        listRes <- lapply(names(listCols), function(y) {
+          listCols[[y]]$DatesR <- DatesR
+          if (is.na(ConvertFun)) {
+            SeriesAggreg.list(listCols[[y]],
                        Format = Format,
-                       except = except,
-                       ConvertFun = ConvertFun[x],
                        recursive = NULL,
-                       ...)
+                       ...,
+                       ConvertFun = listConvertFun[y])
+          } else {
+            SeriesAggreg.list(listCols[[y]],
+                              Format = Format,
+                              recursive = NULL,
+                              ...)
+          }
         })
         names(listRes) <- names(listCols)
         if (is.null(res$DatesR)) {
@@ -133,10 +149,14 @@ SeriesAggreg.list <- function(x,
               "), it will be ignored in the aggregation"
             )
           } else {
-            ConvertFun <- rep(.GetAggregConvertFun(key), ncol(m))
-            res[[key]] <- SeriesAggreg(data.frame(DatesR, m),
+            if (is.na(ConvertFun)) {
+              ConvertFun2 <- rep(.GetAggregConvertFun(key, Format), ncol(m))
+            } else {
+              ConvertFun2 <- rep(ConvertFun, ncol(m))
+            }
+            res[[key]] <- SeriesAggreg.data.frame(data.frame(DatesR, m),
                                        Format = Format,
-                                       ConvertFun = ConvertFun)
+                                       ConvertFun = ConvertFun2)
           }
         }
       }
diff --git a/R/UtilsSeriesAggreg.R b/R/UtilsSeriesAggreg.R
index d610b03b9aa8a3e9ae278087ea8afb59d8b9aee9..9b47efac0bdc70bd1cabbf3cf6db3e6ed1c0603e 100644
--- a/R/UtilsSeriesAggreg.R
+++ b/R/UtilsSeriesAggreg.R
@@ -36,7 +36,7 @@
          Y = "yearly")
 }
 
-.GetAggregConvertFun <- function(x) {
+.GetAggregConvertFun <- function(x, Format) {
   AggregConvertFunTable <- rbind(
     data.frame(ConvertFun = "mean",
                x = c("Prod", "Rout", "Exp", "SnowPack", "ThermalState",
@@ -53,5 +53,8 @@
     iRes <- AggregConvertFunTable$ConvertFun[AggregConvertFunTable$x == iX]
     iRes <- ifelse(test = any(is.na(iRes)), yes = NA, no = iRes) # R < 4.0 compatibility
   })
+  if(Format %in% c("%d", "%m")) {
+    res <- rep("mean", length(res))
+  }
   return(res)
 }
diff --git a/man/SeriesAggreg.Rd b/man/SeriesAggreg.Rd
index 2c095b545f4e10a5fb985e31a4e55c35fec73d69..e5d871b75a1850cf96fa77ad037bbd41a3544028 100644
--- a/man/SeriesAggreg.Rd
+++ b/man/SeriesAggreg.Rd
@@ -25,6 +25,11 @@ Warning: on the aggregated outputs, the dates correspond to the beginning of the
   \code{\link{SeriesAggreg.InputsModel}} and \code{\link{SeriesAggreg.OutputsModel}}
   call \code{\link{SeriesAggreg.list}} which itself calls \code{\link{SeriesAggreg.data.frame}}.
   So, all arguments passed to any \code{\link{SeriesAggreg}} method will be passed to \code{\link{SeriesAggreg.data.frame}}.
+
+  Argument \code{ConvertFun} also supports quantile calculation by using the syntax "Q[nn]" with [nn] the requested percentile.
+  E.g. use "Q90" for calculating 90th percentile in the aggregation.
+  The formula used is: \code{quantile(x, probs = perc / 100, type = 8, na.rm = TRUE)}.
+
 }
 
 
@@ -47,9 +52,9 @@ Warning: on the aggregated outputs, the dates correspond to the beginning of the
              recursive = TRUE,
              \dots)
 
-\method{SeriesAggreg}{InputsModel}(x, \dots)
+\method{SeriesAggreg}{InputsModel}(x, Format, \dots)
 
-\method{SeriesAggreg}{OutputsModel}(x, \dots)
+\method{SeriesAggreg}{OutputsModel}(x, Format, \dots)
 }
 
 
@@ -62,7 +67,7 @@ Warning: on the aggregated outputs, the dates correspond to the beginning of the
 
 \item{NewTimeFormat}{(deprecated) [character] output time step format (i.e. \code{"hourly"}, \code{"daily"}, \code{"monthly"} or \code{"yearly"}). Use the \code{Format} argument instead}
 
-\item{ConvertFun}{[character] names of aggregation functions (e.g. for P[mm], T[degC], Q[mm]: \code{ConvertFun = c("sum", "mean", "sum"})) (default: use the name of the column (see details) or \code{"mean"} for regime calculation)}
+\item{ConvertFun}{[character] names of aggregation functions (e.g. for P[mm], T[degC], Q[mm]: \code{ConvertFun = c("sum", "mean", "sum"})) or name of aggregation function to apply to all elements if the parameter 'x' is a [list] (See details)}
 
 \item{YearFirstMonth}{(optional) [numeric] integer used when \code{Format = "\%Y"} to set when the starting month of the year (e.g. 01 for calendar year or 09 for hydrological year starting in September)}
 
@@ -77,7 +82,6 @@ Warning: on the aggregated outputs, the dates correspond to the beginning of the
 \item{\dots}{Arguments passed to \code{\link{SeriesAggreg.list}} and then to \code{\link{SeriesAggreg.data.frame}}}
 }
 
-
 \value{
 [POSIXct+numeric] data.frame containing a vector of aggregated dates (POSIXct) and time series values numeric)
 }
diff --git a/tests/testthat/test-SeriesAggreg.R b/tests/testthat/test-SeriesAggreg.R
index f71db81857c4531093ed3c6d123dcbea053dfe45..c9edc42ee8d8ba01b8a4bb7f517f096dd0e5ade7 100644
--- a/tests/testthat/test-SeriesAggreg.R
+++ b/tests/testthat/test-SeriesAggreg.R
@@ -71,6 +71,18 @@ test_that("Check SeriesAggreg output values on yearly aggregation", {
   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, "%m")$Precip,
+               SeriesAggreg(BasinObs[, c("DatesR", "P")], "%m", ConvertFun = "mean")$P)
+})
+
 test_that("No DatesR should warning", {
   TabSeries <- list(
     Dates = BasinObs$DatesR,
@@ -78,7 +90,10 @@ test_that("No DatesR should warning", {
     E = BasinObs$E,
     Qmm = BasinObs$Qmm
   )
-  expect_warning(SeriesAggreg(TabSeries, "%Y%m"), regexp = "has been automatically chosen")
+  expect_warning(
+    SeriesAggreg(TabSeries, "%Y%m", ConvertFun = "sum"),
+    regexp = "has been automatically chosen"
+  )
 })
 
 test_that("Check SeriesAggreg.list 'DatesR' argument", {
@@ -187,6 +202,7 @@ test_that("Check data.frame handling in SeriesAggreg.list", {
   expect_warning(SeriesAggreg(InputsModelDown1, "%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,
@@ -198,6 +214,7 @@ test_that("SeriesAggreg from and to the same time step should return initial tim
   expect_warning(SeriesAggreg(I2, "%Y%m"), regexp = "No time-step conversion was performed")
   expect_equal(I2, suppressWarnings(SeriesAggreg(I2, "%Y%m")))
 })
+
 test_that("SeriesAggreg.data.frame with first column not named DatesR should work",
           {
             expect_warning(SeriesAggreg(
@@ -207,3 +224,28 @@ test_that("SeriesAggreg.data.frame with first column not named DatesR should wor
             ),
             regexp = NA)
           })
+
+test_that("SeriesAggreg should work with ConvertFun 'min', 'max' and 'median'", {
+  Qls <- BasinObs[, c("DatesR", "Qls")]
+  test_ConvertFunRegime <- function(x, ConvertFun, TimeFormat) {
+    expect_equal(nrow(SeriesAggreg(x, TimeFormat, ConvertFun = ConvertFun)),
+                  length(unique(format(BasinObs$DatesR, "%Y"))))
+  }
+  lapply(c("max", "min", "median"), function(x) {test_ConvertFunRegime(Qls, x, "%Y")})
+})
+
+test_that("Error on convertFun Q without 0-100", {
+  Qls <- BasinObs[, c("DatesR", "Qls")]
+  expect_error(SeriesAggreg(Qls, "%Y", "q101"))
+  expect_error(SeriesAggreg(Qls, "%Y", "q-2"))
+  expect_error(SeriesAggreg(Qls, "%Y", "q12.5"))
+})
+
+test_that("ConvertFun q50 should be equal to median", {
+  Qls <- BasinObs[, c("DatesR", "Qls")]
+  expect_equal(SeriesAggreg(Qls, "%Y", "q50"),
+               SeriesAggreg(Qls, "%Y", "median"))
+  expect_equal(SeriesAggreg(Qls, "%Y", "q50"),
+               SeriesAggreg(Qls, "%Y", "q050"))
+})
+