airGRdatassim (version 0.1.4)

RunModel_DA: Run discharge simulations with ensemble-based data assimilation

Description

Function which performs discharge ensemble simulations with the assimilation of observed discharges through the Ensemble Kalman filter (EnKF) or the Particle filter (PF) schemes. More information about the efficiency of these data assimilation schemes with GR5J can be found in Piazzi et al. (2021).

Usage

RunModel_DA(InputsModel, InputsPert = NULL, Qobs = NULL,
            IndRun,
            FUN_MOD, Param,
            DaMethod = c("EnKF", "PF", "none"), NbMbr = NULL,
            StateEnKF = NULL, StatePert = NULL,
            Seed = NULL)
# S3 method for OutputsModelDA
[(x, i)
# S3 method for OutputsModelDA
plot(x, Qobs = NULL, main = NULL,
     ColSim = "orangered", ColObs = par("fg"), ...)

Value

[list] 'OutputsModel' containing the outputs organised as follows:

DatesR

[POSIXlt] series of dates (length of IndRun)

QsimEns

[numeric] matrix (dim(NbTime, NbMbr)) of ensemble discharge simulations

EnsStateBkg

[numeric] array (dim(NbTime, NbMbr, Nstate)) of ensemble values of background model states (before the filter update)

EnsStateA

[numeric] array (dim(NbTime, NbMbr, Nstate)) of ensemble values of analysis model states (after the filter update)

NbTime

[integer] atomic vector of length of IndRun

NbMbr

[integer] atomic vector of number of ensemble members

NbState

[integer] atomic vector of number of of model states

On the graphical outputs:

- solid line: medians of the input values

- polygon: minima and maxima of the input values

Arguments

InputsModel

[list] object of class InputsModel containing the data required to run the model (see CreateInputsModel for details)

InputsPert

(optional) [list] object of class InputsModel containing the ensembles of perturbed data required to evaluate the model probabilistic outputs, if the uncertainty in meteorological forcings is taken into account. The variables data proposed in InputsPert will erase the variables data given in InputsModel

Qobs

(optional) [numeric] time series of observed flow [mm/d]

IndRun

[numeric] index of period to be used for the model run [-]

FUN_MOD

[function] daily hydrological model functions (GR4J, GR5J or GR6J; e.g. RunModel_GR5J, RunModel_CemaNeigeGR5J)

Param

[numeric] vector of model parameters (number of parameters depends on the used hydrological model, from 4 parameters in GR4J up to 10 parameters in GR6J with CemaNeige)

DaMethod

[character] name of the data assimilation technique:

"EnKF"	discharge assimilation via Ensemble Kalman filter
"PF"	discharge assimilation via Particle filter
"none"	discharge assimilation is not performed (i.e. open-loop simulation)

NbMbr

(optional) [numeric] number of ensemble members (minimum of 20 recommanded for the EnKF scheme and of 30 for the PF scheme; by default=NULL: 50 if InputsPert is not set or InputsPert$NbMbr otherwise)

StateEnKF

[character] vector of the names of state variables to be updated via EnKF:

"Prod"	level of the production store [mm]
"Rout"	level of the routing store [mm]
"UH1"	unit hydrograph 1 levels [mm] (not defined for the GR5J model)
"UH2"	unit hydrograph 2 levels [mm]

StatePert

[character] vector of the names of state variables to be perturbed via EnKF or PF:

"Prod"	level of the production store [mm]
"Rout"	level of the routing store [mm]
"UH1	unit hydrograph 1 levels [mm] (not defined for the GR5J model)
"UH2"	unit hydrograph 2 levels [mm]

Seed

(optional) [numeric] seed of random number generator

x

[OutputsModelDA] containing the vector of dates (POSIXt) and the time series of numeric values DA model outputs

i

[integer] of the indices to subset a time series or [character] names of the elements to extract

main

(optional) [character] an overall title for the plot (see title)

ColSim, ColObs

(optional) [character] color to be used for simulated flow and observed flow (in any format that col2rgb accepts)

...

other parameters to be passed through to plotting functions

Author

Gaia Piazzi, Olivier Delaigue

Details

Discharge observations are sequentially assimilated at each time step (if DaMethod != "none") via the EnKF or PF schemes. Because of the sequential approach, the analysis states resulting from the assimilation procedure at time step t are used as initial states at the following prediction time step t + 1.

When accounting for the uncertainty in model inputs (see InputsPert), the ensemble discharge simulations are driven by perfect meteorological forecasts, which are generated according to the methodology proposed by Clark et al. (2008) (see CreateInputsPert for details).

It is possible to enable/disable the update of specific state variables via EnKF by defining/not defining their names in StateEnKF. Please note that StateEnKF is therefore required only when using the EnKF.

Both the PF and EnKF can account for uncertainty in model states by perturbing the state variables specified in StatePert (Salamon and Feyen, 2009). It is noteworthy that, when using the EnKF, the perturbation is allowed only for the state variables included in StateEnKF. If the state uncertainty is taken into account, the initial states at the prediction time step t + 1 are the perturbed analysis states resulting from the assimilation and perturbation procedures at time step t.

In order to ensure reproducible results, Seed can be set to fix the randomness in the generation of perturbations.

For further details and guidelines on the choice of the DA technique, see the references section.

Nota: The function can be applied when using GR4J, GR5J and GR6J models (i.e. daily model time step), with or withouth the CemaNeige module (see airGR package).

References

- Clark, M. P., Rupp, D. E., Woods, R. A., Zheng, X., Ibbitt, R. P., Slater, A. G. et al. (2008). Hydrological data assimilation with the ensemble Kalman filter: Use of streamflow observations to update states in a distributed hydrological model. Advances in Water Resources, 31(10), 1309-1324, tools:::Rd_expr_doi("10.1016/j.advwatres.2008.06.005")

- Piazzi, G., Thirel, G., Perrin, C. and Delaigue, O. (2021). Sequential data assimilation for streamflow forecasting: assessing the sensitivity to uncertainties and updated variables of a conceptual hydrological model at basin scale. Water Resources Research, 57, tools:::Rd_expr_doi("10.1029/2020WR028390").

- Salamon, P. and Feyen, L. (2009). Assessing parameter, precipitation, and predictive uncertainty in a distributed hydrological model using sequential data assimilation with the particle filter. Journal of Hydrology, 376(3-4), 428-442, tools:::Rd_expr_doi("10.1016/j.jhydrol.2009.07.051")

See Also

CreateInputsPert

Examples

library(airGRdatassim)

## loading catchment data
data(L0123001, package = "airGR")
Param <- c(X1 = 194.243, X2 = -0.088, X3 = 117.740, X4 = 1.680, X5 = 0.000)

## run period selection
IndRun <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="2006-01-01"),
              which(format(BasinObs$DatesR, format = "%Y-%m-%d")=="2006-01-31"))

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

## number of ensemble members
## minimum of 20 recommanded for the EnKF scheme
## minimum of 30 recommanded for the PF scheme
NbMbr <- 20L

## preparation of perturbed meteorological ensemble
InputsPert <- CreateInputsPert(FUN_MOD = RunModel_GR5J,
                               DatesR = BasinObs$DatesR,
                               Precip = BasinObs$P, PotEvap = BasinObs$E,
                               NbMbr = NbMbr)

## simulation with DA via EnKF
OutputsModelDA <- RunModel_DA(InputsModel = InputsModel,
                              InputsPert = InputsPert,
                              Qobs = BasinObs$Qmm,
                              IndRun = IndRun,
                              FUN_MOD = RunModel_GR5J, Param = Param,
                              DaMethod = "EnKF", NbMbr = NbMbr,
                              StateEnKF = c("Prod", "Rout"),
                              StatePert = c("Prod", "Rout"))

## results preview
plot(OutputsModelDA, Qobs = BasinObs$Qmm[IndRun])

## results preview on a subset
plot(OutputsModelDA[1:10], Qobs = BasinObs$Qmm[IndRun][1:10])