Understanding rainfall-runoff modelling
The role of model components and parameters
Rainfall-runoff models are composed of different elements, e.g. conceptual reservoirs or unit hydrographs, whose behaviour is defined by equations and parameters. While equations are rarely modified by users, parameter values are a key lever that users can activate to taylor the models to their case studies. Understanding the role of model components and parameters is therefore an unavoidable preliminary step to performing hydrological modelling. The airGRteaching package allows such manipulations.
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 airGRteaching and to produce rainfall-runoff transformations considering different model parameter values.
The GR4J model (Perrin, Michel, and Andréassian 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.
# Package loading
library(airGRdatasets)
# Catchment data loading
data("B222001001", package = "airGRdatasets")
# Observed daily time series
ts_obs <- B222001001$TS
# Data processing for GR4J (without Q)
prep_no_q <- PrepGR(DatesR = ts_obs$Date,
Precip = ts_obs$Ptot,
PotEvap = ts_obs$Evap,
HydroModel = "GR4J",
CemaNeige = FALSE)
# Simulation period
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.
# 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)
param_gr4j <- expand.grid(X1 = 350,
X2 = param_x2,
X3 = 90,
X4 = 1.4)
# Streamflow simulations using parameter sets
sim_x2 <- apply(param_gr4j, MARGIN = 1, FUN = function(i_param_gr4j) {
i_sim <- SimGR(PrepGR = prep_no_q,
Param = i_param_gr4j,
SimPer = per_sim,
verbose = FALSE)
i_sim$OutputsModel$Qsim
})
# Graphical comparison
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 [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/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.
# 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)
param_gr4j <- expand.grid(X1 = 350,
X2 = 0,
X3 = 90,
X4 = param_x4)
# Streamflow simulations using parameter sets
sim_x4 <- apply(param_gr4j, MARGIN = 1, FUN = function(i_param_gr4j) {
i_sim <- SimGR(PrepGR = prep_no_q,
Param = i_param_gr4j,
SimPer = per_sim,
verbose = FALSE)
i_sim$OutputsModel$Qsim
})
# Graphical comparison
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 [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 [d]",
bg = "white")
The relative importance of the production and routing functions depends on the time step considered for the rainfall-runoff simulation. The production processes are more important for the larger time steps (e.g. month or year) since it controls the catchment water balance. This can be easily illustrated by aggregating at an annual scale simulations previously performed at a daily time step. The following commands and figure compare, at an annual scale, the GR4J daily simulations previously performed with different X2 parameter values with the simulations performed with different X4 parameter sets. We can observe that at this annual scale, the impact of considering different X4 parameter values is limited compared to the use of different X2 parameter values.
# Aggregation of the simulated streamflow at the yearly time step
sim_x2_y <- cbind(DatesR = as.POSIXct(prep_no_q$InputsModel$DatesR),
as.data.frame(sim_x2))
sim_x2_y <- SeriesAggreg(x = sim_x2_y,
Format = "%Y",
ConvertFun = rep("sum", ncol(sim_x2_y) - 1))
sim_x4_y <- cbind(DatesR = as.POSIXct(prep_no_q$InputsModel$DatesR),
as.data.frame(sim_x4))
sim_x4_y <- SeriesAggreg(x = sim_x4_y,
Format = "%Y",
ConvertFun = rep("sum", ncol(sim_x4_y) - 1))
# 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 [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",
legend = c("X2", "X4"),
lwd = 2, col = c(median(col_param_x2), median(col_param_x4)),
bg = "white")
Finally, the simulation sensitivity on one given parameter can be
also addressed using the
airGRteaching graphical
user interface launched by the ShinyGR() function. It is
very easy by moving the slider of each parameter to change its value,
and to see the effect in real time on the hydrograph of the simulated
streamflows.
On the need to perform a model warm-up
The conceptual reservoirs that compose rainfall-runoff models must be assigned with initial levels when a simulation is done. The way initial levels are defined can lead to potentially significant model errors. The most convenient way that modelers use to initialize rainfall-runoff models is to perform a warm-up run of the model in order to limit the impact of this unknown.
This issue can be illustrated with airGRteaching by considering different warm-up period lengths. The following commands and figure illustrate a portion of the streamflow simulations obtained considering (i) no warm-up period, (ii) a one-month warm-up period, and (iii) a a four-year warm-up period of the GR4J model.
Note that if WupPer is not set by the user, by default
WupPer = NULL, and in this case it ensures a one-year
warm-up using the time steps preceding the SimPer (if data
are available).
# Warm-up and simulation periods
per_wup1m <- c("2002-12-01", "2002-12-31")
per_wup1y <- c("2002-01-01", "2002-12-31")
per_sim <- c("2003-01-01", "2006-12-31")
# Parameter set
param_gr4j <- c(X1 = 350, X32 = 0, X3 = 90, X4 = 1.4)
# Simulation without warm-up period
sim_wup0d <- SimGR(PrepGR = prep_no_q,
Param = param_gr4j,
WupPer = 0L,
SimPer = per_sim)## No warm up period is used## Warning in SimGR(PrepGR = prep_no_q, Param = param_gr4j, WupPer = 0L, SimPer =
## per_sim): "PrepGR" does not contain any Qobs values. The efficiency criterion
## is not computed# Simulation with a 1-month warm-up period
sim_wup1m <- SimGR(PrepGR = prep_no_q,
Param = param_gr4j,
WupPer = per_wup1m,
SimPer = per_sim)## Warning in SimGR(PrepGR = prep_no_q, Param = param_gr4j, WupPer = per_wup1m, :
## "PrepGR" does not contain any Qobs values. The efficiency criterion is not
## computed# Simulation with a 4-year warm-up period
sim_wup1y <- SimGR(PrepGR = prep_no_q,
Param = param_gr4j,
WupPer = per_wup1y,
SimPer = per_sim)## Warning in SimGR(PrepGR = prep_no_q, Param = param_gr4j, WupPer = per_wup1y, :
## "PrepGR" does not contain any Qobs values. The efficiency criterion is not
## computed# Graphical comparison
col_wup <- c("orchid", "orange2", "green3")
matplot(x = as.POSIXct(sim_wup0d$OutputsModel$DatesR),
y = cbind(sim_wup0d$OutputsModel$Qsim,
sim_wup1m$OutputsModel$Qsim,
sim_wup1y$OutputsModel$Qsim),
xlab = "time [d]", ylab = "flow [mm/d]",
type = "l", lty = 1, lwd = 2, col = col_wup,
xlim = as.POSIXct(x = c("2003-01-01", "2003-09-01"), tz = "UTC"))
legend("topright",
legend = c("no warm-up", "1-month warm-up", "1-year warm-up"),
col = col_wup, lwd = 2,
bg = "white")
