Run a hydrological model (minimalistic example)

This vignette shows how to run a hydrological model as a minimalistic example.

Preparations

Load packages

library(openQUARREL)
library(dplyr)

Create a minimum input

# input_data is an example input loaded with openQUARREL
minimum_input <- input_data %>% 
  # filter to a single catchment, i.e. Thur at Jonschwil
  filter(HSU_ID == "2303") %>% 
  # select only a few columns
  select(DatesR, P, T, E)

# convert Date to POSIXct
minimum_input$DatesR <- as.POSIXct(minimum_input$DatesR)
attr(minimum_input$DatesR, "tzone") <- "UTC"

The data frame minimum_input contains columns P for precipitation in mm/d, T for air temperature in degrees Celsius and E for potential evapotranspiration in mm/d. See also the great airGR get started documentation. If you need to calculate E, there are functions available, e.g. in airGR::PE_Oudin.

DatesR	P	T	E
1981-01-01	7.08	-0.42	0.60
1981-01-02	7.82	-2.55	0.32
1981-01-03	27.37	2.43	0.98
1981-01-04	22.56	-1.16	0.51
1981-01-05	8.68	-4.86	0.02
1981-01-06	14.99	-4.00	0.13

Choose a hydrological model

# default_cal_par are default calibration parameter loaded with openQUARREL
# todo: cal_par needs to be a global parameter currently, needs to be changed in future releases
cal_par <- default_cal_par 
# todo: topmodel is not running right now as the required BasinInfo is not provided in example data
model <- "TUW"

Set model parameters

Create a parameter set, which is here the middle points of the provided intervals. See TUWmodel documentation for parameter explanation. The list default_cal_par is loaded with library(openQUARREL) and holds model and calibration specific settings.

param <- (cal_par[[model]]$upper + default_cal_par[[model]]$lower) / 2
print(param)
#>    SCF    DDF     Tr     Ts     Tm  LPrat     FC   BETA     k0     k1     k2 
#>    1.2    2.5    2.0   -1.0    0.0    0.5  300.0   10.0    1.0   16.0  140.0 
#>   lsuz  cperc   bmax croute 
#>   50.5    4.0   15.0   25.0

Create input

This is done specifically for each model/package requirement, here TUWmodel.

# todo: the last input can be empty, will be not required in a future release 
input <- create_input(model, minimum_input, list())

Run the model

sim <- simulate_model(model, param, input)

The output is a list with fields date, Qsim, Qobs, SWE, psolid, pliquid, melt, more_info. You can access typical components as simulated discharge Qsim, and if provided in the input also the observed discharge Qobs. If a snow module is integrated in the model or chosen to be added (see vignette("include_snow_module")) snow water equivalent SWE, the partitioned precipitation fluxes pliquid and psolid, and the snow melt flux. The element more_info is specific to the chosen model and package, in this case TUWmodel.

str(sim)
#> List of 8
#>  $ date     : Date[1:14853], format: "1981-01-01" "1981-01-02" ...
#>  $ Qsim     : num [1:14853] 0.000362 0.001447 0.003247 0.005758 0.008975 ...
#>  $ Qobs     : NULL
#>  $ SWE      : num [1:14853] 6.85 16.24 10.16 37.23 47.65 ...
#>  $ psolid   : num [1:14853] 5.71 7.82 0 22.56 8.68 ...
#>  $ pliquid  : num [1:14853] 1.37 0 27.37 0 0 ...
#>  $ melt     : num [1:14853] 0 0 6.08 0 0 ...
#>  $ more_info:List of 1
#>   ..$ output_model:List of 22
#>   .. ..$ itsteps: int 14853
#>   .. ..$ nzones : int 1
#>   .. ..$ area   : num 1
#>   .. ..$ param  : num [1, 1:15] 1.2 2.5 2 -1 0 0.5 300 10 1 16 ...
#>   .. .. ..- attr(*, "dimnames")=List of 2
#>   .. .. .. ..$ : NULL
#>   .. .. .. ..$ : chr [1:15] "SCF" "DDF" "Tr" "Ts" ...
#>   .. .. ..- attr(*, "names")= chr [1:15] "SCF" "DDF" "Tr" "Ts" ...
#>   .. ..$ incon  : num [1, 1:4] 50 0 2.5 2.5
#>   .. .. ..- attr(*, "names")= chr [1:4] "SSM0" "SWE0" "SUZ0" "SLZ0"
#>   .. ..$ prec   : num [1:14853] 7.08 7.82 27.37 22.56 8.68 ...
#>   .. ..$ airt   : num [1:14853] -0.42 -2.55 2.43 -1.16 -4.86 ...
#>   .. ..$ ep     : num [1:14853] 0.6 0.32 0.98 0.51 0.02 0.13 0 0 0 0.28 ...
#>   .. ..$ output : num [1, 1:20, 1:14853] 3.62e-04 6.85 5.14e+01 1.37 5.71 ...
#>   .. ..$ qzones : num [1, 1:14853] 0.000362 0.001447 0.003247 0.005758 0.008975 ...
#>   .. ..$ q      : num [1, 1:14853] 0.000362 0.001447 0.003247 0.005758 0.008975 ...
#>   .. ..$ swe    : num [1, 1:14853] 6.85 16.24 10.16 37.23 47.65 ...
#>   .. ..$ melt   : num [1, 1:14853] 0 0 6.08 0 0 ...
#>   .. ..$ q0     : num [1, 1:14853] 0 0 0 0 0 0 0 0 0 0 ...
#>   .. ..$ q1     : num [1, 1:14853] 0 0 0 0 0 0 0 0 0 0 ...
#>   .. ..$ q2     : num [1, 1:14853] 0.0355 0.0352 0.035 0.0347 0.0345 ...
#>   .. ..$ moist  : num [1, 1:14853] 51.4 51.4 84.3 84.3 84.3 ...
#>   .. ..$ rain   : num [1, 1:14853] 1.37 0 27.37 0 0 ...
#>   .. ..$ snow   : num [1, 1:14853] 5.71 7.82 0 22.56 8.68 ...
#>   .. ..$ eta    : num [1, 1:14853] 0 0 0.554 0 0 ...
#>   .. ..$ suz    : num [1, 1:14853] 0 0 0 0 0 0 0 0 0 0 ...
#>   .. ..$ slz    : num [1, 1:14853] 4.96 4.93 4.89 4.86 4.83 ...

Plot results

Here is a simple plot for simulated discharge with skipping the first ~3 years (warm-up):

plot(sim$date[3*365:length(sim$date)], sim$Qsim[3*365:length(sim$date)], type = "l", main = "Simulated Discharge", ylab = "Q [mm/d]", xlab = "Date")

Next steps

Now you do this again with a different hydrological model (see table in vignette("model_overview")), or proceed with vignette("include_snow_module").

Michael Schirmer