knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )

The following example is what I consider the first step to understand
how **HBV.IANIGLA** works. As a first approach, I propose you to fit
the streamflow discharge in a synthetic lumped catchment case (in my opinion this
the most simple hydrological modeling case).
The objectives of this excersice are:

- to introduce the building-up of the model.
- to get a first approach to parameter calibration.
- to fit the streamflow discharge in a perfect match case (this does not occur in real world hydrology).
- to make a sensitivity analysis.

Initially, we load the dataset containing all the necessary information to run the model: air temperature, precipitation, potential evapotranspiration and the outflow of our synthetic catchment:

library(HBV.IANIGLA) data("lumped_hbv") head(lumped_hbv) summary(lumped_hbv)

This case is about a catchment without glaciers. In the next lines I will help
you to build-up the model line by line. Because is our first excersice, I will provide
you with the correct *initial conditions* and *parameters* in all modules except for
the **Routing_HBV** function.

Remember that hydrological models are generally build in a top down direction (from precipitation to streamflow routing - note that most hydrological books are structured in the same way).

# we first consider the SnowGlacier module to take into account # precipitation partioning and the snow accumulation/melting. snow_module <- SnowGlacier_HBV(model = 1, inputData = as.matrix( lumped_hbv[ , c('T(ºC)', 'P(mm/d)')] ), initCond = c(20, 2), param = c(1.20, 1.00, 0.00, 2.5) ) # is always advisable to take a look at the output (unless until you get familiarised # with the model) head(snow_module) # now see the function documentation...you should be able to relate the model option # with the matrix output!

Now we will pass the total water produced (rainfall plus snowmelt) to the soil routine. Note that we will use the potential evapotranspiration series,

soil_module <- Soil_HBV(model = 1, inputData = cbind(snow_module[ , "Total"], lumped_hbv$`PET(mm/d)`), initCond = c(100, 1), param = c(200, 0.8, 1.15) ) head(soil_module)

from this module we can get the actual evapotranspiration, soil moisture and recharge
series. In the next step we incorporate the last variable to the routing function.
Remember that the module's parameters (*param* argument) are not calibrated!

routing_module <- Routing_HBV(model = 1, lake = F, inputData = as.matrix(soil_module[ , "Rech"]), initCond = c(0, 0, 0), param = c(0.9, 0.01, 0.001, 0.5, 0.01) ) head(routing_module)

Finally we will apply the transfer function in order to adjust the hydrograph timing,

tf_module <- round( UH(model = 1, Qg = routing_module[ , "Qg"], param = c(1.5) ), 2) # let's plot the "true" and simulated hydrographs plot( x = lumped_hbv[ , "Date"], y = tf_module, type = "l", col = "dodgerblue", xlab = "Date", ylab = "q(mm/d)") lines(x = lumped_hbv[ , "Date"], y = lumped_hbv[ , "qout(mm/d)"], col = "red") # you can get a goodness of fit measurement by, cor(x = lumped_hbv[ , "qout(mm/d)"], y = tf_module) # hint: to use hydrological goodness of fit functions (e.g.: NSE - Nash-Sutcliffe # Efficiency) you can install the hydroGOF package.

As you will note, our HBV model needs some calibration (in he **Routing_HBV**
parameters).

**Now is your turn**

change (manually) the corresponding parameter values to get a

decentstreamflow simulation.Help: the order of magnitude of the actual parameters is correct (e.g.: the real value ofK0is between 0.9 and 0.1 and so on).

**Hints**

When you are changing the parameters 'by hand' is always advisable to plot the results.

Pay attention on how sensible is the result to changes in a single parameter.

In my view, when building an HBV model is always a better idea to construct a
succinct version of it by defining a function. If your are new in the **R**
language I recommend to visit one of the following web pages,

Now lets recycle the previous code to create a general purpose function for lumped basin cases,

## brief arguments desription # basin: data frame with the same structure of the data("lumped_hbv") (colnames included). # param_snow: numeric vector with snow module parameters. # param_soil: numeric vector with soil moisture parameters. # param_routing: numeric vector with the routing parameters. # param_tf: numeric vector with the transfer function parameter. # init_snow: numeric value with initial snow water equivalent. Default value being 20 mm. # init_soil: numeric value with initial soil moisture content. Default value being 100 mm. # init_routing: numeric vector with bucket water initial values. Default values are 0 mm. ## output # simulated streamflow series. hbv_lumped <- function(basin, param_snow, param_soil, param_routing, param_tf, init_snow = 20, init_soil = 100, init_routing = c(0, 0, 0) ){ snow_module <- SnowGlacier_HBV(model = 1, inputData = as.matrix( basin[ , c('T(ºC)', 'P(mm/d)')] ), initCond = c(init_snow, 2), param = param_snow ) soil_module <- Soil_HBV(model = 1, inputData = cbind(snow_module[ , "Total"], basin$`PET(mm/d)`), initCond = c(init_soil, 1), param = param_soil ) routing_module <- Routing_HBV(model = 1, lake = F, inputData = as.matrix(soil_module[ , "Rech"]), initCond = init_routing, param = param_routing ) tf_module <- UH(model = 1, Qg = routing_module[ , "Qg"], param = param_tf ) out <- round(tf_module, 2) return(out) }

After running this chunk of code you should see the function name in the **Global Environment**, now you can use it.

**Advice**

When I use my own functions I usually write them in a separate R file. This helps me to keep my main script clean. To call the function you just run the single code line

source(file = "my_hbv_model.R")andvoila!!!...the function is ready to use.

Note that we have defined default values for all the initial conditions. This means
that when running the function, you don't have to give them a value. For the sake
of simplicity I gave them the correct values...so as this is our first case example
forget about them and **focus on the other function arguments**.

**Now is your turn**

- If you didn't understand the logic behind the function is a good idea to take a break (and
maybe a breath) and revisit the (first) recommended web page R bloggers. It is advisable
to reproduce and
*play*a little bit with the web page examples.**Programming is a necessary condition to learn how to program!**. - Analyze carefully how the function arguments and modules'
outputs are used by the subsequent HBV functions, and how the function's returning value
is set. Without looking at the
*hbv_lumped*example, can you do it by your own? - Now suppose that you want to use the
**Routing_HBV**model number two: what do you need to modify? is the*param_routing*argument of the same length? what about the*init_routing*argument?

In this section we are going to make something that looks more closely to what we should do in real-world problems: work with the arguments of an HBV function. Imagine what will happen if you decide to build-up a model every time you work in a modeling problem...it will be too prone error! To avoid pitfall it would be a better idea to construct a general purpose function.

In the following lines we are going to use our *hbv_lumped* model to calibrate
the routing parameters in order to fit the simulated streamflow discharge.

# one of the ideas behind a function is to succinctly express # operations that will look complicated or too large. In our # case, by modifying (in the same line of arguments) just some # parameters we obtain the basin discharge streamflow <- hbv_lumped(basin = lumped_hbv, param_snow = c(1.20, 1.00, 0.00, 2.5), param_soil = c(200, 0.8, 1.15), param_routing = c(0.2, 0.01, 0.005, 0.75, 0.15), param_tf = c(1.5)) plot( x = lumped_hbv[ , "Date"], y = streamflow, type = "l", col = "dodgerblue", xlab = "Date", ylab = "q(mm/d)") lines(x = lumped_hbv[ , "Date"], y = lumped_hbv[ , "qout(mm/d)"], col = "red") cor(x = streamflow, y = lumped_hbv[ , "qout(mm/d)"] )

Now you have a compact expression for your modeling problem. Recycling this kind of functions
will also save you a lot of time. Be also aware that you only have to build them **once**!

**Now is your turn**

- if you haven't, get the perfect (or unless an almost perfect) streamflow fit.
- change another parameter(s) and see what happen with the hydrograph.

According to @beven:2008 a sensitivity analysis is:

the response of a model to a change in a parameter or input variable. This reaction can either be assesed locally in the model space or over a global sample of points.

In this part of our **lumped model vignette**, I will give you the answer to our calibration
problem and we are going to make a sensitivity analysis on them.

**Rounting_HBV parameters:** `c(0.1, 0.05, 0.002, 0.9, 0.1)`

One way to look at the sensitivity of a single parameter, is to change its value in a reasonable range to evaluate some model response(s) (in our case, the basin discharge) [@pianosi:2016]. In this example our evaluation measurement will be the Pearson correlation coefficient.

# first we are going to generate a vector with the possible values # of our parameter. Let's say that we want to evaluate the effect of # changing the K0 value. target_param <- seq(from = 0.1, to = 0.9, by = 0.01) # now we get the number fo required iterations and we create an empty # vector to store the evaluation function results n_it <- length(target_param) response <- c() # finally run our function n_it times and we plot the results for(i in 1:n_it){ streamflow <- hbv_lumped(basin = lumped_hbv, param_snow = c(1.20, 1.00, 0.00, 2.5), param_soil = c(200, 0.8, 1.15), param_routing = c(target_param[i], 0.05, 0.002, 0.9, 0.1), param_tf = c(1.5)) response[i] <- cor(x = streamflow, y = lumped_hbv[ , "qout(mm/d)"] ) } plot(x = target_param, y = response, type = "p", ylim = c(0, 1), pch = 20, cex = 0.5)

**Now is your turn**

- explore what happen if the target value is another parameter.
- try not only with the
**Routing_HBV**parameters, but with the other modules. - which model parameters are the most sensible for streamflow discharge?

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.