V04_cemaneige_hysteresis.Rmd 10.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
---
title: "Using satellite snow data for calibrating an improved version of CemaNeige"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Improved version of CemaNeige}
  %\VignetteEncoding{UTF-8}
---



# Introduction

## Scope

Rainfall-runoff models that include a snow accumulation and melt module are still often calibrated using only discharge observations.
  
After the work of Riboust et al. (2019), we propose now in airGR an improved version of the degree-day CemaNeige snow and accumulation module. This new version is based on a more accurate representation of the relationship that exists at the basin scale between the Snow Water Equivalent (SWE) and the Snow Cover Area (SCA). To do so, a linear SWE-SCA hysteresis, which represents the fact that snow accumulation is rather homogeneous and snow melt is more heterogeneous, was implemented.
  
This new CemaNeige version presents two more parameters to calibrate. It also presents the advantage of allowing using satellite snow data to constrain the calibration in addition to discharge. 
Riboust et al. (2019) show that while the simulated discharge is not significantly improved, the snow simulation is much improved. In addition, they show that the model is more robust (i.e. transferable) in terms of discharge, which has many implications for climate change impact studies.
  
The configuration that was identified as optimal by Riboust et al. (2019) includes a CemaNeige module run on 5 elevation bands and an objective function determine by a composite function of KGE' calculated on discharge (75% weight) and KGE' calculated on each elevation band (5 % for each).
  
In this page, we show how to use and calibrate this new CameNeige version. 
  
  
## Data preparation
  
We load an example data set from the package. Please note that this data set includes MODIS data that was pre-calculated for 5 elevation bands and for which days with few data (more than 40% cloud coverage) were assigned as missing values. 
  
  
## loading catchment data
```{r, warning=FALSE}
data(X0310010)
summary(BasinObs)
```

## Object model preparation

We assume that the R global environment contains data and functions from the Get Started page.

The calibration period has been defined from 2000-09-01 to 2005-08-31, and the validation period from 2005-09-01 to 2010-07-31.


```{r, warning=FALSE}
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR4J, DatesR = BasinObs$DatesR, 
                                 Precip = BasinObs$P, PotEvap = BasinObs$E)

## ---- calibration step

## short calibration period selection (< 6 months)
Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/09/2000 00:00"), 
               which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="31/08/2005 00:00"))


# ---- validation step

## validation period selection
Ind_Val <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/09/2005 00:00"), 
               which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="31/07/2010 00:00"))
```



# Calibration and evaluation of the new CemaNeige module

In order to use the hysteresis, a new argument (IsHyst) is added in the CreateRunOptions and CreateCalibOptions functions and has to be set to TRUE. 

```{r, warning=FALSE}
## preparation of the RunOptions object for the calibration period
RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Cal, IsHyst = TRUE)

## preparation of the RunOptions object for the validation period
RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
                                   InputsModel = InputsModel, IndPeriod_Run = Ind_Val, IsHyst = TRUE)

## preparation of the CalibOptions object
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J, FUN_CALIB = Calibration_Michel, 
                                   IsHyst = TRUE)
```

In order to calibrate and assess the model performance, we will follow the recommendations of Riboust et al. (2019). This is now possible in airGR with the added functionality that permits calculated composite criteria by combining different metrics. 


```{r, warning=FALSE}
## efficiency criteria: 75 % KGE'(Q) + 5 % KGE'(SCA) on each of the 5 layers
InputsCrit_Cal  <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2),
                                    InputsModel = InputsModel, RunOptions = RunOptions_Cal,
                                    obs = list(BasinObs$Qmm[Ind_Cal], BasinObs$SCA1[Ind_Cal], BasinObs$SCA2[Ind_Cal],
                                               BasinObs$SCA3[Ind_Cal], BasinObs$SCA4[Ind_Cal], BasinObs$SCA5[Ind_Cal]),
                                    varObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
                                    weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))

InputsCrit_Val  <- CreateInputsCrit(FUN_CRIT = list(ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2, ErrorCrit_KGE2),
                                    InputsModel = InputsModel, RunOptions = RunOptions_Val,
                                    obs = list(BasinObs$Qmm[Ind_Val], BasinObs$SCA1[Ind_Val], BasinObs$SCA2[Ind_Val],
                                               BasinObs$SCA3[Ind_Val], BasinObs$SCA4[Ind_Val], BasinObs$SCA5[Ind_Val]),
                                    varObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
                                    weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))
``` 

We can now calibrate the model.

```{r, warning=FALSE}
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
                            InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
                            FUN_MOD = RunModel_CemaNeigeGR4J,
                            FUN_CALIB = Calibration_Michel)
```


Now we can run it on the calibration and validation periods and assess it. 

```{r, warning=FALSE}
## run on the calibration period
OutputsModel_Val <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel, RunOptions = RunOptions, 
                                           Param = OutputsCalib$ParamFinalR)

## evaluation 
OutputsCrit_Cal <- ErrorCrit(InputsCrit = InputsCrit, OutputsModel = OutputsModel_Cal)
```

Find below the performance of the model over the calibration period. 

```{r, warning=FALSE}
str(OutputsCrit_Cal)



## run on the validation period
OutputsModel_Val <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel, RunOptions = RunOptions, 
                                           Param = OutputsCalib$ParamFinalR)

## evaluation 
OutputsCrit_Val <- ErrorCrit(InputsCrit = InputsCrit, OutputsModel = OutputsModel_Val)
```

Find below the performance of the model over the validation period. 

```{r, warning=FALSE}
str(OutputsCrit_Val)
```


# Comparison with the performance of the initial CemaNeige version 

Here we use the same InputsModel object and calibration and validation periods. However, we have to redefine the way we run the model (RunOptions), calibrate and assess it (InputsCrit). The objective function is only based on KGE'(Q). 

```{r, warning=FALSE}
## preparation of RunOptions object
RunOptions_Cal_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
                                          InputsModel = InputsModel, IndPeriod_Run = Ind_Cal, IsHyst = FALSE)

RunOptions_Val_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J, InputsModel = InputsModel, 
                                          IndPeriod_Run = Ind_Val, IsHyst = FALSE)

InputsCrit_Cal_NoHyst <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
                                          InputsModel = InputsModel, RunOptions = RunOptions_Cal_NoHyst,
                                          obs = BasinData$Q[Ind_Cal], varObs = "Q")

## preparation of CalibOptions object
CalibOptions_NoHyst <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
                                          FUN_CALIB = Calibration_Michel, IsHyst = FALSE)
```

We can now calibrate it. 

```{r, warning=FALSE}
## calibration
OutputsCalib_NoHyst <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions,
                                   InputsCrit = InputsCrit_Cal_NoHyst, CalibOptions = CalibOptions_NoHyst,
                                   FUN_MOD = RunModel_CemaNeigeGR4J,
                                   FUN_CALIB = Calibration_Michel)
```

And run it over the calibration and validation periods. 

```{r, warning=FALSE}
OutputsModel_Cal_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel, RunOptions = RunOptions_Cal_NoHyst,  
                                                  Param = OutputsCalib_NoHyst$ParamFinalR)

OutputsModel_Val_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val_NoHyst, 
                                                  Param = OutputsCalib$ParamFinalR)
```

In order to assess the model performance over the two periods, we will use the InputsCrit objects prepared before, which allow assessing also the performance in terms of snow simulation.


```{r, warning=FALSE}
OutputsCrit_Cal_NoHyst <- ErrorCrit(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal_NoHyst)

OutputsCrit_Val_NoHyst <- ErrorCrit(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val_NoHyst)
```

We can see below that the performance of the initial model is similar to the new one over the calibration period in terms of discharge, but also that the new version calibrated using SCA provides better performance in terms of snow. 
However, over the validation, we see that the discharge simulated by the new version brings better performance (in addition to improved SCA also). This shows the interests of the combined use of an hysteresis and of SCA data for calibration in CemaNeige. 


Reference

Riboust, P., Thirel, G., Le Moine, N., and Ribstein, P.: Revisiting a simple degree-day model for integrating satellite data: implementation of SWE-SCA hystereses. Journal of Hydrology and Hydromechanics, DOI: 10.2478/johh-2018-0004, 67, 1, 70–81, 2019.