Scope

Discharge simulations are affected by several sources of uncertainty (e.g. errors in meteorological forcings, sub-optimal parameter estimates).

Data assimilation (DA) allows to combine measurements and model simulations under the assumption that both supply useful information to obtain the most likely estimate of the system state. Among the sequential methods, the ensemble-based techniques allow to explicitly handle different sources of uncertainty and to quantify the unknown errors of both model states and observations.

After the study of Piazzi et al. (accepted), airGRdatassim package (Piazzi and Delaigue 2021) provides functions allowing to perform DA-based ensemble discharge simulations. airGRdatassim is a package based on the airGR hydrological modeling package (Coron et al. 2017, 2023). This package is designed to flexibly use the Ensemble Kalman filter (EnKF) or Particle filter (PF) technique to assimilate daily discharge observations into GR4J, GR5J, and GR6J models (use the command ?airGR to get the references of the GR models).

The uncertainty in both the meteorological forcings (i.e. precipitation and potential evapotranspiration) and model states (i.e. production store level, routing store level and unit hydrographs) can be easily taken into account by enabling specific perturbation procedures.

In this user guide we go through the assimilation procedure with the aim of explaining how to perform DA-based discharge simulations. Along with the description of the modelling system, key recommendations are provided to support the definition of the operational configuration.

Loading catchment data

airGRdatassim is a package based on the airGR hydrological modeling package (Coron et al. 2017, 2023).

All the examples rely on the following data set, which is available in airGR package:

library(airGRdatassim)

## Loading required package: airGR

data(L0123001, package = "airGR") 
head(BasinObs)

##       DatesR    P    T   E   Qls    Qmm
## 1 1984-01-01  4.1  0.5 0.2  2640 0.6336
## 2 1984-01-02 15.9  0.2 0.2  3440 0.8256
## 3 1984-01-03  0.8  0.9 0.3 12200 2.9280
## 4 1984-01-04  0.0  0.5 0.3  7600 1.8240
## 5 1984-01-05  0.0 -1.6 0.1  6250 1.5000
## 6 1984-01-06  0.0  0.9 0.3  5650 1.3560

Ensemble size

The ensemble size can be defined by assigning NbMbr with the number of ensemble members.

It is noteworthy that the PF technique is more sensitive to the ensemble size than the EnKF scheme. Indeed, compared to the EnKF scheme, the PF generally needs a larger ensemble size to efficiently perform the weighting and resampling of particles. It is recommended to generate and use an ensemble of at least 20 or 30 members, when using the EnKF or the PF scheme, respectively.

Enable/disable the assimilation via EnKF or PF

The DA methodology can be chosen between EnKF and PF by defining DaMethod. To perform the open-loop simulation (i.e., without DA), the assimilation of discharge observations can be disabled.

If DaMethod = "PF", the PF performs the joint update of all the state variables.

If DaMethod = "EnKF", the EnKF performs the update of the state variables specified in StateEnKF.

If DaMethod = "none", discharge assimilation is not performed (i.e., open-loop simulation).

Among the model states, the update of the routing store level ensures the most benefit from the assimilation of observed discharges. Conversely, the update of the unit hydrograph only has a slight impact on the accuracy of the DA-based discharge simulations (Piazzi et al. accepted).

Model uncertainties

Both the uncetainties in meteorological forcings and state variables can be taken into account in the assimilation scheme.

For both DA techniques, the state variables defined in StatePert are perturbed.

The perturbation of model states is performed through a normally-distributed null-mean noise. The noise variance is assumed equal to the variance of the state variables resulting from the analysis procedure and it is restricted between upper and lower limits (Salamon and Feyen 2009).

When considering the uncertainty in the forcings, the CreateInputsPert() function generates an ensemble for precipitation or/and potential evapotranspiration, provided as function argument/s. It is noteworthy that a significant reduction in the ensemble spread may occur especially in no-rain periods, when accounting only for meteorological uncertainty.

Even though a more comprehensive representation of the state uncertainty can prevent the ensemble shrinkage, it has a different impact on the performance of the PF and the EnKF schemes (Piazzi et al. accepted). Indeed, while a higher model error covariance can lead to an overweighting of observations in the EnKF-based analysis procedure, a larger spread of the ensemble simulations allows for a more straightforward weighting and thus a more efficient resampling of particles in the PF scheme.

Definition of DA settings

The example below considers EnKF-based discharge simulations relying on a 100-member ensemble, when accounting for the uncertainty in model state variables (see StatePert). More specifically, only the levels of the production ("Prod") and the routing ("Rout") stores are updated and perturbed (see StatePert).

It is of key importance to consider that StateEnKF allows to select the state variables to be updated and StatePert identifies the state variables to be perturbed in the EnKF scheme. Conversely, in the PF scheme StatePert identifies the state variables of interest only within the perturbation procedure, as this DA technique jointly updates all the model states.

# ensemble size 
NbMbr <- 100L

# "EnKF" ; "PF" ; "none"
DaMethod <- "EnKF"

# "Prod": production store level
# "Rout": routing store level;
# "UH1": unit hydrograph 1 levels (not defined in GR5J model)
# "UH2": unit hydrograph 2 levels
StateEnKF <- c("Prod", "Rout")
StatePert <- c("Prod", "Rout") 

EnKF-based discharge simulations

The following example supplies probabilistic discharge simulations relying on the EnKF-based assimilation of daily discharge observations.

Param <- c(X1 = 194.243, X2 = -0.088, X3 = 117.740, X4 = 1.680, X5 = 0.000)
ResEnKF <- RunModel_DA(InputsModel = InputsModel,
                       InputsPert = InputsPert,
                       Qobs = BasinObs$Qmm,
                       IndRun = IndRun,
                       FUN_MOD = RunModel_GR5J, Param = Param,
                       DaMethod = "EnKF", NbMbr = NbMbr,
                       StateEnKF = StateEnKF, StatePert = StatePert)

As defined in StateEnKF, in this example the EnKF scheme updates only the levels of the production ("Prod") and the routing ("Rout") stores through the assimilation of discharge observations. Along with the uncertainty in meterological forcings (InputsPert is supplied), also the uncertainty in model states (StatePert) is taken into account. Therefore, the selected state variables (i.e. the store levels) are consistently perturbed after their update.

The EnKF scheme relies on mass-constraints to guarantee the non-negativity of state variables by ensuring minimum state values. It is noteworthy that the level of production store cannot drop below the 5 percent of its capacity. A further constraint prevents the levels of both the production and the routing stores from exceeding their maximum capacities.

PF-based discharge simulations

The example below provides ensemble discharge simulations resulting from the daily PF-based assimilation of discharge observations.

The PF scheme relies on the Sequential Importance Resampling (SIR) technique (Gordon, Salmond, and Smith 1993). According to the SIR approach, the importance weights associated to the particles are updated according to the likelihood of particle states with respect to the observed one. A resampling procedure allows to discard particles having a low probability and to replicate those having a high importance weight. If the ensemble is too squeezed to properly evaluate the importance weights, all particles are assigned equal weights.

ResPF <- RunModel_DA(InputsModel = InputsModel,
                     InputsPert = InputsPert,
                     Qobs = BasinObs$Qmm,
                     IndRun = IndRun,
                     FUN_MOD = RunModel_GR5J, Param = Param, 
                     DaMethod = "PF", NbMbr = NbMbr,
                     StatePert = "Rout")

Unlike the previous example, here only the uncertainty in the level of the routing store (see StatePert) is taken into account. Therefore, the routing store level is the only state variable to be perturbed after the updating procedure.

After the state perturbation, a mass-constraint prevents the levels of both the production and the routing stores from exceeding their maximum capacities.