---
title: "Chapter 2: A lumped model"
output: rmarkdown::html_vignette
author: "Ezequiel Toum"
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{Chapter 2: A lumped model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## The problem

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. 

## Data

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:

```{r}
library(HBV.IANIGLA)

data("lumped_hbv")

head(lumped_hbv)
summary(lumped_hbv)

```

## Building the model

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). 

```{r}
# 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,

```{r}
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!

```{r}
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,

```{r fig.width = 7 }
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 *decent* streamflow simulation. **Help**: the order of magnitude of the actual
>parameters is correct (e.g.: the real value of *K0* is 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.


## The HBV model as a function

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,

* [R bloggers](https://www.r-bloggers.com/2016/05/how-to-make-a-function-in-r/)
* [Advanced R](https://adv-r.hadley.nz/functions.html?q=func#functions)

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

```{r definition, echo = TRUE}
## 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")* and *voila*!!!...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](https://www.r-bloggers.com/2016/05/how-to-make-a-function-in-r/). 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?

## Using the function to run the model

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.

```{r fig.width = 7}
# 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. 

## Sensitivity analysis

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. 

```{r}
# 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?

## References