V04_cemaneige_hysteresis.Rmd 11.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
---
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}
---


12
13
14
15
```{r, echo=FALSE}
library(airGR)
```

16
17
18
19
20
21
22

# Introduction

## Scope

Rainfall-runoff models that include a snow accumulation and melt module are still often calibrated using only discharge observations.
  
23
After the work of @riboust_revisiting_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.
24
25
  
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. 
26
@riboust_revisiting_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.
27
  
28
The configuration that was identified as optimal by @riboust_revisiting_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).
29
30
31
32
33
34
  
In this page, we show how to use and calibrate this new CameNeige version. 
  
  
## Data preparation
  
35
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. 
36
37
38
39
40
41
42
43
  
  
## loading catchment data
```{r, warning=FALSE}
data(X0310010)
summary(BasinObs)
```

44

45
46
## Object model preparation

47
We assume that the R global environment contains data and functions from the [Get Started](V01_get_started.html) vignette.
48
49
50
51
52
53

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
54
55
56
57
58
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_CemaNeigeGR4J,
                                 DatesR = BasinObs$DatesR, Precip = BasinObs$P,
                                 PotEvap = BasinObs$E, TempMean = BasinObs$T,
                                 ZInputs = median(BasinInfo$HypsoData),
                                 HypsoData = BasinInfo$HypsoData, NLayers = 5)
59
60
61
62

## ---- calibration step

## short calibration period selection (< 6 months)
63
64
Ind_Cal <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2000-09-01"), 
               which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2005-08-31"))
65
66
67
68
69


# ---- validation step

## validation period selection
70
71
Ind_Val <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2005-09-01"), 
               which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "2010-07-31"))
72
73
74
75
76
77
```



# Calibration and evaluation of the new CemaNeige module

78
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`. 
79
80
81
82

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

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

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

97
In order to calibrate and assess the model performance, we will follow the recommendations of @riboust_revisiting_2019. This is now possible in **airGR** with the added functionality that permits calculated composite criteria by combining different metrics. 
98
99
100


```{r, warning=FALSE}
101
## efficiency criterion: 75 % KGE'(Q) + 5 % KGE'(SCA) on each of the 5 layers
102
InputsCrit_Cal  <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
103
                                    InputsModel = InputsModel, RunOptions = RunOptions_Cal,
104
                                    Obs = list(BasinObs$Qmm[Ind_Cal],
105
106
107
108
109
                                               BasinObs$SCA1[Ind_Cal],
                                               BasinObs$SCA2[Ind_Cal],
                                               BasinObs$SCA3[Ind_Cal],
                                               BasinObs$SCA4[Ind_Cal],
                                               BasinObs$SCA5[Ind_Cal]),
110
111
                                    VarObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
                                    Weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))
112

113
InputsCrit_Val  <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
114
                                    InputsModel = InputsModel, RunOptions = RunOptions_Val,
115
                                    Obs = list(BasinObs$Qmm[Ind_Val],
116
117
118
119
120
                                               BasinObs$SCA1[Ind_Val],
                                               BasinObs$SCA2[Ind_Val],
                                               BasinObs$SCA3[Ind_Val],
                                               BasinObs$SCA4[Ind_Val],
                                               BasinObs$SCA5[Ind_Val]),
121
122
                                    VarObs = list("Q", "SCA", "SCA", "SCA", "SCA", "SCA"),
                                    Weights = list(0.75, 0.05, 0.05, 0.05, 0.05, 0.05))
123
124
125
126
127
128
129
130
131
132
133
134
135
``` 

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


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

138
```{r, warning=FALSE, message=FALSE}
139
## run on the calibration period
140
141
OutputsModel_Cal <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
                                           RunOptions = RunOptions_Cal, 
142
143
144
                                           Param = OutputsCalib$ParamFinalR)

## evaluation 
145
OutputsCrit_Cal <- ErrorCrit(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
146
147
```

148

149
150
151
Find below the performance of the model over the calibration period. 

```{r, warning=FALSE}
152
153
str(OutputsCrit_Cal, max.level = 2)
```
154

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

157
```{r, warning=FALSE}
158
## run on the validation period
159
160
OutputsModel_Val <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
                                           RunOptions = RunOptions_Val, 
161
162
163
                                           Param = OutputsCalib$ParamFinalR)

## evaluation 
164
OutputsCrit_Val <- ErrorCrit(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
165
166
167
168
169
```

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

```{r, warning=FALSE}
170
str(OutputsCrit_Val, max.level = 2)
171
172
173
```


174

175
176
# Comparison with the performance of the initial CemaNeige version 

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

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

186
187
188
189
RunOptions_Val_NoHyst <- CreateRunOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
                                          InputsModel = InputsModel, 
                                          IndPeriod_Run = Ind_Val,
                                          IsHyst = FALSE)
190
191

InputsCrit_Cal_NoHyst <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
192
193
                                          InputsModel = InputsModel,
                                          RunOptions = RunOptions_Cal_NoHyst,
194
                                          Obs = BasinObs$Qmm[Ind_Cal], VarObs = "Q")
195
196
197

## preparation of CalibOptions object
CalibOptions_NoHyst <- CreateCalibOptions(FUN_MOD = RunModel_CemaNeigeGR4J,
198
199
                                          FUN_CALIB = Calibration_Michel,
                                          IsHyst = FALSE)
200
201
202
203
204
205
```

We can now calibrate it. 

```{r, warning=FALSE}
## calibration
206
207
208
OutputsCalib_NoHyst <- Calibration(InputsModel = InputsModel, InputsCrit = InputsCrit_Cal_NoHyst,
                                   RunOptions = RunOptions_Cal_NoHyst, CalibOptions = CalibOptions_NoHyst,
                                   FUN_MOD = RunModel_CemaNeigeGR4J, FUN_CALIB = Calibration_Michel)
209
210
211
212
213
```

And run it over the calibration and validation periods. 

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

218
219
220
OutputsModel_Val_NoHyst <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
                                                  RunOptions = RunOptions_Val_NoHyst, 
                                                  Param = OutputsCalib_NoHyst$ParamFinalR)
221
222
223
224
225
```

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.


226
227
228
229
230
231
232
233
234
```{r, warning=FALSE, message=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 check the performance over the calibration and the validation period.
235

236
237
238
```{r, warning=FALSE}
str(OutputsCrit_Cal_NoHyst, max.level = 2)
str(OutputsCrit_Val_NoHyst, max.level = 2)
239
240
241
```

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. 
242
However, over the validation period, 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. 
243
244
245



246
# References