Scope

In the airGR package, the classical way to calibrate a model is to use the Michel’s algorithm (see the Calibration_Michel() function).

The Michel’s algorithm combines a global and a local approach. A screening is first performed either based on a rough predefined grid (considering various initial values for each parameter) or from a list of initial parameter sets. The best set identified in this screening is then used as a starting point for the steepest descent local search algorithm.

In some specific situations, for example if the calibration period is too short and by consequence non representative of the catchment behaviour, a local calibration algorithm can give poor results.

In this vignette, we show a method using a known parameter set that can be used as an alternative for the grid-screening calibration procedure, and we well compare to two methods using the Calibration_Michel() function. The generalist parameters sets introduced here are taken from Andréassian et al. (2014).

Data preparation

We load an example data set from the package and the GR4J parameter sets.

## loading catchment data
data(L0123001)

## loading generalist parameter sets
data(Param_Sets_GR4J)

The given GR4J X4u variable does not correspond to the actual GR4J X4 parameter. As explained in Andréassian et al. (2014, sec. 2.1), the given GR4J X4u value has to be adjusted (rescaled) using catchment area (S) [km2] as follows: X4 = X4u / 5.995 * S^0.3 (please= note that the formula is erroneous in the publication).

It means we need first to transform the X4 parameter.

Param_Sets_GR4J$X4 <- Param_Sets_GR4J$X4u / 5.995 * BasinInfo$BasinArea^0.3
Param_Sets_GR4J$X4u <- NULL
Param_Sets_GR4J <- as.matrix(Param_Sets_GR4J)

Please, find below the summary of the 27 sets of the 4 parameters.

##        X1             X2              X3            X4      
##  Min.   : 126   Min.   :-54.5   Min.   :  8   Min.   :1.21  
##  1st Qu.: 208   1st Qu.: -2.0   1st Qu.: 35   1st Qu.:1.75  
##  Median : 291   Median : -1.1   Median : 76   Median :2.10  
##  Mean   : 471   Mean   : -3.4   Mean   : 90   Mean   :2.09  
##  3rd Qu.: 359   3rd Qu.: -0.6   3rd Qu.:106   3rd Qu.:2.45  
##  Max.   :4006   Max.   :  0.8   Max.   :318   Max.   :3.47

Object model preparation

We assume that the R global environment contains data and functions from the Get Started vignette.

The calibration period has been defined from 1990-01-01 to 1990-02-28, and the validation period from 1990-03-01 to 1999-12-31.

As a consequence, in this example the calibration period is very short, less than 6 months.

## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR, 
                                 Precip = BasinObs$P, PotEvap = BasinObs$E)

## ---- calibration step

## short calibration period selection (< 6 months)
Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/1990 00:00"), 
               which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="28/02/1990 00:00"))

## preparation of the RunOptions object for the calibration period
RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Cal)

## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit_Cal  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, 
                                    RunOptions = RunOptions_Cal, Obs = BasinObs$Qmm[Ind_Cal])


## ---- validation step

## validation period selection
Ind_Val <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/03/1990 00:00"), 
               which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="31/12/1999 00:00"))

## preparation of the RunOptions object for the validation period
RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Val)

## efficiency criterion (Nash-Sutcliffe Efficiency) on the validation period
InputsCrit_Val  <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, 
                                    RunOptions = RunOptions_Val, Obs = BasinObs$Qmm[Ind_Val])

GR4J parameter distributions quantiles used in the grid-screening step

By default, the predefined grid used by the Calibration_Michel() function contains parameters quantiles computed after recursive calibrations on 1200 basins (from Australia, France and USA).

CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)

The parameter set corresponding to the best criterion is the following:

##     X1     X2     X3     X4 
## 165.05   6.46 422.50   2.45

The Nash-Sutcliffe Efficiency criterion computed on the calibration period is better (0.495) than with the generalist parameter sets, but the one computed on the validation period is lower (0.57). This shows that the generalist parameter sets give more robust model in this case.

GR4J parameter sets used in the grid-screening step

It is also possible to give to the CreateCalibOptions() function a matrix of parameter sets used for the grid-screening calibration procedure. So, it possible is to use by this way the GR4J generalist parameter sets.

CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel,
                                   StartParamList = Param_Sets_GR4J)

Here is the parameter set corresponding to the best criteria found.

##     X1     X2     X3     X4 
## 161.62   6.53 429.00   2.45

The results are the same here. The Nash-Sutcliffe Efficiency criterion computed on the calibration period is better (0.495), but the one computed on the validation period is just a little bit lower (0.568) than the classical calibration.

Generally, the advantage of using GR4J parameter sets rather than the GR4J generalist parameter quantiles is that they make more sense than a simple exploration of combinations of quantiles of parameter distributions (each parameter set represents a consistent ensemble). In addition, for the first step, the number of iterations is smaller (27 runs instead of 81), which can save time if one wants to run a very large number of calibrations.