---
title: "Severn_05: Modeling ungauged stations"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Severn_05: Modeling ungauged stations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: airGRiwrm.bib
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.asp = 0.68,
  out.width = "70%",
  fig.align = "center"
)
```

```{r setup}
library(airGRiwrm)
```

An *Ungauged* station is a virtual hydrometric station where no observed flows are
available for calibration.

## Why modeling *Ungauged* station in a semi-distributed model?

*Ungauged* nodes in a semi-distributed model can be used to reach two different goals:

- increase spatial resolution of the rain fall to improve streamflow simulation [@lobligeoisWhenDoesHigher2014].
- simulate streamflows in location of interest for management purpose

This vignette introduces the implementation in airGRiwrm of the method developed by @lobligeoisMieuxConnaitreDistribution2014 for calibrating *Ungauged* nodes in a
semi-distributed model.

## Presentation of the study case

Using the study case of the vignette #1 and #2, we consider this time that nodes `54001` and
`54029` are *Ungauged*. We simulate the streamflow at these locations by sharing
hydrological parameters of the gauged node `54032`.

```{r network, echo = FALSE}
plot.mermaid("
graph LR
id95[54095]
id01[54001]
id29[54029]

id95 -->| 42 km| id01

subgraph Shared parameters from node 54032
id01 -->| 45 km| 54032
id29 -->| 32 km| 54032
end

54032 -->| 15 km| 54057
54002 -->| 43 km| 54057

classDef UpUng fill:#eef
classDef UpGau fill:#aaf
classDef IntUng fill:#efe
classDef IntGau fill:#afa
classDef DirInj fill:#faa

class id29 UpUng
class 54057,54032 IntGau
class id01 IntUng
class id95,54002 UpGau
")
```

Hydrological parameters at the ungauged nodes will be the same as the one at the
gauged node `54032` except for the unit hydrograph parameter which depend on the
area of the sub-basin. @lobligeoisMieuxConnaitreDistribution2014 provides the
following conversion formula for this parameter:

$$
x_{4i} = \left( \dfrac{S_i}{S_{BV}} \right) ^ {0.3} X_4
$$
With $X_4$ the unit hydrograph parameter for the entire basin at `54032` which
as an area of $S_{BV}$; $S_i$ the area and $x_{4i}$ the parameter for the
sub-basin $i$.

## Using *Ungauged* stations in the airGRiwrm model

*Ungauged* stations are specified by using the model `"Ungauged"` in the `model`
column provided in the `CreateGRiwrm` function:

```{r griwrm}
data(Severn)
nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
nodes$model <- "RunModel_GR4J"
nodes$model[nodes$gauge_id %in% c("54029", "54001")] <- "Ungauged"
griwrmV05 <- CreateGRiwrm(
  nodes,
  list(id = "gauge_id", down = "downstream_id", length = "distance_downstream")
)
griwrmV05
```

It should be noted that the `GRiwrm` object includes a column which automatically
define the first gauged station at downstream for each *Ungauged* node.
It is also possible to manually define the donor node of an *Ungauged* node,
which may be upstream or in a parallel sub-basin. Type `?CreateGRiwrm` for more
details.

On the following network scheme, the *Ungauged* nodes are clearer than gauged ones
with the same color (blue for upstream nodes and green for intermediate and
downstream nodes)

```{r plot_network}
plot(griwrmV05)
```


## Generation of the GRiwrmInputsModel object

The formatting of the input data is described in the vignette "V01_Structure_SD_model". The following code chunk resumes this formatting procedure:

```{r obs}
BasinsObs <- Severn$BasinsObs
DatesR <- BasinsObs[[1]]$DatesR
PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
Precip <- ConvertMeteoSD(griwrmV05, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrmV05, PotEvapTot)
```

Then, the `GRiwrmInputsModel` object can be generated taking into account the new `GRiwrm` object:

```{r InputsModel}
IM_U <- CreateInputsModel(griwrmV05, DatesR, Precip, PotEvap)
```
## Calibration of the model integrating ungauged nodes

Calibration options is detailed in vignette "V02_Calibration_SD_model".
We also apply a parameter regularization here but only where an upstream simulated catchment is available.

The following code chunk resumes this procedure:

```{r RunOptions}
IndPeriod_Run <- seq(
  which(DatesR == (DatesR[1] + 365*24*60*60)), # Set aside warm-up period
  length(DatesR) # Until the end of the time series
)
IndPeriod_WarmUp = seq(1,IndPeriod_Run[1]-1)
RunOptions <- CreateRunOptions(IM_U,
                               IndPeriod_WarmUp = IndPeriod_WarmUp,
                               IndPeriod_Run = IndPeriod_Run)
Qobs <- cbind(sapply(BasinsObs, function(x) {x$discharge_spec}))
InputsCrit <- CreateInputsCrit(IM_U,
                               FUN_CRIT = ErrorCrit_KGE2,
                               RunOptions = RunOptions, Obs = Qobs[IndPeriod_Run,],
                               AprioriIds = c("54057" = "54032", "54032" = "54095"),
                               transfo = "sqrt", k = 0.15
)
CalibOptions <- CreateCalibOptions(IM_U)
```

The **airGR** calibration process is applied on each hydrological node of the `GRiwrm` network from upstream nodes to downstream nodes but this time the calibration of the sub-basin `54032` invokes a semi-distributed model composed of the nodes `54029`, `54001` and `54032` sharing the same parameters.

```{r Calibration}
OC_U <- suppressWarnings(
  Calibration(IM_U, RunOptions, InputsCrit, CalibOptions))
```

Hydrological parameters for sub-basins

## Run the model with the optimized model parameters

The hydrological model uses parameters inherited from the calibration of the gauged sub-basin `54032` for the *Ungauged* nodes `54001` and `54029`:

```{r param}
ParamV05 <- sapply(griwrmV05$id, function(x) {OC_U[[x]]$Param})
dfParam <- do.call(
  rbind,
  lapply(ParamV05, function(x)
    if (length(x)==4) {return(c(NA, x))} else return(x))
)
colnames(dfParam) <- c("velocity", paste0("X", 1:4))
knitr::kable(round(dfParam, 3))
```

We can run the model with these calibrated parameters:

```{r RunModel}
OutputsModels <- RunModel(
  IM_U,
  RunOptions = RunOptions,
  Param = ParamV05
)
```

and plot the comparison of the modeled and the observed flows including the so-called *Ungauged* stations :

```{r plot, fig.height = 5, fig.width = 8}
plot(OutputsModels, Qobs = Qobs[IndPeriod_Run,], which = c("Regime", "CumFreq"))
```


# References