The HOOPLA toolbox: a HydrOlOgical Prediction LAboratory to explore ensemble rainfall-runoff modeling

This technical report introduces the HydrOlOgical Prediction LAboratory (HOOPLA) developed at Université Laval for ensemble lumped hydrological modelling. HOOPLA includes functionalities to perform calibration, simulation, and forecast for multiple hydrological models and various time steps. It includes a range of hydrometeorological tools such as calibration algorithms, data assimilation techniques, potential evapotranspiration formulas and a snow accounting routine. HOOPLA is a flexible framework coded in MATLAB that allows easy integration of user-defined hydrometeorological tools. This report 5 also illustrates HOOPLA’s functionalities using a set of 31 Canadian catchments.


Introduction
In the 70s, meteorologists realized the importance of addressing uncertainty in forecasting to achieve more accurate predictions, but it is only in the 90s that ensemble meteorological forecasting was more broadly developed to account for the chaotic nature of the atmospheric processes (Molteni et al., 1996;Buizza et al., 1999). Following on the success, hydrological forecasters 10 adopted the concept of ensemble and gradually integrated probabilistic meteorological forecasts into hydrological prediction systems (e.g. Pappenberger et al., 2005;Thirel et al., 2008;Cloke and Pappenberger, 2009;He et al., 2009). Yet, it appears that consideration of meteorological forcing uncertainty remains insufficient to adequately describe the hydrological uncertainty and that additional (hydrologic) sources of errors need to be taken into account (Vrugt and Robinson, 2007;Salamon and Feyen, 2010). 15 In this respect, hydrological model initial conditions received increasing attention by, in particular, the means of data assimilation. By updating model states according to streamflow measurements, various data assimilation schemes proved to be able to increase forecast performance (Weerts and El Serafy, 2006;DeChant and Moradkhani, 2012;Lee et al., 2011, e.g.). Yet, Salamon and Feyen (2009) pointed out that data assimilation underestimates predictive uncertainty if the model structural error is not explicitly accounted for, which remains a challenge in traditional data assimilation schemes. 20 In parallel, multimodel and model intercomparison experiments indicated that it is unlikely that a single structure may correctly describe rainfall-runoff uncertainty (Clark et al., 2008;Pechlivanidis et al., 2011) and that the combination of several models is often superior to individual ones, in addition to providing a proxy for assessing related uncertainty (Ajami et al., 2006; the procedure to implement new tools to HOOPLA. Thanks to these attributes, HOOPLA is polyvalent and may be used for various purposes and to address different specialties in hydrology. It can be, for example, used for short-range forecasting, data assimilation, or multimodel studies. Because HOOPLA was not conceived for a specific application, it only includes basic functions for result visualization in order not to orient users toward a certain interpretation of the results. As a counterpart, the format of the input and outputs are heavily 75 documented in the user manual. Finally, to address the different sources of uncertainty located along the modelling hydrometeorological chain, HOOPLA adopts an approach based on ensembles. Because ensembles can be computationally expensive, a parallel computing feature that allows to run simultaneously several models and catchments has been implemented (see Appendix B).

80
A collection of 20 lumped and conceptual hydrological models is provided in HOOPLA. These models perform rainfall-runoff transformation by the means of interconnected reservoirs (buckets) that fill and empty over time to simulate the hydrological processes. Additionally, the models assume that a catchment has homogeneous properties over its domain. Models were selected based on their: -Complexity: only lowly to moderately complex models were selected providing that they perform well. Therefore, the 85 retained models possess a rather low number of free parameters and reservoirs (see Table 1).
-Diversity: the models were developed by different research and engineering teams in various contexts and for different purposes. Differences in model structure is therefore fostered and should be reflected in dissimilar behaviors.
The selection originates from the work of Perrin (2000), where 35 models were tested in a context of model comparison to investigate the influence of the number of optimized parameters and model structure on the streamflow simulation performance.
For example, these changes may include a modification of the catchment spatial representation (all distributed models were converted to lumped models), a reduction of the number of free parameters, or discarding elements of the structure that did not address rainfall-runoff modelling (like potential evapotranspiration and snow accounting routine). Table 1 provides a qualitative list of modifications is provided in the HOOPLA user manual (Thiboult et al., 2018).

Snow accounting routine
The CemaNeige (Valéry et al., 2014) snow accounting routine (SAR) is implemented in HOOPLA to simulate snow accumulation and melt. It performs the partitioning of precipitation between liquid and solid phases by using a simple and parsimonious approach based on 1) a spatial discretization of the catchment in altitudinal bands, 2) a transition temperature range that con-105 trols the liquid/solid fraction, 3) a term that describes the snowpack thermal inertia, and 4) a degree-day formula to estimate the snowmelt. CemaNeige results are subsequently used as inputs for the hydrological models, if necessary. CemaNeige is currently the only snow accounting routine included in HOOPLA as other tested routines did not bring satisfactory performance at this time.

Data assimilation schemes 110
Data assimilation aims at integrating information contained in observations to improve simulation accuracy while accounting for uncertainties in the measurements and model. More precisely, data assimilation schemes typically adjust the model states based on the latest available observations to provide better initial conditions for the next modelling step.
Two probabilistic, sequential data assimilation schemes are provided in HOOPLA: -The Sequential Importance Resampling filter (SIR). It belongs to the class of particle filters (PF) and is also referred to 115 as Sequential Monte Carlo or Bootstrap filter. More details about Bayesian inference and the PF derivation can be found in Arulampalam et al. (2002). The latter was the main inspiration for the SIR code implemented in HOOPLA.
-The Ensemble Kalman Filter (EnKF). A sound description of the EnKF is presented in Evensen (2003) and the EnKF code in HOOPLA was written according to Mandel (2006) recommendations.
These two sequential data assimilation schemes were selected to be included in HOOPLA because they have been extensively 120 tested in hydrological sciences and proved to be efficient, which is reflected by their frequent use in the literature (Liu and Gupta, 2007). Even if both the PF and the EnKF adopt a Monte Carlo estimation method, they use a different approach and resort to different approximations.

Calibration algorithms
Calibration aims at identifying an optimal parameter set with regards to a chosen cost function (or efficiency criterion).

125
HOOPLA provides two calibration algorithms: the Shuffle Complex Evolution (SCE, Duan et al., 1992) and the Dynamically Dimensioned Search algorithm (DDS, Tolson and Shoemaker, 2007). Both methods are iterative and global (they seek to find the optimal set among the entire parameter space). They mainly differ in the way the parameter space is sampled.
SCE is based on several groups of points (parameter sets called complexes) spanning the parameter space. At each iteration, these complexes evolve along several mechanisms (namely reflection, contraction, and mutation) and are subsequently recom-130 5 https://doi.org/10.5194/hess-2020-6 Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License.
bined. These steps are iterated until the algorithm converges or a preset maximum of iterations is reached. The SCE is often considered as a state-of-the-art algorithm for high dimensional problems but remains computationally expensive.
DDS is inspired by the manual calibration procedure. The DDS explores the entire parameter space during the first optimization stages, and then progressively fixes parameter values, reducing the parameter space dimension. The parameter space is sampled via random perturbations and reflections. A distinctive feature of the DDS is that it automatically scales the parameter 135 search to find good solutions within a specified maximum number of iterations. As a result, it is a computationally efficient method that provides a better solution than the SCE when the number of model evaluations is limited. Therefore DDS is well suited for calibrating models with a large number of free parameters or that require important computing resources.

Potential evapotranspiration formulas
Potential evapotranspiration is used to quantify the atmospheric demand in water vapour and constitutes one of the hydrological 140 model inputs (along precipitation). A large number of potential evapotranspiration formula exist. Only three were implemented in HOOPLA (see Table 2) because hydrological models usually express little sensitivity towards the PET input (Seiller and Anctil, 2016). Kharrufa and Oudin are parsimonious formulas in the sense that they both rely on an energy balance rationale that only requires mean air temperature as input while the Penman formula is more complex and uses a combinational approach (energy 145 balance and aerodynamic) that necessitates more data.

An application of HOOPLA
This section presents typical results obtained with HOOPLA in calibration, simulation, and forecast mode. Computations are performed with a 3-hour time step on 31 catchments situated in the southern part of the Province of Québec, whose rivers flow into the Saint Lawrence river. The catchments area and the mean annual specific discharge range from 515 to 6839 km 2 and 150 from 1.1 to 3.4 mm/day respectively. The streamflow data were collected by the Direction de l'Expertise Hydrique du Québec Because the watersheds undergo a nival influence, CemaNeige, the snow accounting routine included in HOOPLA, is activated. The Oudin formula (Oudin et al., 2005) is chosen among the available potential evapotranspiration formulas because of its low requirement of input data and known good performance when combined with the 20 aforementioned hydrological models at those sites. Finally, all 20 hydrological models provided with HOOPLA are used.

160
In this example, the 20 hydrological models are calibrated individually with the Shuffle Complex Evolution algorithm (SCE, Duan et al., 1992) to maximize the modified Kling Gupta Efficiency score (Kling et al., 2012, KGE'). The SCE is parameterized such that the maximum number of iterations and the number of complexes are respectively set to 100 and 25. CemaNeige free parameters are not calibrated; the default HOOPLA three-hour time step values are used instead (snowmelt factor θ G1 = 0.4 mm/ o C/3h and cold-content factor θ G2 = 0.93).
165 Figure 1 illustrates the performance values of the 20 models over the 31 catchments with the KGE'. Models perform generally well with an overall median KGE' close to 0.85. Performance varies more across catchments, which is expected considering that some are intrinsically harder to model. On this topic, the two catchments that exhibit a median performance below 0.6 are among the smallest ones (smaller than 800 km 2 ), which may indicate a possible limitation of HOOPLA regarding this aspect.

Simulation
Simulation is carried out in combination with one of the data assimilation techniques provided in HOOPLA, namely the Ensemble Kalman Filter. HOOPLA parameters are defined in order to issue 50 members and to update the hydrological models every eight time steps, i.e. once a day. Following Thiboult et al. (2016), the EnKF hyperparameters are defined so that the temperature uncertainty is assumed Gaussian with a standard deviation equal to 2 o C, the precipitation uncertainty 175 is approximated by a Gamma distribution, which standard deviation is equal to 50 % of the observed precipitation, and the streamflow uncertainty is assumed Gaussian as well, with a standard deviation equal to 10 % of the observed streamflow. These hyperparameters are shared by every model and catchment. Figure 2 confirms that HOOPLA provides a high level of performance over the catchments with an overall median KGE' greater than 0.90. Data assimilation successfully reduces the simulation bias as shown by the global improvement of perfor-180 mance over the calibration. This suggests that an adequate data assimilation schemes may help compensate for a calibration that sought minimizing the error over an entire chronicle. This gain varies according to catchments, but also to models to a lesser extent (Thiboult and Anctil, 2015).

Forecast
Forecast resorts to the same data assimilation setup than simulation. The 50-member meteorological ensemble forecast issued by the ECMWF is used as input for lead times up to 10 days. Therefore, each hydrological model ensemble is composed of 2500 members (as in Thiboult et al., 2016). Forecast performance decreases with increasing lead time (Figure 3, 4, and 5). The loss of accuracy mainly comes from the decrease of the meteorological predictive skills and the fading influence of data assimilation. Yet, for the 80 th lead time (10 days ahead), the median KGE' still lies around 0.55, indicating that the system remain an insightful predictor.

195
The multimodel approach thrives in short range forecasting and stands out from the other models as its median performance is greater for most lead times ( Figure 5). It may be occasionally outperformed by individual models, but the latter differ from one catchment to another.  10 https://doi.org/10.5194/hess-2020-6 Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License.
collection of tools and options to perform hydrological computations that fits one's needs. The toolbox includes among others 20 dissimilar lumped models, a snow accounting routine, three potential evapotranspiration formulas, two data assimilation schemes, two calibration algorithms, and parallel computing functionalities. Special attention was paid to create a framework that can host additional tools, making HOOPLA a suitable platform to develop and test tools designed by the users. Future 205 development will be driven by users' feedback and their contribution.
Code and data availability. The HOOPLA software, its manual, and a set of demo hydrometerological data are freely available under the BSD 2-Clause License at: https://github.com/AntoineThiboult/HOOPLA.

13
https://doi.org/10.5194/hess-2020-6 Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License. Table B1 provides an insight about computational time with a 3h time step. Tests were carried out over a 10-year period (29225 265 time steps) on a computer equipped with a processor Intel(R) Core(TM) i5-4210M clocked at 2.60 GHz and 8 Go of RAM.

Appendix B: Computation times
The experiment was repeated 100 times to minimize the influence of computer perturbations. As a reminder, HOOPLA offers a parallel computing option that may significantly decrease the time when multiple tasks are performed.  l'Information sur le Milieu Atmosphérique for their collaboration and for allowing the diffusion of the demo hydrometeorological observation dataset provided with HOOPLA. We are also thankful to the TIGGE initiative for granting us permission to distribute a subset of their meteorological forecast database as demo forcing in HOOPLA. Finally, we express our great appreciation to the IRSTEA Antony (HYCAR team), and in particular Charles Perrin, for sharing over the years, their extensive knowledge about the hydrological models.