Skip to content
GitLab

Source

Target

Hide whitespace changes
Inline Side-by-side
Showing
with 1747 additions and 0 deletions
+1747 -0
vignettes/V00_airgr_ref.bib 0 → 100644
View file @ 4bccfee2
@phdthesis{ficchi_adaptive_2017,
title = {An adaptive hydrological model for multiple time-steps: diagnostics and improvements based on fluxes consistency},
shorttitle = {An adaptive hydrological model for multiple time-steps},
url = {https://www.theses.fr/2017PA066097},
abstract = {Cette thèse vise à explorer la question du changement d'échelle temporelle en modélisation hydrologique conceptuelle. Les principaux objectifs sont : (i) étudier les effets du changement du pas de temps sur les performances, les paramètres et la structure des modèles hydrologiques ; (ii) mettre au point un modèle pluie-débit applicable à différents pas de temps. Notre point de départ est le modèle global journalier GR4J, développé à Irstea. Ce modèle a été choisi comme le modèle de référence à adapter à d'autres résolutions plus fines, jusqu'à des pas de temps infra-horaires, en suivant une approche descendante. Pour nos tests, nous avons construit une base de données de 240 bassins versants non influencés en France, à différents pas de temps allant de 6 minutes à 1 jour, en utilisant: (i) les données pluviométriques à 6 minutes et la réanalyse des lames d'eau journalières à plus haute résolution spatiale ; (ii) les données de température journalière pour le calcul de l'évapotranspiration potentielle ; (iii) les données hydrométriques à pas de temps variable. Nous avons étudié l'impact de la distribution temporelle des entrées sur les performances du modèle en se focalisant sur la simulation de crue, sur la base de 2400 événements. Ensuite, notre évaluation du modèle a porté sur l'analyse de la cohérence des flux internes du modèle à différents pas de temps, afin d'assurer une performance satisfaisante à travers un fonctionnement du modèle cohérent. Notre diagnostic du modèle nous a permis d'identifier une amélioration de la structure du modèle à différents pas de temps infra-journaliers basée sur la complexification de la composante d'interception du modèle.},
urldate = {2017-11-24},
school = {Université Pierre et Marie Curie, Paris 6},
author = {Ficchi, Andrea},
month = feb,
year = {2017},
keywords = {Rainfall-runoff modelling, airGRcite, Cohérence structurelle, Événements de crue, Interception, Modèles GR, Modélisation pluie-Débit, Pas de temps fin, Short time step}
}
@article{ficchi_hydrological_2019,
title = {Hydrological modelling at multiple sub-daily time steps: model improvement via flux-matching},
issn = {00221694},
shorttitle = {Hydrological modelling at multiple sub-daily time steps},
url = {https://linkinghub.elsevier.com/retrieve/pii/S0022169419305281},
doi = {10.1016/j.jhydrol.2019.05.084},
language = {en},
urldate = {2019-06-11},
journal = {Journal of Hydrology},
author = {Ficchì, Andrea and Perrin, Charles and Andréassian, Vazken},
month = jun,
year = {2019},
keywords = {airGR, airGRcite}
}
@phdthesis{mathevet_quels_2005,
address = {Paris},
title = {Quels modèles pluie-débit globaux au pas de temps horaire ? {Développements} empiriques et comparaison de modèles sur un large échantillon de bassins versants},
url = {http://webgr.irstea.fr/wp-content/uploads/2012/07/2005-MATHEVET-THESE.pdf},
school = {ENGREF},
author = {Mathevet, Thibault},
year = {2005},
keywords = {airGR, IRSTEA-HBAN-HYDRO},
file = {Mathevet_2005_Quels_modèles_pluie-débit_globaux_au_pas_de_temps.pdf:E\:\\Data\\boulot\\misc\\scientific_doc\\_bibliographie\\zotero\\Thèse\\Mathevet_2005_Quels_modèles_pluie-débit_globaux_au_pas_de_temps.pdf:application/pdf}
}
@phdthesis{mouelhi_vers_2003,
title = {Vers une chaîne cohérente de modèles pluie-débit conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et journalier},
url = {http://webgr.irstea.fr/wp-content/uploads/2012/07/2003-MOUELHI-THESE.pdf},
abstract = {,},
urldate = {2016-03-22},
school = {Paris, ENGREF},
author = {Mouelhi, Safouane},
month = jan,
year = {2003},
keywords = {airGR, IRSTEA-HBAN-HYDRO},
file = {Mouelhi_2003_Vers_une_chaîne_cohérente_de_modèles_pluie-débit.pdf:E\:\\Data\\boulot\\misc\\scientific_doc\\_bibliographie\\zotero\\Thèse\\Mouelhi_2003_Vers_une_chaîne_cohérente_de_modèles_pluie-débit.pdf:application/pdf}
}
@phdthesis{le_moine_bassin_2008,
title = {Le bassin versant de surface vu par le souterrain : une voie d’amélioration des performances et du réalisme des modèles pluie-débit ?},
url = {http://webgr.irstea.fr/wp-content/uploads/2012/07/2008-LE_MOINE-THESE.pdf},
school = {Université Pierre et Marie Curie, Paris 6},
author = {Le Moine, Nicolas},
year = {2008},
keywords = {airGR, IRSTEA-HBAN-HYDRO},
file = {Le_Moine_2008_Le_bassin_versant_de_surface_vu_par_le_souterrain.pdf:E\:\\Data\\boulot\\misc\\scientific_doc\\_bibliographie\\zotero\\Thèse\\Le_Moine_2008_Le_bassin_versant_de_surface_vu_par_le_souterrain.pdf:application/pdf}
}
@article{mouelhi_stepwise_2006,
title = {Stepwise development of a two-parameter monthly water balance model},
volume = {318},
issn = {00221694},
url = {http://linkinghub.elsevier.com/retrieve/pii/S0022169405003033},
doi = {10.1016/j.jhydrol.2005.06.014},
language = {en},
number = {1-4},
urldate = {2016-03-22},
journal = {Journal of Hydrology},
author = {Mouelhi, Safouane and Michel, Claude and Perrin, Charles and Andréassian, Vazken},
month = mar,
year = {2006},
keywords = {airGR, IRSTEA-HBAN-HYDRO},
pages = {200--214}
}
@article{mouelhi_linking_2006,
title = {Linking stream flow to rainfall at the annual time step: {The} {Manabe} bucket model revisited},
volume = {328},
issn = {0022-1694},
shorttitle = {Linking stream flow to rainfall at the annual time step},
doi = {10.1016/j.jhydrol.2005.12.022},
abstract = {Trying to determine annual stream flow only from annual rainfall and possibly annual potential evapotranspiration is a challenge that forces hydrologists to focus on very basic questions such as: Is stream flow of the current year only dependent on this year's rainfall? Does previous years' rainfall play a role? Is the role, if any, of antecedent rainfall confined to the transfer function, to the production function or both? How much complexity (i.e. how many free parameters) is required to describe the rainfall-runoff transformation at the annual time step? Our analysis starts with the bucket model of Manabe [Manabe, S., 1969. Climate and the ocean circulation. 1. The atmospheric circulation and the hydrology of the Earth's surface. Mon. Weather Rev. 97(11), 739-774]. We assess this model using data from a large set of basins located in five different countries. We follow a new approach to compare the Manabe model performance to that of other models. The poor relevance of this model at the annual time step is demonstrated, and we put forward a new model that is both simple and clearly better than Manabe bucket model. This model provides new insight on crucial questions raised in the area of annual rainfall-runoff modelling.},
number = {1-2},
urldate = {2011-06-06},
journal = {Journal of Hydrology},
author = {Mouelhi, Safouane and Michel, Claude and Perrin, Charles and Andréassian, Vazken},
month = aug,
year = {2006},
keywords = {IRSTEA-HBAN-HYDRO},
pages = {283--296}
}
@article{perrin_improvement_2003,
title = {Improvement of a parsimonious model for streamflow simulation},
volume = {279},
issn = {0022-1694},
url = {http://www.sciencedirect.com/science/article/B6V6C-49507JG-1/2/c9654f1418557e3106055320abca339a},
doi = {10.1016/S0022-1694(03)00225-7},
abstract = {Hydrologists have been struggling over the past decades to improve rainfall-runoff models. As a consequence, models proposed 20-30 years ago still keep evolving as progress is made in the understanding of catchment hydrological behaviour. Here we present the GR4J model, a daily lumped rainfall-runoff model which is the result of a continuous improvement process over the last 15 years. The article provides the mathematical formulation of a new four-parameter version of the model. Model performance is assessed on a large sample of catchments: compared to other rainfall-runoff models, the GR4J performance is among the best ones. It also gives better results than the previous three-parameter model version, especially in the simulation of low flows. The tests indicate that a four-parameter structure corresponds to the maximum level of complexity that could be afforded in the model. Adding more free parameters did not bring significant improvements. The gain in model robustness with this new version should enhance the confidence in the practical use of this simple model for water engineering and resource management. The discussion underlines the potential limits introduced in the modelling process when one relies on a priori concepts in building a model structure and it stresses the value of large catchment samples to assess models.},
number = {1-4},
urldate = {2010-07-13},
journal = {Journal of Hydrology},
author = {Perrin, Charles and Michel, Claude and Andréassian, Vazken},
year = {2003},
keywords = {airGR, GR4J model, IRSTEA-HBAN-HYDRO, Model complexity, Model improvement, Model robustness, Parsimony, Rainfall-runoff modelling},
pages = {275--289},
file = {Perrin_et_al_2003_Improvement_of_a_parsimonious_model_for.pdf:E\:\\Data\\boulot\\misc\\scientific_doc\\_bibliographie\\zotero\\Article de revue\\Perrin_et_al_2003_Improvement_of_a_parsimonious_model_for.pdf:application/pdf}
}
@article{pushpalatha_downward_2011,
title = {A downward structural sensitivity analysis of hydrological models to improve low-flow simulation},
volume = {411},
issn = {0022-1694},
url = {http://www.sciencedirect.com/science/article/pii/S0022169411006846},
doi = {10.1016/j.jhydrol.2011.09.034},
abstract = {Better simulation and earlier prediction of river low flows are needed for improved water management. Here, a top–down structural analysis to improve a hydrological model in a low-flow simulation perspective is presented. Starting from a simple but efficient rainfall–runoff model (GR5J), we analyse the sensitivity of low-flow simulations to progressive modifications of the model’s structure. These modifications correspond to the introduction of more complex routing schemes and/or the addition of simple representations of groundwater–surface water exchanges. In these tests, we wished to improve low-flow simulation while avoiding performance losses in high-flow conditions, i.e. keeping a general model. In a typical downward modelling perspective, over 60 versions of the model were tested on a large set of French catchments corresponding to various low-flow conditions, and performance was evaluated using criteria emphasising errors in low-flow conditions. The results indicate that several best performing structures yielded quite similar levels of efficiency. The addition of a new flow component to the routing part of the model yielded the most significant improvement. In spite of the close performance of several model structures, we conclude by proposing a modified model version of GR5J with a single additional parameter.},
number = {1–2},
urldate = {2014-04-16},
journal = {Journal of Hydrology},
author = {Pushpalatha, Raji and Perrin, Charles and Le Moine, Nicolas and Mathevet, Thibault and Andréassian, Vazken},
month = dec,
year = {2011},
keywords = {airGR, Downward approach, GR6J model, IRSTEA-HBAN-HYDRO, low flows, Lumped model, Model efficiency, Simulation, Uncertainty},
pages = {66--76},
file = {Pushpalatha_et_al_2011_A_downward_structural_sensitivity_analysis_of.pdf:E\:\\Data\\boulot\\misc\\scientific_doc\\_bibliographie\\zotero\\Article de revue\\Pushpalatha_et_al_2011_A_downward_structural_sensitivity_analysis_of.pdf:application/pdf}
}
@article{valery_as_2014,
title = {'{As} simple as possible but not simpler': what is useful in a temperature-based snow-accounting routine? {Part} 2 - {Sensitivity} analysis of the {Cemaneige} snow accounting routine on 380 catchments},
issn = {0022-1694},
shorttitle = {'{As} simple as possible but not simpler'},
url = {http://www.sciencedirect.com/science/article/pii/S0022169414003321},
doi = {10.1016/j.jhydrol.2014.04.058},
abstract = {This paper investigates the degree of complexity required in a snow accounting routine to ultimately simulate flows at the catchment outlet. We present a simple, parsimonious and general snow accounting routine (SAR), called Cemaneige, that can be associated with any precipitation–runoff model to simulate discharge at the catchment scale. To get results of general applicability, this SAR was tested on a large set of 380 catchments from four countries (France, Switzerland, Sweden and Canada) and combined with four different hydrological models. Our results show that five basic features provide a good reliability and robustness to the SAR, namely considering: (1) a transition range of temperature for the determination of the solid fraction of precipitation; (2) five altitudinal bands of equal area for snow accumulation; (3) the cold-content of the snowpack (with a parameter controlling snowpack inertia); (4) a degree-day factor controlling snowmelt; (5) uneven snow distribution in each band. This general SAR includes two internal states (the snowpack and its cold-content). Results also indicate that only two free parameters (here snowmelt factor and cold-content factor) are warranted in a SAR at the daily time step and that further complexity is not supported by improvements in flow simulation efficiency. To justify the reasons for considering the five features above, a sensitivity analysis comparing Cemaneige with other SAR versions is performed. It analyses the snow processes which should be selected or not to bring significant improvement in model performances. Compared with the six existing SARs presented in the companion article (Valéry et al., this issue) on the 380 catchments set, Cemaneige shows better performance on average than five of these six SARs. It provides performance similar to the sixth SAR (MORD4) but with only half its number of free parameters. However, CemaNeige still appears perfectible on mountainous catchments (France and Switzerland) where the lumped SAR, MORD4, outperforms Cemaneige. Cemaneige can easily be adapted for simulation on ungauged catchments: fixing its two parameters to default values much less degrades performances than the other best performing SAR. This may partly due to the Cemaneige parsimony.},
number = {517(0)},
urldate = {2014-05-14},
journal = {Journal of Hydrology},
author = {Valéry, Audrey and Andréassian, Vazken and Perrin, Charles},
year = {2014},
keywords = {airGR, Degree-day approach, IRSTEA-HBAN-HYDRO, Precipitation-runoff models, Snow accounting Routine (SAR), Snow accumulation, Snowmelt, snowpack variability},
pages = {1176--1187}
}
@article{andreassian_seeking_2014,
title = {Seeking genericity in the selection of parameter sets: {Impact} on hydrological model efficiency},
volume = {50},
issn = {00431397},
shorttitle = {Seeking genericity in the selection of parameter sets},
url = {http://doi.wiley.com/10.1002/2013WR014761},
doi = {10.1002/2013WR014761},
language = {en},
number = {10},
urldate = {2017-10-25},
journal = {Water Resources Research},
author = {Andréassian, Vazken and Bourgin, François and Oudin, Ludovic and Mathevet, Thibault and Perrin, Charles and Lerat, Julien and Coron, Laurent and Berthet, Lionel},
month = oct,
year = {2014},
pages = {8356--8366}
}
@article{riboust_revisiting_2019,
title = {Revisiting a {Simple} {Degree}-{Day} {Model} for {Integrating} {Satellite} {Data}: {Implementation} of {Swe}-{Sca} {Hystereses}},
volume = {67},
issn = {0042-790X},
shorttitle = {Revisiting a {Simple} {Degree}-{Day} {Model} for {Integrating} {Satellite} {Data}},
url = {http://content.sciendo.com/view/journals/johh/67/1/article-p70.xml},
doi = {10.2478/johh-2018-0004},
number = {1},
urldate = {2019-02-18},
journal = {Journal of Hydrology and Hydromechanics},
author = {Riboust, Philippe and Thirel, Guillaume and Moine, Nicolas Le and Ribstein, Pierre},
month = mar,
year = {2019},
keywords = {airGRcite},
pages = {70--81}
}
@article{delavenne_regularization_2019,
ids = {lavenneRegularizationApproachImprove2019a},
title = {A {{Regularization Approach}} to {{Improve}} the {{Sequential Calibration}} of a {{Semidistributed Hydrological Model}}},
author = {de Lavenne, Alban and Andréassian, Vazken and Thirel, Guillaume and Ramos, Maria-Helena and Perrin, Charles},
date = {2019},
journaltitle = {Water Resources Research},
volume = {55},
pages = {8821--8839},
issn = {1944-7973},
doi = {10.1029/2018WR024266},
url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018WR024266},
urldate = {2021-06-07},
abstract = {In semidistributed hydrological modeling, sequential calibration usually refers to the calibration of a model by considering not only the flows observed at the outlet of a catchment but also the different gauging points inside the catchment from upstream to downstream. While sequential calibration aims to optimize the performance at these interior gauged points, we show that it generally fails to improve performance at ungauged points. In this paper, we propose a regularization approach for the sequential calibration of semidistributed hydrological models. It consists in adding a priori information on optimal parameter sets for each modeling unit of the semidistributed model. Calibration iterations are then performed by jointly maximizing simulation performance and minimizing drifts from the a priori parameter sets. The combination of these two sources of information is handled by a parameter k to which the method is quite sensitive. The method is applied to 1,305 catchments in France over 30 years. The leave-one-out validation shows that, at locations considered as ungauged, model simulations are significantly improved (over all the catchments, the median KGE criterion is increased from 0.75 to 0.83 and the first quartile from 0.35 to 0.66), while model performance at gauged points is not significantly impacted by the use of the regularization approach. Small catchments benefit most from this calibration strategy. These performances are, however, very similar to the performances obtained with a lumped model based on similar conceptualization.},
annotation = {\_eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2018WR024266},
file = {C\:\\Users\\david.dorchies\\Zotero\\storage\\QNCGSJXB\\Lavenne et al. - 2019 - A Regularization Approach to Improve the Sequentia.pdf;C\:\\Users\\david.dorchies\\Zotero\\storage\\5DMRIV8Y\\2018WR024266.html},
keywords = {regularization,semidistributed model,stepwise calibration,ungauged basins},
langid = {english},
number = {11},
options = {useprefix=true}
}
vignettes/airGR.Rmd → vignettes/V01_get_started.Rmd
View file @ 4bccfee2
---
title: "airGR -- Overview"
title: "Get Started with airGR"
author: "Guillaume Thirel, Olivier Delaigue, Laurent Coron"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{airGR}
%\VignetteIndexEntry{Get Started with airGR}
%\VignetteEncoding{UTF-8}
---
# Introduction
**airGR** is a package which brings into the [**R software**](https://cran.r-project.org/) the hydrological modelling tools used and developed at the [Catchment Hydrology Team](https://webgr.irstea.fr/?lang=en) at [Irstea (France)](http://www.irstea.fr/en/), including the [**GR rainfall-runoff models**](https://webgr.irstea.fr/modeles/?lang=en) and a snowmelt and accumulation model, **CemaNeige**. Each model core is coded in **Fortran** to ensure low computational time. The other package functions (i.e. mainly the calibration algorithm and the efficiency criteria calculation) are coded in **R**.
**airGR** is a package that brings into the [**R software**](https://cran.r-project.org/) the hydrological modelling tools used and developed at the [Catchment Hydrology Research Group](https://webgr.inrae.fr/home/) at [INRAE (France)](https://www.inrae.fr/en), including the [**GR rainfall-runoff models**](https://webgr.inrae.fr/models/) that can be applied either on a **lumped** or **semi-distributed** way. A snow accumulation and melt model ([**CemaNeige**](https://webgr.inrae.fr/models/snow-model/)) and the associated functions for the calibration and evaluation of models are also included. Each model core is coded in **Fortran** to ensure low computational time. The other package functions (i.e. mainly the calibration algorithm and the efficiency criteria calculation) are coded in **R**.
The **airGR** package has been designed to fulfil two major requirements: to facilitate the use by non-expert users and to allow flexibility regarding the addition of external criteria, models or calibration algorithms. The names of the functions and their arguments were chosen to this end. **airGR** also contains basics plotting facilities.
The **airGR** package has been designed to fulfill two major requirements: to facilitate the use by non-expert users and to allow flexibility regarding the addition of external criteria, models or calibration algorithms. The names of the functions and their arguments were chosen to this end. **airGR** also contains basics plotting facilities.
Six hydrological models and one snowmelt and accumulation model are implemented in **airGR**. The snow model can be used alone or together with the daily hydrological models.
Seven hydrological models and one snow melt and accumulation model are implemented in **airGR**. The hydrological models can be applied either on a lumped way or on a semi-distributed way (on sub-catchments). The snow model can either be used alone or with the daily or hourly hydrological models. Naturally each hydrological model can also be used alone.
The models can be called within **airGR** using the following functions:
* `RunModel_GR4H()`: four-parameter hourly lumped hydrological model [@mathevet_quels_2005]
* `RunModel_GR5H()`: five-parameter hourly lumped hydrological model [@ficchi_adaptive_2017; @ficchi_hydrological_2019]
* `RunModel_GR4J()`: four-parameter daily lumped hydrological model [@perrin_improvement_2003]
* `RunModel_GR5J()`: five-parameter daily lumped hydrological model [@le_moine_bassin_2008]
* `RunModel_GR6J()`: six-parameter daily lumped hydrological model [@pushpalatha_downward_2011]
* `RunModel_GR2M()`: two-parameter monthly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_stepwise_2006]
* `RunModel_GR1A()`: one-parameter yearly lumped hydrological model [@mouelhi_vers_2003; @mouelhi_linking_2006]
* `RunModel_CemaNeige()`: two-parameter degree-day snowmelt and accumulation model [@valery_as_2014]
* `RunModel_CemaNeige()`: two-parameter degree-day snowmelt and accumulation model [@valery_as_2014; @riboust_revisiting_2019]
* `RunModel_CemaNeigeGR4H()`: combined use of **GR4H** and **CemaNeige**
* `RunModel_CemaNeigeGR5H()`: combined use of **GR5H** and **CemaNeige**
* `RunModel_CemaNeigeGR4J()`: combined use of **GR4J** and **CemaNeige**
* `RunModel_CemaNeigeGR5J()`: combined use of **GR5J** and **CemaNeige**
* `RunModel_CemaNeigeGR6J()`: combined use of **GR6J** and **CemaNeige**
The [**GRP**](https://webgr.irstea.fr/modeles/modele-de-prevision-grp/?lang=en) forecasting model and the [**Otamin**](https://webgr.irstea.fr/modeles/otamin/?lang=en) predictive uncertainty tool are not availabe in **airGR**.
The [**GRP**](https://webgr.inrae.fr/models/hydrological-forecasting-model-grp/) forecasting model and the [**Otamin**](https://webgr.inrae.fr/software/otamin/) predictive uncertainty tool are not available in **airGR**.
In this vignette, we show how to prepare and run a calibration and a simulation with airGR hydrological models.
# Loading data
......@@ -42,11 +48,11 @@ In the following example, we use a data sample contained in the package. For rea
First, it is necessary to load the **airGR** package:
```{r}
```{r, warning=FALSE}
library(airGR)
```
This is an example of a `data.frame` of hydrometeorological observations time series for a fictional catchment included in the **airGR** package that contains:
Below is presented an example of a `data.frame` of daily hydrometeorological observations time series for a fictional catchment included in the **airGR** package that contains:
* *DatesR*: dates in the POSIXt format
* *P*: average precipitation [mm/day]
......@@ -57,9 +63,9 @@ This is an example of a `data.frame` of hydrometeorological observations time se
```{r}
data(L0123001)
summary(BasinObs)
summary(BasinObs, digits = 2)
```
The usual functions (e.g. `read.table()`) can be used to load real-case datasets.
The usual functions (e.g. `read.table()`) can be used to load real-case data sets.
......@@ -69,18 +75,18 @@ To run a model, the functions of the **airGR** package (e.g. the models, calibra
To facilitate the use of the package, there are several functions dedicated to the creation of these objects:
* `CreateInputsModel()`: prepares the inputs for the different hydrological models (times series of dates, rainfall, observed streamflow, etc.)
* `CreateRunOptions()`: prepares the options for the hydrological model run (warm-up period, calibration period, etc.)
* `CreateInputsCrit()`: prepares the options in order to compute efficiency criterions (choice of the criterion, choice of the transformation on streamflows: "log", "root", etc.)
* `CreateInputsModel()`: prepares the inputs for the different hydrological models (times series of dates, precipitation, observed discharge, etc.)
* `CreateRunOptions()`: prepares the options for the hydrological model run (warm up period, calibration period, etc.)
* `CreateInputsCrit()`: prepares the options in order to compute the efficiency criterion (choice of the criterion, choice of the transformation on discharge: "log", "sqrt", etc.)
* `CreateCalibOptions()`: prepares the options for the hydrological model calibration algorithm (choice of parameters to optimize, predefined values for uncalibrated parameters, etc.)
## InputsModel object
To run a GR hydrologic model, the user has to prepare the input data with the `CreateInputsModel()` function.
To run a GR hydrological model or CemaNeige, the user has to prepare the input data with the `CreateInputsModel()` function.
As arguments, this function needs the function name corresponding to the model the user wants to run, a vector of dates, a vector of precipitation and a vector of potential evapotranspiration.
In the example below, we already have the potential evapotranspiration. If the user doesn't have these data, it is possible to compute it with the [Oudin's formula](http://dx.doi.org/10.1016/j.jhydrol.2004.08.026) with the `PEdaily_Oudin()` function (this function only needs julian days, daily average air temperature and latitude).
In the example below, we already have the potential evapotranspiration. If the user does not have these data, it is possible to compute it with the [Oudin's formula](http://dx.doi.org/10.1016/j.jhydrol.2004.08.026) with the `PE_Oudin()` function (this function only needs Julian days, daily average air temperature and latitude).
Missing values (`NA`) of precipitation (or potential evapotranspiration) are **not allowed**.
......@@ -92,7 +98,6 @@ str(InputsModel)
```
## RunOptions object
The `CreateRunOptions()` function allows to prepare the options required to the `RunModel*()` functions, which are the actual models functions.
......@@ -100,24 +105,24 @@ The `CreateRunOptions()` function allows to prepare the options required to the
The user must at least define the following arguments:
* `FUN_MOD`: the name of the model function to run
* `InputsModel`: the associated inputs data
* `InputsModel`: the associated input data
* `IndPeriod_Run`: the period on which the model is run
To select a period for which the user wants to run the model, select the corresponding indexes for different time periods (not the POSIXt dates), as follows:
```{r}
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="31/12/1999 00:00"))
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
str(Ind_Run)
```
The intialization of hydrological models is of the utmost importance. Indeed, an inaccurate intialisation causes poor quality streamflow simulations during the earliest stages of the running period. For example, in the GR models, the reservoirs levels are by default set to 50 % of their capacity, which may be far from their ideal value. Two solutions are offered to accurately initialize the GR models in **airGR**: manually predefining the initial states or running the models during a warm up period before the actual running period. It is generally advised to set up this warm up period to be equal or superior to one year.
The initialization of hydrological models is of the utmost importance. Indeed, an inaccurate initialization causes poor quality discharge simulations during the earliest stages of the running period. For example, in the GR models, by default, the production and the routing store levels store level are respectively set to 30 % and 50 % of their capacity, which may be far from their ideal value. Two solutions are offered to accurately initialize the GR models in **airGR**: manually predefining the initial states (e.g. from a previous run) or running the models during a warm up period before the actual running period. It is generally advised to set up this warm up period to be equal or superior to one year.
As a consequence, it is possible to define in `CreateRunOptions()` the following arguments:
* `IniStates`: the initial states of the 2 unit hydrographs (20 + 40 = 60 units)
* `IniResLevels`: the initial levels of the production and routing stores
* `IndPeriod_WarmUp`: the warm-up period used to run the model, to be defined in the same format as `IndPeriod_Run`
* `IndPeriod_WarmUp`: the warm up period used to run the model, to be defined in the same format as `IndPeriod_Run`
```{r}
......@@ -128,31 +133,34 @@ str(RunOptions)
```
The `CreateRunOptions()` function returns warnings if the default initialization options are used:
* `IniStates` and `IniResLevels` are automatically set to initialize all the model states at 0, except for the production and routing stores which are initialised at 50 % of their capacity
* `IndPeriod_WarmUp` default setting ensures a one-year warm-up using the time steps preceding the `IndPeriod_Run`, if available
* `IniStates` and `IniResLevels` are automatically set to initialize all the model states at 0, except for the production and routing stores, which are initialized at respectively 30 % and 50 % of their capacity
* `IndPeriod_WarmUp` default setting ensures a one-year warm up using the time steps preceding the `IndPeriod_Run`, if available
## InputsCrit object
The `CreateInputsCrit()` function allows to prepare the input in order to calculate a criterion. It is possible to define the following arguments:
* `FUN_CRIT`: the name of the error criterion function (they are introduced later on)
* `InputsModel`: the inputs of the hydrological model previously prepared by the `CeateInputsModel()` function
* `FUN_CRIT`: the name of the error criterion function (the available functions are introduced later on)
* `InputsModel`: the inputs of the hydrological model previously prepared by the `CreateInputsModel()` function
* `RunOptions`: the options of the hydrological model previously prepared by the `CreateRunOptions()` function
* `Qobs`: the observed streamflows expressed in *mm/time step*
* `VarObs`: the name of the considered variable (by default `"Q"` for the discharge)
* `Obs`: the observed variable time serie (e.g. the discharge expressed in *mm/time step*)
Missing values (`NA`) are **allowed** for observed streamflows.
Missing values (`NA`) are **allowed** for observed discharge.
It is possible to compute a composite criterion (e.g. the average between NSE computed on discharge and NSE computed on log of discharge). In this case, users have to provide lists to the following arguments (some of the are optional): `FUN_CRIT`, `Obs`, `VarObs`, `BoolCrit`, `transfo`, `Weights.`
```{r}
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions, Qobs = BasinObs$Qmm[Ind_Run])
RunOptions = RunOptions, VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run])
str(InputsCrit)
```
## CalibOptions object
The user needs to prepare the calibration options with the `CreateCalibOptions()`function. For that, it is necessary to define the following arguments:
Before using the automatic calibration tool, the user needs to prepare the calibration options with the `CreateCalibOptions()` function. For that, it is necessary to define the following arguments:
* `FUN_MOD`: the name of the model function
* `FUN_CALIB`: the name of the calibration algorithm
......@@ -165,7 +173,7 @@ str(CalibOptions)
# Criteria
The evaluation of the quality of a simulation is estimated through the calculation of criteria. These criteria can be used both as objective-functions during the calibration of the model, or as a measure for evaluating its control performance.
The evaluation of the quality of a simulation is estimated through the calculation of criteria. These criteria can be used both as objective-functions during the calibration of the model, or as a measure for evaluating its performance on a control period.
The package offers the possibility to use different criteria:
......@@ -175,7 +183,7 @@ The package offers the possibility to use different criteria:
* `ErrorCrit_KGE()`: Kling-Gupta efficiency criterion (KGE)
* `ErrorCrit_KGE2()`: modified Kling-Gupta efficiency criterion (KGE')
It is also possible to create user-defined criteria. For doing that, it is only necessary to define the function in **R** following the same syntax that the criteria functions included in **airGR**.
It is also possible to create user-defined criteria. For doing that, it is only necessary to define the function in **R** following the same syntax as the criteria functions included in **airGR**.
# Calibration
......@@ -193,33 +201,32 @@ The calibration algorithm optimizes the error criterion selected as objective-fu
```{r}
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions,
InputsCrit = InputsCrit, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J, FUN_CRIT = ErrorCrit_NSE)
FUN_MOD = RunModel_GR4J)
Param <- OutputsCalib$ParamFinalR
Param
```
The `Calibration_Michel()` function is the only one implemented in the **airGR** package to calibrate the model, but the user can implement its own calibration function.
The `Calibration_Michel()` function is the only one implemented in the **airGR** package to calibrate the model, but the user can implement its own calibration function. Two vignettes explain how it can be done ([2.1 Plugging in new calibration](V02.1_param_optim.html) and [2.2 MCMC parameter estimation](V02.2_param_mcmc.html)).
This function returns a vector with the parameters of the chosen model, which means that the number of values can differ depending on the model that is used. It is possible to use the `Calibration_Michel()` function with user-implemented hydrological models.
The `Calibration_Michel()` function returns a vector with the parameters of the chosen model, which means that the number of values can differ depending on the model that is used. It is possible to use the `Calibration_Michel()` function with user-implemented hydrological models.
# Validation
# Control
This step assess the predictive capacity of the model. Validation is defined as the estimation of the accuracy of the model on datasets that are not used in its construction, and in particular its calibration.
The classical way to perform a validation is to keep data from a period separated from the calibration period. If possible, this control period should correspond to climatic situations rather that differ from those of the calibration period in order to better point out the qualities and weakness of the model. This exercise is necessary for assessing the robustness of the model, that is to say its ability to keep stable performances outside the calibration conditions.
This step assesses the predictive capacity of the model. Control is defined as the estimation of the accuracy of the model on data sets that are not used in its construction, and in particular its calibration.
The classical way to perform a control is to keep data from a period separated from the calibration period. If possible, this control period should correspond to climatic situations that differ from those of the calibration period in order to better point out the qualities and weaknesses of the model. This exercise is necessary for assessing the robustness of the model, that is to say its ability to keep stable performances outside of the calibration conditions.
Performing a model validation with **airGR** is similar to running a simulation (see below).
Performing a model control with **airGR** is similar to running a simulation (see below), followed by the computation of one or several performance criteria.
# Simulation
## Simulation run
To run a model, the user has to use the `RunModel*()` functions (`InputsModel`, `RunOptions` and parameters).
All the data needed have already been prepared.
All the data needed have already been prepared in the previous steps defined in this guide.
```{r}
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions, Param = Param)
......@@ -227,31 +234,41 @@ str(OutputsModel)
```
## Results preview
Although it is possible for the user to design its own graphics from the outputs of the `RunModel*()` functions, the **airGR** package offers the possibility to make use of the `plot.OutputsModel()` function (or `plotl()` with a `OutputsModel` object). This function returns a dashboard of results including various graphs (depending on the model used):
Although it is possible for the user to design its own graphics from the outputs of the `RunModel*()` functions, the **airGR** package offers the possibility to make use of the `plot()` function. This function returns a dashboard of results including various graphs (depending on the model used):
* time series of total precipitation and simulated streamflows (and observed streamflows if provided)
* interannual median monthly simulated streamflows (and monthly observed streamflows if provided)
* cumulative frequency plot for simulated streamflows (and for observed streamflows if provided)
* correlation plot between simulated and observed streamflows (if observed streamflows provided)
* time series of total precipitation and simulated discharge (and observed discharge if provided)
* interannual average daily simulated discharge (and daily observed discharge if provided) and interannual average monthly precipitation
* cumulative frequency plot for simulated discharge (and for observed discharge if provided)
* correlation plot between simulated and observed discharge (if observed discharge provided)
```{r,eval=F}
plot(OutputsModel, Qobs = BasinObs$Qmm[Ind_Run])
```
Moreover, if the CemaNeige model is used, the temperature and the simulated snowpack time series are plotted.
Moreover, if the CemaNeige model is used, the air temperature and the simulated snowpack water equivalent time series are plotted.
## Efficiency criterion
To evaluate the efficiency of the model, it is possible to use the same criterion as defined at the calibration step or to use an other criterion.
To evaluate the efficiency of the model, it is possible to use the same criterion as defined at the calibration step or to use another criterion.
```{r}
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
str(OutputsCrit)
```
```{r}
OutputsCrit <- ErrorCrit_KGE(InputsCrit = InputsCrit, OutputsModel = OutputsModel)
str(OutputsCrit)
```
## Features diagram
The figure below presents the implementation flow of the **airGR** functions.
<br><br>
<img src="fig/airGR_features_diagram.svg" alt="airGR features diagram" width=75% />
# References
vignettes/V02.1_param_optim.Rmd 0 → 100644
View file @ 4bccfee2
---
title: "Plugging in new calibration algorithms in airGR"
author: "François Bourgin, Guillaume Thirel"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Plugging in new calibration algorithms}
%\VignetteEncoding{UTF-8}
---
```{r setup, warning=FALSE, include=FALSE, fig.keep='none', results='hide'}
library(airGR)
library(DEoptim)
# library(hydroPSO) # Needs R version >= 3.6 or latticeExtra <= 0.6-28 on R 3.5. Archived on 2023-10-16 as requires archived packages 'hydroTSM' and 'hydroGOF'.
library(Rmalschains)
library(caRamel)
library(ggplot2)
library(GGally)
# source("airGR.R")
set.seed(321)
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
load(system.file("vignettesData/vignetteParamOptimCaramel.rda", package = "airGR"))
```
# Introduction
## Scope
The Michel's calibration strategy (`Calibration_Michel()` function) is the calibration algorithm proposed in **airGR**. However, other optimization methods can be used in combination with **airGR**.
We show here how to use different R packages to perform parameter estimation.
In this vignette, we use the **GR4J** model to illustrate the different optimization strategies.
In particular, we assume that the R global environment contains input climate data, observed discharge and functions from the [Get Started](V01_get_started.html) vignette, as shown below.
Please note that the calibration period is defined in the `CreateRunOptions()` function .
<!-- ```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=FALSE} -->
<!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
<!-- ``` -->
```{r Calibration_Michel, echo=TRUE, eval=FALSE}
example("Calibration_Michel")
```
In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
```{r RunOptions, results='hide', eval=FALSE}
RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
Outputs_Sim = "Qsim")
```
Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used.
Please note that this vignette is only for illustration purposes and does not provide any guidance about which optimization strategies is recommended for the family of the **GR** models.
## Definition of the necessary function
Parameter estimation can be performed by defining a function that takes a parameter set as input and returns the value of the performance criterion.
There are two important steps: the transformation of parameters to real space and the computation of the value of the performance criterion.
Here we choose to minimize the root mean square error.
The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model.
```{r OptimGR4J, warning=FALSE, results='hide', eval=FALSE}
OptimGR4J <- function(ParamOptim) {
## Transformation of the parameter set to real space
RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
Direction = "TR")
## Simulation given a parameter set
OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = RawParamOptim)
## Computation of the value of the performance criteria
OutputsCrit <- airGR::ErrorCrit_RMSE(InputsCrit = InputsCrit,
OutputsModel = OutputsModel,
verbose = FALSE)
return(OutputsCrit$CritValue)
}
```
In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space:
```{r boundsGR4J, warning=FALSE, results='hide', eval=FALSE}
lowerGR4J <- rep(-9.99, times = 4)
upperGR4J <- rep(+9.99, times = 4)
```
# Local optimization
We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space:
```{r local1, warning=FALSE, results='hide', eval=FALSE}
startGR4J <- c(4.1, 3.9, -0.9, -8.7)
optPORT <- stats::nlminb(start = startGR4J,
objective = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
control = list(trace = 1))
```
The RMSE value reaches a local minimum value after 35 iterations.
We can also try a multi-start approach to test the consistency of the local optimization.
Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function).
For each starting point, a local optimization is performed.
```{r local2, warning=FALSE, results='hide', eval=FALSE}
startGR4JDistrib <- TransfoParam_GR4J(ParamIn = CalibOptions$StartParamDistrib,
Direction = "RT")
startGR4J <- expand.grid(data.frame(startGR4JDistrib))
optPORT_ <- function(x) {
opt <- stats::nlminb(start = x,
objective = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
control = list(trace = 1))
}
listOptPORT <- apply(startGR4J, MARGIN = 1, FUN = optPORT_)
```
We can then extract the best parameter sets and the value of the performance criteria:
```{r local3, warning=FALSE, results='hide', eval=FALSE}
parPORT <- t(sapply(listOptPORT, function(x) x$par))
objPORT <- sapply(listOptPORT, function(x) x$objective)
resPORT <- data.frame(parPORT, RMSE = objPORT)
```
As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets.
```{r local4, warning=FALSE}
summary(resPORT)
```
The existence of several local minima illustrates the importance of defining an appropriate starting point or of using a multi-start strategy or a global optimization strategy.
# Global optimization
Global optimization is most often used when facing a complex response surface, with multiple local mimina.
Here we use the following R implementation of some popular strategies:
* [DEoptim: differential evolution](https://cran.r-project.org/package=DEoptim)
* [hydroPSO: particle swarm](https://cran.r-project.org/package=hydroPSO)
* [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains)
## Differential Evolution
```{r optDE, warning=FALSE, results='hide', eval=FALSE}
optDE <- DEoptim::DEoptim(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
control = DEoptim::DEoptim.control(NP = 40, trace = 10))
```
## Particle Swarm
```{r hydroPSO1, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
# to install the package temporary removed from CRAN
# Rtools needed (windows : https://cran.r-project.org/bin/windows/Rtools/)
# install.packages("https://cran.r-project.org/src/contrib/Archive/hydroPSO/hydroPSO_0.5-1.tar.gz",
# repos = NULL, type = "source", dependencies = TRUE)
```
```{r hydroPSO2, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
control = list(write2disk = FALSE, verbose = FALSE))
```
## MA-LS-Chains
```{r optMALS, warning=FALSE, results='hide', eval=FALSE}
optMALS <- Rmalschains::malschains(fn = OptimGR4J,
lower = lowerGR4J, upper = upperGR4J,
maxEvals = 2000)
```
# Results
As it can be seen in the table below, the four additional optimization strategies tested lead to very close optima.
```{r resGLOB1, warning=FALSE, echo=FALSE, eval=FALSE}
resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
round(rbind(
OutputsCalib$ParamFinalR,
airGR::TransfoParam_GR4J(ParamIn = optPORT$par , Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = as.numeric(optPSO$par) , Direction = "TR"),
airGR::TransfoParam_GR4J(ParamIn = optMALS$sol , Direction = "TR")),
digits = 3))
```
```{r resGLOB2, warning=FALSE, echo=FALSE}
resGLOB
```
<!-- This is an expected result because the response surface for quadratic performance criteria of the **GR4J** model is generally sufficiently well defined in the transformed parameter space to allow using a local optimization strategy instead of a more time consuming global optimization strategy. -->
# Multiobjective optimization
Multiobjective optimization is used to explore possible trade-offs between different performances criteria.
Here we use the following R implementation of an efficient strategy:
* [caRamel: Automatic Calibration by Evolutionary Multi Objective Algorithm](https://cran.r-project.org/package=caRamel)
Motivated by using the rainfall-runoff model for low flow simulation, we explore the trade-offs between the KGE values obtained without any data transformation and with the inverse transformation.
First, the `OptimGR4J()` function previously used is modified to return two values.
```{r, warning=FALSE, results='hide', eval=FALSE}
InputsCrit_inv <- InputsCrit
InputsCrit_inv$transfo <- "inv"
MOptimGR4J <- function(i) {
if (algo == "caRamel") {
ParamOptim <- x[i, ]
}
## Transformation of the parameter set to real space
RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
Direction = "TR")
## Simulation given a parameter set
OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = RawParamOptim)
## Computation of the value of the performance criteria
OutputsCrit1 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit,
OutputsModel = OutputsModel,
verbose = FALSE)
## Computation of the value of the performance criteria
OutputsCrit2 <- airGR::ErrorCrit_KGE(InputsCrit = InputsCrit_inv,
OutputsModel = OutputsModel,
verbose = FALSE)
return(c(OutputsCrit1$CritValue, OutputsCrit2$CritValue))
}
```
## caRamel
caRamel is a multiobjective evolutionary algorithm combining the MEAS algorithm and the NGSA-II algorithm.
```{r, warning=FALSE, results='hide', eval=FALSE}
algo <- "caRamel"
optMO <- caRamel::caRamel(nobj = 2,
nvar = 4,
minmax = rep(TRUE, 2),
bounds = matrix(c(lowerGR4J, upperGR4J), ncol = 2),
func = MOptimGR4J,
popsize = 100,
archsize = 100,
maxrun = 15000,
prec = rep(1.e-3, 2),
carallel = FALSE,
graph = FALSE)
```
The algorithm returns parameter sets that describe the pareto front, illustrating the trade-off between overall good performance and good performance for low flow.
```{r, fig.width=6, fig.height=6, warning=FALSE}
ggplot() +
geom_point(aes(optMO$objectives[, 1], optMO$objectives[, 2])) +
coord_equal(xlim = c(0.4, 0.9), ylim = c(0.4, 0.9)) +
xlab("KGE") + ylab("KGE [1/Q]") +
theme_bw()
```
The parameter sets can be viewed in the parameter space, illustrating different populations.
```{r fig.height=6, fig.width=6, message=FALSE, warning=FALSE}
param_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
airGR::TransfoParam_GR4J(x, Direction = "TR")
})
GGally::ggpairs(data.frame(t(param_optMO)), diag = NULL) + theme_bw()
```
```{r fig.height=6, fig.width=12, message=FALSE, warning=FALSE}
RunOptions$Outputs_Sim <- "Qsim"
run_optMO <- apply(optMO$parameters, MARGIN = 1, FUN = function(x) {
airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = x)
}$Qsim)
run_optMO <- data.frame(run_optMO)
ggplot() +
geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
y = run_optMO$X1)) +
geom_line(aes(x = as.POSIXct(InputsModel$DatesR[Ind_Run]),
y = run_optMO$X54),
colour = "darkred") +
scale_x_datetime(limits = c(as.POSIXct("1998-01-01"), NA)) +
ylab("Discharge [mm/d]") + xlab("Date") +
scale_y_sqrt() +
theme_bw()
```
vignettes/V02.2_param_mcmc.Rmd 0 → 100644
View file @ 4bccfee2
---
title: "Parameter estimation within a Bayesian MCMC framework"
author: "François Bourgin"
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{MCMC Parameter estimation}
%\VignetteEncoding{UTF-8}
---
```{r, include=FALSE, fig.keep='none', results='hide'}
library(airGR)
library(coda)
library(FME)
library(ggmcmc)
set.seed(123)
load(system.file("vignettesData/vignetteParamMCMC.rda", package = "airGR"))
```
# Introduction
## Scope
In this vignette, we give an example of parameter estimation within a Bayesian MCMC approach.
We use the **GR4J** model and we assume that the R global environment contains data and functions from the **airGR** [Get Started](V01_get_started.html) vignette.
<!-- ```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=FALSE} -->
<!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
<!-- ``` -->
```{r, echo=TRUE, eval=FALSE, eval=FALSE}
example("Calibration_Michel")
```
In order for the `RunModel_*()` functions to run faster during the parameter estimation process, it is recommended that the outputs contain only the simulated flows (see the `Outputs_Sim` argument in the `CreateRunOptions()` help page).
```{r, results='hide', eval=FALSE}
RunOptions <- airGR::CreateRunOptions(FUN_MOD = RunModel_GR4J, InputsModel = InputsModel,
IndPeriod_Run = Ind_Run,
Outputs_Sim = "Qsim")
```
Please refer to the [2.1 Plugging in new calibration](V02.1_param_optim.html) vignette for explanations on how to plug in a parameter estimation algorithm to **airGR**.
Please note that this vignette is only for illustration purposes and does not provide any guidance about which parameter inference strategy is recommended for the family of the GR models.
## Standard Least Squares (SLS) Bayesian inference
We show how to use the DRAM algorithm for SLS Bayesian inference, with the `modMCMC()` function of the [FME](https://cran.r-project.org/package=FME) package.
First, we need to define a function that returns twice the opposite of the log-likelihood for a given parameter set.
Nota: in the `LogLikeGR4J()` function, the computation of the log-likelihood is simplified in order to ensure a good computing performance. It corresponds to a translation of the two following lines.
```{r, echo=TRUE, eval=FALSE, purl=FALSE}
Likelihood <- sum((ObsY - ModY)^2, na.rm = TRUE)^(-sum(!is.na(ObsY)) / 2)
LogLike <- -2 * log(Likelihood)
```
In our simplified setting of Gaussian likelihood with measurement error integrated out, the log of the sum of squared error is related to the log-likelihood.
Note that we do not use here any discharge transformation, although applying Box-Cox transformation is quite popular in hydrological modelling.
```{r, results='hide', eval=FALSE}
LogLikeGR4J <- function(ParamOptim) {
## Transformation to real space
RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
Direction = "TR")
## Simulation given a parameter set
OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
RunOptions = RunOptions,
Param = RawParamOptim)
## Computation of the log-likelihood: N * log(SS)
ObsY <- InputsCrit$Obs
ModY <- OutputsModel$Qsim
LogLike <- sum(!is.na(ObsY)) * log(sum((ObsY - ModY)^2, na.rm = TRUE))
}
```
# MCMC algorithm for Bayesian inference
## Estimation of the best-fit parameters as a starting point
We start by using the PORT optimization routine to estimate the best-fit parameters.
```{r, results='hide', eval=FALSE}
optPORT <- stats::nlminb(start = c(4.1, 3.9, -0.9, -8.7),
objective = LogLikeGR4J,
lower = rep(-9.9, times = 4), upper = rep(9.9, times = 4),
control = list(trace = 1))
iniParPORT <- optPORT$par
```
## Running 3 chains for convergence assessment
We run 3 chains with different initial values to assess the convergence of the Markov chains.
The number of iterations is fixed to 2000 with a burning length of 0.
Nota: in this example, there are relatively few iterations (2000), in order to limit the running time of this vignette. In addition, the burning length has been set to zero in order to show the convergence process but, in a true exercise, it is better to define more iterations (5000) and to burn the first iterations.
With the DRAM algorithm, the covariance of the proposal is updated every 100 runs and delayed rejection is applied.
```{r, results='hide', eval=FALSE}
iniParPORT <- data.frame(Chain1 = iniParPORT,
Chain2 = iniParPORT,
Chain3 = iniParPORT,
row.names = paste0("X", 1:4))
iniParPORT <- sweep(iniParPORT, MARGIN = 2, STATS = c(1, 0.9, 1.1), FUN = "*")
iniParPORT[iniParPORT < -9.9] <- -9.9
iniParPORT[iniParPORT > +9.9] <- +9.9
mcmcDRAM <- apply(iniParPORT, MARGIN = 2, FUN = function(iIniParPORT) {
FME::modMCMC(f = LogLikeGR4J,
p = iIniParPORT,
lower = rep(-9.9, times = 4), ## lower bounds for GR4J
upper = rep(+9.9, times = 4), ## upper bounds for GR4J
niter = 2000,
jump = 0.01,
outputlength = 2000,
burninlength = 0,
updatecov = 100, ## adaptative Metropolis (AM)
ntrydr = 2) ## delayed rejection (RD)
})
```
## MCMC diagnostics and visualisation tools
There are several diagnostics that can be used to check the convergence of the chains.
The R package [coda](https://cran.r-project.org/package=coda) provides several diagnostic tests.
Among others, the Gelman and Rubin's convergence can be used. A value close to 1 suggests acceptable convergence.
The result will be better with more iterations than 2000. As we kept the iterations during the convergence process, we have to set the `autoburnin` argument to `TRUE` in order to consider only the second half of the series.
Note that we rescale model parameter sets of the GR4J model from the transformed space to the real space.
```{r, results='hide', eval=FALSE}
multDRAM <- coda::as.mcmc.list(lapply(mcmcDRAM, FUN = function(x) {
coda::as.mcmc(airGR::TransfoParam_GR4J(as.matrix(x$pars), Direction = "TR"))
}))
gelRub <- coda::gelman.diag(multDRAM, autoburnin = TRUE)$mpsrf
```
```{r}
gelRub
```
In addition, graphical tools can be used, with for example the [ggmcmc](https://cran.r-project.org/package=ggmcmc) package.
First, the evolution of the Markov chains can be seen with a traceplot:
```{r, fig.width=6, fig.height=9, warning=FALSE}
parDRAM <- ggmcmc::ggs(multDRAM) ## to convert object for using by all ggs_* graphical functions
ggmcmc::ggs_traceplot(parDRAM)
```
The posterior density for each parameter can then be visualised:
```{r, fig.width=6, fig.height=9, warning=FALSE}
burnParDRAM <- parDRAM[parDRAM$Iteration > 1000, ] # to keep only the second half of the series
ggmcmc::ggs_density(burnParDRAM)
```
Finally, a paired plot can be useful to look at the correlation between parameters.
Highly correlated parameters can be seen as an indication for a structural deficit of the model used.
```{r, fig.width=6, fig.height=6, results='hide'}
ggmcmc::ggs_pairs(burnParDRAM, lower = list(continuous = "density"))
```
## Exploring further possibilities
We only presented one MCMC algorithm and one parameter inference setting.
Nonetheless, other approaches can be explored within the same framework.
One can for example use different data transformations to deal with the limitations of the Gaussian error model.
<!-- Another extension is to infer the parameters of more advanced error model in addition of the **GR4J** model parameters. -->
vignettes/V03_param_sets_GR4J.Rmd 0 → 100644
View file @ 4bccfee2
---
title: "Generalist parameter sets for the GR4J model"
author: "Olivier Delaigue, Guillaume Thirel"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{GR4J parameter sets}
%\VignetteEncoding{UTF-8}
---
# Introduction
## Scope
```{r, warning=FALSE, include=FALSE}
library(airGR)
options(digits = 3)
```
In the **airGR** package, the classical way to calibrate a model is to use the Michel's algorithm (see the `Calibration_Michel()` function).
The Michel's algorithm combines a global and a local approach. A screening is first performed either based on a rough predefined grid (considering various initial values for each parameter) or from a list of initial parameter sets.
The best set identified in this screening is then used as a starting point for the steepest descent local search algorithm.
In some specific situations, for example if the calibration period is too short and by consequence non representative of the catchment behaviour, a local calibration algorithm can give poor results.
In this vignette, we show a method using a known parameter set that can be used as an alternative for the grid-screening calibration procedure, and we well compare to two methods using the `Calibration_Michel()` function. The generalist parameters sets introduced here are taken from @andreassian_seeking_2014.
## Data preparation
We load an example data set from the package and the GR4J parameter sets.
```{r, warning=FALSE}
## loading catchment data
data(L0123001)
## loading generalist parameter sets
data(Param_Sets_GR4J)
```
The given GR4J **X4u** variable does not correspond to the actual GR4J **X4** parameter. As explained in @andreassian_seeking_2014 [section 2.1], the given GR4J **X4u** value has to be adjusted (rescaled) using catchment area (S) [km2] as follows: `X4 = X4u / 5.995 * S^0.3` (please= note that **the formula is erroneous in the publication**).
It means we need first to transform the **X4** parameter.
```{r, warning=FALSE}
Param_Sets_GR4J$X4 <- Param_Sets_GR4J$X4u / 5.995 * BasinInfo$BasinArea^0.3
Param_Sets_GR4J$X4u <- NULL
Param_Sets_GR4J <- as.matrix(Param_Sets_GR4J)
```
Please, find below the summary of the `r nrow(Param_Sets_GR4J)` sets of the `r ncol(Param_Sets_GR4J)` parameters.
```{r, warning=FALSE, echo=FALSE}
summary(Param_Sets_GR4J)
```
## Object model preparation
We assume that the R global environment contains data and functions from the [Get Started](V01_get_started.html) vignette.
The calibration period has been defined from **1990-01-01** to **1990-02-28**, and the validation period from **1990-03-01** to **1999-12-31**.
As a consequence, in this example the calibration period is very short, less than 6 months.
```{r, warning=FALSE, include=TRUE}
## preparation of the InputsModel object
InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J, 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/01/1990 00:00"),
which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="28/02/1990 00:00"))
## preparation of the RunOptions object for the calibration period
RunOptions_Cal <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Cal)
## efficiency criterion: Nash-Sutcliffe Efficiency
InputsCrit_Cal <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions_Cal, Obs = BasinObs$Qmm[Ind_Cal])
## ---- validation step
## validation period selection
Ind_Val <- seq(which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="01/03/1990 00:00"),
which(format(BasinObs$DatesR, format = "%d/%m/%Y %H:%M")=="31/12/1999 00:00"))
## preparation of the RunOptions object for the validation period
RunOptions_Val <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModel, IndPeriod_Run = Ind_Val)
## efficiency criterion (Nash-Sutcliffe Efficiency) on the validation period
InputsCrit_Val <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel,
RunOptions = RunOptions_Val, Obs = BasinObs$Qmm[Ind_Val])
```
# Calibration of the GR4J model with the generalist parameter sets
It is recommended to use the generalist parameter sets when the calibration period is less than 6 months.
As shown in @andreassian_seeking_2014 [figure 4], a recommended way to use the `Param_Sets_GR4J` `data.frame` is to run the GR4J model with each parameter set and to select the best one according to an objective function (here we use the Nash-Sutcliffe Efficiency criterion).
```{r, warning=FALSE, message=FALSE}
OutputsCrit_Loop <- apply(Param_Sets_GR4J, 1, function(iParam) {
OutputsModel_Cal <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Param = iParam)
OutputsCrit <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
return(OutputsCrit$CritValue)
})
```
Find below the `r nrow(Param_Sets_GR4J)` criteria corresponding to the different parameter sets.
The criterion values are quite low (from `r OutputsCrit_Loop[which.min(OutputsCrit_Loop)]` to `r OutputsCrit_Loop[which.max(OutputsCrit_Loop)]`), which can be expected as this does not represents an actual calibration.
```{r, echo=FALSE}
OutputsCrit_Loop
```
The parameter set corresponding to the best criterion is the following:
```{r, warning=FALSE, message=FALSE, echo=FALSE}
Param_Best <- unlist(Param_Sets_GR4J[which.max(OutputsCrit_Loop), ])
Param_Best
## validation
OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val,
Param = Param_Best)
OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
```
Now we can compute the Nash-Sutcliffe Efficiency criterion on the validation period. A quite good value (`r OutputsCrit_Val$CritValue`) is found.
# Calibration of the GR4J model with the built-in Calibration_Michel function
As seen above, the Michel's calibration algorithm is based on a local search procedure.
It is **not recommanded** to use the `Calibration_Michel()` function when the **calibration period is less than 6 month**.
We will show below its application on the same short period for two configurations of the grid-screening step to demonstrate that it is less efficient than the generalist parameters sets calibration.
## GR4J parameter distributions quantiles used in the grid-screening step
By default, the predefined grid used by the `Calibration_Michel()` function contains parameters quantiles computed after recursive calibrations on 1200 basins (from Australia, France and USA).
```{r, warning=FALSE, message=FALSE}
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
```
```{r, warning=FALSE, message=FALSE, include=FALSE}
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel)
OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR, FUN_MOD = RunModel_GR4J)
OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
## validation
OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val,
Param = OutputsCalib$ParamFinalR)
OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
```
The parameter set corresponding to the best criterion is the following:
```{r, warning=FALSE, message=FALSE, echo=FALSE}
names(OutputsCalib$ParamFinalR) <- paste0("X", 1:4)
OutputsCalib$ParamFinalR
```
The Nash-Sutcliffe Efficiency criterion computed on the calibration period is better (`r OutputsCrit_Cal$CritValue`) than with the generalist parameter sets, but the one computed on the validation period is lower (`r OutputsCrit_Val$CritValue`).
This shows that the generalist parameter sets give more robust model in this case.
## GR4J parameter sets used in the grid-screening step
It is also possible to give to the `CreateCalibOptions()` function a matrix of parameter sets used for the grid-screening calibration procedure.
So, it possible is to use by this way the GR4J generalist parameter sets.
```{r, warning=FALSE, message=FALSE}
CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel,
StartParamList = Param_Sets_GR4J)
```
```{r, warning=FALSE, message=FALSE, include=FALSE}
## calibration
OutputsCalib <- Calibration(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
InputsCrit = InputsCrit_Cal, CalibOptions = CalibOptions,
FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel)
OutputsModel_Cal <- RunModel(InputsModel = InputsModel, RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR, FUN_MOD = RunModel_GR4J)
OutputsCrit_Cal <- ErrorCrit_NSE(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
## validation
OutputsModel_Val <- RunModel_GR4J(InputsModel = InputsModel, RunOptions = RunOptions_Val, Param = OutputsCalib$ParamFinalR)
OutputsCrit_Val <- ErrorCrit_NSE(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
```
Here is the parameter set corresponding to the best criteria found.
```{r, warning=FALSE, message=FALSE, echo=FALSE}
names(OutputsCalib$ParamFinalR) <- paste0("X", 1:4)
OutputsCalib$ParamFinalR
```
The results are the same here. The Nash-Sutcliffe Efficiency criterion computed on the calibration period is better (`r OutputsCrit_Cal$CritValue`), but the one computed on the validation period is just a little bit lower (`r OutputsCrit_Val$CritValue`) than the classical calibration.
Generally, the advantage of using GR4J parameter sets rather than the GR4J generalist parameter quantiles is that they make more sense than a simple exploration of combinations of quantiles of parameter distributions (each parameter set represents a consistent ensemble). In addition, for the first step, the number of iterations is smaller (27 runs instead of 81), which can save time if one wants to run a very large number of calibrations.
# References
vignettes/V04_cemaneige_hysteresis.Rmd 0 → 100644
View file @ 4bccfee2
---
title: "Using satellite snow cover area data for calibrating and improving CemaNeige"
author: "Guillaume Thirel, Olivier Delaigue"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Using CemaNeige with Hysteresis}
%\VignetteEncoding{UTF-8}
---
```{r, warning=FALSE, include=FALSE}
library(airGR)
load(system.file("vignettesData/vignetteCNHysteresis.rda", package = "airGR"))
```
# Introduction
## Scope
Rainfall-runoff models that include a snow accumulation and melt module are still often calibrated using only discharge observations. This can result in poor snow modeling as the swnow module parameters can distorted in order to allow skilful discharge simulation.
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 on the basin and snow melt is more heterogeneous, was implemented.
This new CemaNeige version has 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_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 time) in terms of discharge, which has many implications for climate change impact studies.
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).
In this page, we show how to use and calibrate this new CemaNeige 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, eval=FALSE}
data(X0310010)
summary(BasinObs)
```
## Object model preparation
We assume that the R global environment contains data and functions from the [Get Started](V01_get_started.html) vignette.
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. CemaNeige will be used in coupling with GR4J in this vignette.
```{r, warning=FALSE, eval=FALSE}
## preparation of the InputsModel object
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)
## ---- calibration step
## calibration period selection
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"))
## ---- validation step
## validation period selection
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"))
```
# Calibration and evaluation of the new CemaNeige module
In order to use the Linear Hysteresis, a new argument (`IsHyst`) is added in the `CreateRunOptions()` and `CreateCalibOptions()` functions and has to be set to `TRUE`.
```{r, warning=FALSE, eval=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_revisiting_2019. This is now possible in **airGR** with the added functionality that permits calculating composite criteria by combining different metrics.
```{r, warning=FALSE, eval=FALSE}
## efficiency criterion: 75 % KGE'(Q) + 5 % KGE'(SCA) on each of the 5 layers
InputsCrit_Cal <- CreateInputsCrit(FUN_CRIT = rep("ErrorCrit_KGE2", 6),
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 = rep("ErrorCrit_KGE2", 6),
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, eval=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 period and assess it.
```{r, warning=FALSE, message=FALSE, eval=FALSE}
## run on the calibration period
OutputsModel_Cal <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
RunOptions = RunOptions_Cal,
Param = OutputsCalib$ParamFinalR)
## evaluation
OutputsCrit_Cal <- ErrorCrit(InputsCrit = InputsCrit_Cal, OutputsModel = OutputsModel_Cal)
```
Find below the performance of the model over the calibration period.
```{r, warning=FALSE}
str(OutputsCrit_Cal, max.level = 2)
```
Now we can run the model on the validation period and assess it.
```{r, warning=FALSE, message=FALSE, eval=FALSE}
## run on the validation period
OutputsModel_Val <- RunModel_CemaNeigeGR4J(InputsModel = InputsModel,
RunOptions = RunOptions_Val,
Param = OutputsCalib$ParamFinalR)
## evaluation
OutputsCrit_Val <- ErrorCrit(InputsCrit = InputsCrit_Val, OutputsModel = OutputsModel_Val)
```
Find below the performance of the model over the validation period.
```{r, warning=FALSE}
str(OutputsCrit_Val, max.level = 2)
```
# 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` argument), calibrate and assess it (`InputsCrit` argument). The objective function is only based on KGE'(Q). Note how we set the `IsHyst` argument to `FALSE` in the `CreateRunOptions()` and the `CreateCalibOptions()` functions.
```{r, warning=FALSE, eval=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 = BasinObs$Qmm[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 the model.
```{r, warning=FALSE, eval=FALSE}
## calibration
OutputsCalib_NoHyst <- Calibration(InputsModel = InputsModel,
InputsCrit = InputsCrit_Cal_NoHyst,
RunOptions = RunOptions_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, eval=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_NoHyst$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, message=FALSE, eval=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.
```{r, warning=FALSE}
str(OutputsCrit_Cal_NoHyst, max.level = 2)
str(OutputsCrit_Val_NoHyst, max.level = 2)
```
We can see above that the performance of the initial model is slightly better than 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 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 a linear hysteresis and of SCA data for calibration in CemaNeige.
# References
vignettes/V05_sd_model.Rmd 0 → 100644
View file @ 4bccfee2
---
title: "Simulated vs observed upstream flows in calibration of semi-distributed GR4J model"
author: "David Dorchies"
bibliography: V00_airgr_ref.bib
output: rmarkdown::html_vignette
vignette: >
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Simulated vs observed upstream flows in calibration of semi-distributed GR4J model}
%\VignetteEncoding{UTF-8}
---
```{r, include=FALSE, fig.keep='none', results='hide'}
library(airGR)
options(digits = 3)
```
# Introduction
## Scope
The **airGR** package implements semi-distributed model capabilities using a lag model between subcatchments. It allows to chain together several lumped models as well as integrating anthropogenic influence such as reservoirs or withdrawals.
Here we explain how to implement the semi-distribution with **airGR**. For everyday use, however, it is easier to use the [**airGRiwrm**](https://cran.r-project.org/package=airGRiwrm) package.
`RunModel_Lag` documentation gives an example of simulating the influence of a reservoir in a lumped model. Try `example(RunModel_Lag)` to get it.
In this vignette, we show how to calibrate 2 sub-catchments in series with a semi-distributed model consisting of 2 GR4J models.
For doing this we compare 3 strategies for calibrating the downstream subcatchment:
- using upstream observed flows
- using upstream simulated flows
- using upstream simulated flows and parameter regularisation [@delavenne_regularization_2019]
We finally compare these calibrations with a theoretical set of parameters.
This comparison is based on the Kling-Gupta Efficiency computed on the root-squared discharges as performance criterion.
## Model description
```{r, warning=FALSE, include=FALSE}
library(airGR)
options(digits = 3)
```
We use an example data set from the package that unfortunately contains data for only one catchment.
```{r, warning=FALSE}
## loading catchment data
data(L0123001)
```
Let's imagine that this catchment of 360 km² is divided into 2 subcatchments:
- An upstream subcatchment of 180 km²
- 100 km downstream another subcatchment of 180 km²
We consider that meteorological data are homogeneous on the whole catchment, so we use the same pluviometry `BasinObs$P` and the same evapotranspiration `BasinObs$E` for the 2 subcatchments.
For the observed flow at the downstream outlet, we generate it with the assumption that the upstream flow arrives at downstream with a constant delay of 2 days.
```{r}
QObsDown <- (BasinObs$Qmm + c(0, 0, BasinObs$Qmm[1:(length(BasinObs$Qmm)-2)])) / 2
options(digits = 5)
summary(cbind(QObsUp = BasinObs$Qmm, QObsDown))
options(digits = 3)
```
With a delay of 2 days between the 2 gauging stations, the theoretical Velocity parameter should be equal to:
```{r}
Velocity <- 100 * 1e3 / (2 * 86400)
paste("Velocity: ", format(Velocity), "m/s")
```
# Calibration of the upstream subcatchment
The operations are exactly the same as the ones for a GR4J lumped model. So we do exactly the same operations as in the [Get Started](V01_get_started.html) vignette.
```{r}
InputsModelUp <- CreateInputsModel(FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E)
Ind_Run <- seq(which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1990-01-01"),
which(format(BasinObs$DatesR, format = "%Y-%m-%d") == "1999-12-31"))
RunOptionsUp <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModelUp,
IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
IniStates = NULL, IniResLevels = NULL)
# Error criterion is KGE computed on the root-squared discharges
InputsCritUp <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelUp,
RunOptions = RunOptionsUp,
VarObs = "Q", Obs = BasinObs$Qmm[Ind_Run],
transfo = "sqrt")
CalibOptionsUp <- CreateCalibOptions(FUN_MOD = RunModel_GR4J, FUN_CALIB = Calibration_Michel)
OutputsCalibUp <- Calibration_Michel(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
InputsCrit = InputsCritUp, CalibOptions = CalibOptionsUp,
FUN_MOD = RunModel_GR4J)
```
And see the result of the simulation:
```{r}
OutputsModelUp <- RunModel_GR4J(InputsModel = InputsModelUp, RunOptions = RunOptionsUp,
Param = OutputsCalibUp$ParamFinalR)
```
# Calibration of the downstream subcatchment
## Creation of the InputsModel objects
We need to create `InputsModel` objects completed with upstream information with upstream observed flow for the calibration of first case and upstream simulated flows for the other cases:
```{r}
InputsModelDown1 <- CreateInputsModel(
FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
Qupstream = matrix(BasinObs$Qmm, ncol = 1), # upstream observed flow
LengthHydro = 100, # distance between upstream catchment outlet & the downstream one [km]
BasinAreas = c(180, 180) # upstream and downstream areas [km²]
)
```
For using upstream simulated flows, we should concatenate a vector with the simulated flows for the entire period of simulation (warm-up + run):
```{r}
Qsim_upstream <- rep(NA, length(BasinObs$DatesR))
# Simulated flow during warm-up period (365 days before run period)
Qsim_upstream[Ind_Run[seq_len(365)] - 365] <- OutputsModelUp$RunOptions$WarmUpQsim
# Simulated flow during run period
Qsim_upstream[Ind_Run] <- OutputsModelUp$Qsim
InputsModelDown2 <- CreateInputsModel(
FUN_MOD = RunModel_GR4J, DatesR = BasinObs$DatesR,
Precip = BasinObs$P, PotEvap = BasinObs$E,
Qupstream = matrix(Qsim_upstream, ncol = 1), # upstream observed flow
LengthHydro = 100, # distance between upstream catchment outlet & the downstream one [km]
BasinAreas = c(180, 180) # upstream and downstream areas [km²]
)
```
## Calibration with upstream flow observations
We calibrate the combination of Lag model for upstream flow transfer and GR4J model for the runoff of the downstream subcatchment:
```{r}
RunOptionsDown <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
InputsModel = InputsModelDown1,
IndPeriod_WarmUp = NULL, IndPeriod_Run = Ind_Run,
IniStates = NULL, IniResLevels = NULL)
InputsCritDown <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE, InputsModel = InputsModelDown1,
RunOptions = RunOptionsDown,
VarObs = "Q", Obs = QObsDown[Ind_Run],
transfo = "sqrt")
CalibOptionsDown <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
FUN_CALIB = Calibration_Michel,
IsSD = TRUE) # specify that it's a SD model
OutputsCalibDown1 <- Calibration_Michel(InputsModel = InputsModelDown1,
RunOptions = RunOptionsDown,
InputsCrit = InputsCritDown,
CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
```
`RunModel` is run in order to automatically combine GR4J and Lag models.
```{r}
OutputsModelDown1 <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = OutputsCalibDown1$ParamFinalR,
FUN_MOD = RunModel_GR4J)
```
Performance of the model validation is then:
```{r}
KGE_down1 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown1)
```
## Calibration with upstream simulated flow
We calibrate the model with the `InputsModel` object previously created for substituting the observed upstream flow with the simulated one:
```{r}
OutputsCalibDown2 <- Calibration_Michel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
InputsCrit = InputsCritDown,
CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
ParamDown2 <- OutputsCalibDown2$ParamFinalR
```
## Calibration with upstream simulated flow and parameter regularisation
The regularisation follow the method proposed by @delavenne_regularization_2019.
As a priori parameter set, we use the calibrated parameter set of the upstream catchment and the theoretical velocity:
```{r}
ParamDownTheo <- c(Velocity, OutputsCalibUp$ParamFinalR)
```
The Lavenne criterion is initialised with the a priori parameter set and the value of the KGE of the upstream basin.
```{r}
IC_Lavenne <- CreateInputsCrit_Lavenne(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Obs = QObsDown[Ind_Run],
AprParamR = ParamDownTheo,
AprCrit = OutputsCalibUp$CritFinal)
```
The Lavenne criterion is used instead of the KGE for calibration with regularisation
```{r}
OutputsCalibDown3 <- Calibration_Michel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
InputsCrit = IC_Lavenne,
CalibOptions = CalibOptionsDown,
FUN_MOD = RunModel_GR4J)
```
The KGE is then calculated for performance comparisons:
```{r}
OutputsModelDown3 <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = OutputsCalibDown3$ParamFinalR,
FUN_MOD = RunModel_GR4J)
KGE_down3 <- ErrorCrit_KGE(InputsCritDown, OutputsModelDown3)
```
# Discussion
## Identification of Velocity parameter
Both calibrations overestimate this parameter:
```{r}
mVelocity <- matrix(c(Velocity,
OutputsCalibDown1$ParamFinalR[1],
OutputsCalibDown2$ParamFinalR[1],
OutputsCalibDown3$ParamFinalR[1]),
ncol = 1,
dimnames = list(c("theoretical",
"calibrated with observed upstream flow",
"calibrated with simulated upstream flow",
"calibrated with sim upstream flow and regularisation"),
c("Velocity parameter")))
knitr::kable(mVelocity)
```
## Value of the performance criteria with theoretical calibration
Theoretically, the parameters of the downstream GR4J model should be the same as the upstream one with the velocity as extra parameter:
```{r}
OutputsModelDownTheo <- RunModel(InputsModel = InputsModelDown2,
RunOptions = RunOptionsDown,
Param = ParamDownTheo,
FUN_MOD = RunModel_GR4J)
KGE_downTheo <- ErrorCrit_KGE(InputsCritDown, OutputsModelDownTheo)
```
## Parameters and performance of each subcatchment for all calibrations
```{r}
comp <- matrix(c(0, OutputsCalibUp$ParamFinalR,
rep(OutputsCalibDown1$ParamFinalR, 2),
OutputsCalibDown2$ParamFinalR,
OutputsCalibDown3$ParamFinalR,
ParamDownTheo),
ncol = 5, byrow = TRUE)
comp <- cbind(comp, c(OutputsCalibUp$CritFinal,
OutputsCalibDown1$CritFinal,
KGE_down1$CritValue,
OutputsCalibDown2$CritFinal,
KGE_down3$CritValue,
KGE_downTheo$CritValue))
colnames(comp) <- c("Velocity", paste0("X", 1:4), "KGE(√Q)")
rownames(comp) <- c("Calibration of the upstream subcatchment",
"Calibration 1 with observed upstream flow",
"Validation 1 with simulated upstream flow",
"Calibration 2 with simulated upstream flow",
"Calibration 3 with simulated upstream flow and regularisation",
"Validation theoretical set of parameters")
knitr::kable(comp)
```
Even if calibration with observed upstream flows gives an improved performance criteria, in validation using simulated upstream flows the result is quite similar as the performance obtained with the calibration with upstream simulated flows. The theoretical set of parameters give also an equivalent performance but still underperforming the calibration 2 one. Regularisation allows to get similar performance as the one for calibration with simulated flows but with the big advantage of having parameters closer to the theoretical ones (Especially for the velocity parameter).
# References
vignettes/fig/airGR_features_diagram.svg 0 → 100644
View file @ 4bccfee2
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg xmlns:dc="http://purl.org/dc/elements/1.1/" xmlns:cc="http://creativecommons.org/ns#" xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" xmlns:svg="http://www.w3.org/2000/svg" xmlns="http://www.w3.org/2000/svg" xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd" xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape" id="svg2" version="1.1" inkscape:version="0.91 r13725" xml:space="preserve" width="900" height="675" viewBox="0 0 900 675" sodipodi:docname="airGR_features.svg"><metadata id="metadata8"><rdf:RDF><cc:Work rdf:about=""><dc:format>image/svg+xml</dc:format><dc:type rdf:resource="http://purl.org/dc/dcmitype/StillImage"/><dc:title/></cc:Work></rdf:RDF></metadata><defs id="defs6"><clipPath clipPathUnits="userSpaceOnUse" id="clipPath16"><path d="m 0,1.2207e-4 720,0 0,539.99999793 -720,0 L 0,1.2207e-4 Z" id="path18" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath26"><path d="M 3.08e-7,1.2207e-4 720,1.2207e-4 720,540.00012 3.08e-7,540.00012 l 0,-539.99999793 z" id="path28" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath102"><path d="M -2.645e-6,540 720,540 720,1.2207e-4 -6.582e-6,1.2207e-4" id="path104" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath116"><path d="M -1.161e-6,540 720,540 720,1.2207e-4 l -720.000001161,0" id="path118" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath132"><path d="M -2.728e-6,540 720,540 720,1.2207e-4 l -720.000006665,0" id="path134" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath><clipPath clipPathUnits="userSpaceOnUse" id="clipPath178"><path d="M -1.52e-6,540 720,540 720,6.104e-6 -1.52e-6,6.104e-6" id="path180" inkscape:connector-curvature="0" style="clip-rule:evenodd"/></clipPath></defs><sodipodi:namedview pagecolor="#ffffff" bordercolor="#666666" borderopacity="1" objecttolerance="10" gridtolerance="10" guidetolerance="10" inkscape:pageopacity="0" inkscape:pageshadow="2" inkscape:window-width="1920" inkscape:window-height="1132" id="namedview4" showgrid="false" inkscape:zoom="0.98890193" inkscape:cx="333.64663" inkscape:cy="425.84952" inkscape:window-x="-8" inkscape:window-y="-8" inkscape:window-maximized="1" inkscape:current-layer="g10"/><g id="g10" inkscape:groupmode="layer" inkscape:label="airGR_features" transform="matrix(1.25,0,0,-1.25,0,675)"><g id="g12"><g id="g14" clip-path="url(#clipPath16)"><path d="m 0,6.104e-6 720,0 L 720,540.00001 l -720,0 L 0,6.104e-6 Z" style="fill:#ffffff;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path20" inkscape:connector-curvature="0"/></g></g><g id="g22"><g id="g24" clip-path="url(#clipPath26)"><path d="M 391.5,540 390.85,-2.1562" style="fill:none;stroke:#7f7f7f;stroke-width:2;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:8, 6;stroke-dashoffset:0;stroke-opacity:1" id="path30" inkscape:connector-curvature="0"/></g></g><path d="m 308.97,304.02 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.65 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.97 -6.62,6.62 l 0,26.46 z" style="fill:#00b0f0;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path32" inkscape:connector-curvature="0"/><g id="g34"><path d="m 308.97,304.02 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.65 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.97 -6.62,6.62 l 0,26.46 z" style="fill:none;stroke:#00b0f0;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path36" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,320.57,284.74)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text38"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 52.416 61.866001 71.351997 80.802002 86.832001 93.870003 109.26 118.746 128.196 137.214" y="0" sodipodi:role="line" id="tspan40">CreateInputsModel</tspan></text>
<path d="m 481.47,473.17 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.61,-6.61 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.61 l 0,26.46 z" style="fill:#cc00cc;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path42" inkscape:connector-curvature="0"/><g id="g44"><path d="m 481.47,473.17 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.61,-6.61 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.61 l 0,26.46 z" style="fill:none;stroke:#cc00cc;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path46" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,493.2,453.94)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text48"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 57.438 66.096001 70.181999 74.267998 83.718002 95.634003 104.994 111.024 115.074 124.56 134.00999" y="0" sodipodi:role="line" id="tspan50">CreateCalibOptions</tspan></text>
<path d="m 481.5,400.38 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.62 -6.62,-6.62 l -151.2,0 c -3.65,0 -6.61,2.96 -6.61,6.62 l 0,26.45 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path52" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,494.81,381.12)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text54"><tspan x="0 9.5939999 18.216 22.356001 26.406 35.855999 41.759998 50.273998 56.304001 60.354 69.839996 79.290001 88.290001 103.68 107.73 115.29 124.74 133.758" y="0" sodipodi:role="line" id="tspan56">Calibration_Michel</tspan></text>
<path d="m 308.97,131.92 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.618 -6.62,-6.618 l -151.2,0 c -3.65,0 -6.61,2.958 -6.61,6.618 l 0,26.45 z" style="fill:#92d050;fill-opacity:0.30196001;fill-rule:evenodd;stroke:none" id="path58" inkscape:connector-curvature="0"/><g id="g60"><path d="m 308.97,131.92 c 0,3.66 2.96,6.62 6.61,6.62 l 151.2,0 c 3.66,0 6.62,-2.96 6.62,-6.62 l 0,-26.45 c 0,-3.66 -2.96,-6.618 -6.62,-6.618 l -151.2,0 c -3.65,0 -6.61,2.958 -6.61,6.618 l 0,26.45 z" style="fill:none;stroke:#92d050;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path62" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,331.25,112.61)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text64"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 52.416 61.866001 71.351997 80.802002 86.832001 93.870003 103.464 109.746 113.832" y="0" sodipodi:role="line" id="tspan66">CreateInputsCrit</tspan></text>
<path d="m 31.145,132.03 c 0,3.66 2.961,6.62 6.615,6.62 l 49.139,0 c 3.653,0 6.615,-2.96 6.615,-6.62 l 0,-26.46 c 0,-3.65 -2.962,-6.61 -6.615,-6.61 l -49.139,0 c -3.654,0 -6.615,2.96 -6.615,6.61 l 0,26.46 z" style="fill:#92d050;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path68" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,47.784,112.73)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text70"><tspan x="0 9.4499998 13.59 23.076" y="0" sodipodi:role="line" id="tspan72">plot</tspan></text>
<path d="m 481.5,385.53 -23.9,0 c -0.89,0 -1.62,0.73 -1.62,1.62 l 0,0.03 1.62,-1.63 -15.77,0 0,3.25 15.77,0 c 0.9,0 1.63,-0.72 1.63,-1.62 l 0,-0.03 -1.63,1.63 23.9,0 0,-3.25 z m -38.04,-3.23 -9.75,4.88 9.75,4.87 0,-9.75 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path74" inkscape:connector-curvature="0"/><path d="m 308.97,120.32 -14.17,0 c -0.9,0 -1.63,-0.73 -1.63,-1.63 l 0,-0.09 1.63,1.62 -6.05,0 0,-3.25 6.05,0 c 0.89,0 1.62,0.73 1.62,1.63 l 0,0.09 -1.62,-1.62 14.17,0 0,3.25 z m -18.6,3.15 -9.75,-4.87 9.75,-4.88 0,9.75 z" style="fill:#92d050;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path76" inkscape:connector-curvature="0"/><path d="m 120.22,388.8 -57.891,0 c -0.897,0 -1.625,-0.72 -1.625,-1.62 l 0,-240.41 3.25,0 0,240.41 -1.625,-1.63 57.891,0 0,3.25 z m -62.766,-240.4 4.875,-9.75 4.875,9.75 -9.75,0 z" style="fill:#7030a0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path78" inkscape:connector-curvature="0"/><path d="m 308.97,289.16 -77.96,0 c -0.9,0 -1.62,0.73 -1.62,1.63 l 0,68.39 3.25,0 0,-68.39 -1.63,1.62 77.96,0 0,-3.25 z m -82.83,68.4 4.9,9.74 4.85,-9.76 -9.75,0.02 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path80" inkscape:connector-curvature="0"/><path d="m 308.97,220.52 -78.19,0 c -0.89,0 -1.62,0.72 -1.62,1.62 l 0,137.07 3.25,0 0,-137.07 -1.63,1.63 78.19,0 0,-3.25 z m -83.06,137.06 4.87,9.75 4.88,-9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path82" inkscape:connector-curvature="0"/><path d="m 473.4,220.52 90.31,0 c 0.9,0 1.63,0.72 1.63,1.62 l 0,137.05 -3.25,0 0,-137.05 1.62,1.63 -90.31,0 0,-3.25 z m 95.19,137.04 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path84" inkscape:connector-curvature="0"/><path d="m 365.67,51.027 -303.341,0 c -0.897,0 -1.625,0.728 -1.625,1.625 l 0,38.183 3.25,0 0,-38.183 -1.625,1.625 303.341,0 0,-3.25 z M 57.454,89.21 l 4.875,9.75 4.875,-9.75 -9.75,0 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path86" inkscape:connector-curvature="0"/><path d="m 473.4,289.16 90.31,0 c 0.9,0 1.63,0.73 1.63,1.63 l 0,68.4 -3.25,0 0,-68.4 1.62,1.62 -90.31,0 0,-3.25 z m 95.19,68.4 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path88" inkscape:connector-curvature="0"/><path d="m 342.99,387.18 9.92,19.84 70.88,0 9.92,-19.84 -9.92,-19.85 -70.88,0 -9.92,19.85 z" style="fill:#00508c;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path90" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,365.11,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text92"><tspan x="0 8.8739996 17.496 23.4 32.021999" y="0" sodipodi:role="line" id="tspan94">Param</tspan></text>
<path d="m 564.12,173.36 -172.94,0 c -0.89,0 -1.62,-0.73 -1.62,-1.63 l 0,-25.07 3.25,0 0,25.07 -1.63,-1.62 172.94,0 0,3.25 z m -177.81,-25.07 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path96" inkscape:connector-curvature="0"/><g id="g98"><g id="g100" clip-path="url(#clipPath102)"><path d="m 562.06,440.1 0,-16.55 c 0,-0.9 0.72,-1.63 1.62,-1.63 l 0.03,0 -1.62,1.63 0,-8.43 3.25,0 0,8.43 c 0,0.9 -0.73,1.62 -1.63,1.62 l -0.03,0 1.63,-1.62 0,16.55 -3.25,0 z m -3.22,-23.35 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#cc00cc;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path106" inkscape:connector-curvature="0"/></g></g><path d="m 473.4,117.07 124.74,0 c 0.89,0 1.62,0.73 1.62,1.62 l 0,240.49 -3.25,0 0,-240.49 1.63,1.63 -124.74,0 0,-3.25 z m 129.61,240.48 -4.87,9.75 -4.88,-9.75 9.75,0 z" style="fill:#92d050;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path108" inkscape:connector-curvature="0"/><path d="m 354.33,388.8 -34.84,0 c -0.9,0 -1.63,-0.72 -1.63,-1.62 l 0,0 1.63,1.62 -26.72,0 0,-3.25 26.72,0 c 0.9,0 1.62,0.73 1.62,1.63 l 0,0 -1.62,-1.63 34.84,0 0,3.25 z m -59.93,3.25 -9.75,-4.87 9.75,-4.88 0,9.75 z" style="fill:#00508c;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path110" inkscape:connector-curvature="0"/><g id="g112"><g id="g114" clip-path="url(#clipPath116)"><path d="m 200.03,367.33 0,-114.44 c 0,-0.9 -0.73,-1.63 -1.62,-1.63 l 0,0 1.62,1.63 0,-106.32 -3.25,0 0,106.32 c 0,0.89 0.73,1.62 1.63,1.62 l 0,0 -1.63,-1.62 0,114.44 3.25,0 z m 3.25,-219.14 -4.87,-9.75 -4.88,9.75 9.75,0 z" style="fill:#7030a0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path120" inkscape:connector-curvature="0"/></g></g><path d="m 120.22,400.41 c 0,3.65 2.96,6.61 6.62,6.61 l 151.19,0 c 3.66,0 6.62,-2.96 6.62,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.62 -6.62,-6.62 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.62 l 0,26.46 z" style="fill:#7030a0;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path122" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,164.42,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text124"><tspan x="0 9.7740002 19.224001 28.674 44.063999 53.549999 63 71.963997" y="0" sodipodi:role="line" id="tspan126">RunModel</tspan></text>
<g id="g128"><g id="g130" clip-path="url(#clipPath132)"><path d="m 389.56,270.94 0,-14.47 c 0,-0.9 0.73,-1.63 1.62,-1.63 l 0,0 -1.62,1.63 0,-6.36 3.25,0 0,6.36 c 0,0.89 -0.73,1.62 -1.63,1.62 l 0,0 1.63,-1.62 0,14.47 -3.25,0 z m -3.25,-19.2 4.87,-9.75 4.88,9.75 -9.75,0 z" style="fill:#00b0f0;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path136" inkscape:connector-curvature="0"/></g></g><text transform="matrix(1,0,0,-1,509.9,509.71)" style="font-variant:normal;font-weight:bold;font-size:24px;font-family:Calibri;-inkscape-font-specification:'Calibri Bold';writing-mode:lr-tb;fill:#5f5f5f;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text138"><tspan x="0 12.696 24.552 30.455999 36.360001 49.248001 57.240002 68.879997 77.208 83.112 96.024002" y="0" sodipodi:role="line" id="tspan140">Calibration</tspan></text>
<text transform="matrix(1,0,0,-1,142.75,509.38)" style="font-variant:normal;font-weight:bold;font-size:24px;font-family:Calibri;-inkscape-font-specification:'Calibri Bold';writing-mode:lr-tb;fill:#5f5f5f;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text142"><tspan x="0 11.352 17.256001 36.84 49.728001 55.632 67.199997 75.528 81.431999 94.344002" y="0" sodipodi:role="line" id="tspan144">Simulation</tspan></text>
<path d="m 308.97,235.37 c 0,3.65 2.96,6.62 6.62,6.62 l 151.19,0 c 3.66,0 6.62,-2.97 6.62,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.62,-6.61 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.61 l 0,26.46 z" style="fill:#00b0f0;fill-opacity:0.2;fill-rule:evenodd;stroke:none" id="path146" inkscape:connector-curvature="0"/><g id="g148"><path d="m 308.97,235.37 c 0,3.65 2.96,6.62 6.62,6.62 l 151.19,0 c 3.66,0 6.62,-2.97 6.62,-6.62 l 0,-26.46 c 0,-3.65 -2.96,-6.61 -6.62,-6.61 l -151.19,0 c -3.66,0 -6.62,2.96 -6.62,6.61 l 0,26.46 z" style="fill:none;stroke:#00b0f0;stroke-width:3;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:9, 3;stroke-dashoffset:0;stroke-opacity:1" id="path150" inkscape:connector-curvature="0"/></g><text transform="matrix(1,0,0,-1,324.29,216.07)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#000000;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text152"><tspan x="0 9.5939999 15.606 24.57 33.119999 38.880001 47.844002 57.618 67.068001 76.517998 88.433998 97.793999 103.824 107.874 117.36 126.81" y="0" sodipodi:role="line" id="tspan154">CreateRunOptions</tspan></text>
<path d="m 116.19,400.41 c 0,3.65 2.96,6.61 6.62,6.61 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.62 -6.61,-6.62 l -151.2,0 c -3.66,0 -6.62,2.96 -6.62,6.62 l 0,26.46 z" style="fill:#7030a0;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path156" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,160.39,381.14)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text158"><tspan x="0 9.7740002 19.224001 28.674 44.063999 53.549999 63 71.963997" y="0" sodipodi:role="line" id="tspan160">RunModel</tspan></text>
<path d="m 116.19,131.83 c 0,3.65 2.96,6.61 6.62,6.61 l 151.2,0 c 3.65,0 6.61,-2.96 6.61,-6.61 l 0,-26.46 c 0,-3.66 -2.96,-6.618 -6.61,-6.618 l -151.2,0 c -3.66,0 -6.62,2.958 -6.62,6.618 l 0,26.46 z" style="fill:#92d050;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path162" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,167.11,112.51)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text164"><tspan x="0 8.7840004 14.994 20.988001 30.474001 36.756001 46.349998 52.542 56.627998" y="0" sodipodi:role="line" id="tspan166">ErrorCrit</tspan></text>
<path d="m 351.85,52.652 9.92,19.845 59.53,0 9.93,-19.845 -9.93,-19.844 -59.53,0 -9.92,19.844 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:evenodd;stroke:none" id="path168" inkscape:connector-curvature="0"/><text transform="matrix(1,0,0,-1,372.5,46.56)" style="font-variant:normal;font-weight:normal;font-size:18px;font-family:Calibri;-inkscape-font-specification:Calibri;writing-mode:lr-tb;fill:#ffffff;fill-opacity:1;fill-rule:nonzero;stroke:none" id="text170"><tspan x="0 12.114 21.6 30.959999" y="0" sodipodi:role="line" id="tspan172">Qobs</tspan></text>
<g id="g174"><g id="g176" clip-path="url(#clipPath178)"><path d="m 392.81,71.77 0,13.541 c 0,0.897 -0.73,1.625 -1.63,1.625 l 0,0 1.63,-1.625 0,5.416 -3.25,0 0,-5.416 c 0,-0.897 0.73,-1.625 1.62,-1.625 l 0,0 -1.62,1.625 0,-13.541 3.25,0 z m 3.25,17.332 -4.88,9.75 -4.87,-9.75 9.75,0 z" style="fill:#e46c0a;fill-opacity:1;fill-rule:nonzero;stroke:none" id="path182" inkscape:connector-curvature="0"/></g></g><path d="m 563.71,222.75 0.62,-50.29" style="fill:none;stroke:#00b0f0;stroke-width:3.25;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:10;stroke-dasharray:none;stroke-opacity:1" id="path184" inkscape:connector-curvature="0"/></g></svg>
\ No newline at end of file