Regional multi-objective calibration for distributed hydrological modelling: a decision tree based approach

Large scale modelling is becoming increasingly important in hydrology, particularly to characterize and quantify changes in the hydrological regime, whose drivers are typically large-scale phenomena, up to the global scale (e.g., climate change). This can be done with distributed models by estimating spatially consistent model parameters i.e. parameters having a functional relationship with catchment characteristics. In this study we adopt the newly developed PArameter Set Shuffling (PASS) approach, based on a machine learning decision tree algorithm, for the regional calibration of the TUWmodel over North-Western Italy. The method exploits observed patterns of locally calibrated parameters and catchment (climatic and geomorphological) descriptors, to derive functional relationships between the variables. The calibration procedure is performed by including snow cover information, as captured by MODIS datasets, in the model efficiency function. The results show that the PASS regionalization procedure allows to obtain very good regional model efficiencies, without significant loss of performance when moving from training to test catchments and from calibration to verification period, confirming the robustness of the methodology. We also highlight that using snow information in the calibration procedure is helpful to obtain spatially consistent model parameters for this study area. In the spirit of “obtaining good results for the right reasons”, this should be a preferred approach when performing the regional calibration of distributed hydrologic models over mountainous regions.

1 Introduction

One of the main challenges of hydrological modelling at the large scale is the estimation of spatially distributed model parameters, which are consistent with climatic and geomorphologic features of the territory, according to a process-based approach (Archfield et al., 2015; Clark et al., 2016; Gupta et al., 2014; Mizukami et al., 2017; Paniconi and Putti, 2015). The goal is to obtain parameter sets for each gridded element or hydrologic unit having the same functional relationship with climatic and geomorphological variables. Many studies dealt with parameters regionalization techniques, in particular within the context of runoff prediction in ungauged basins (e.g., Blöschl et al., 2013; Merz and Blöschl, 2004; Parajka et al., 2013; Seibert, 1999; Troch et al., 2003). An innovative approach was introduced by Merz et al. (2020), with the idea of using machine learning based algorithms to derive relationships between locally calibrated parameters and catchment descriptors, in order to predict parameters in a spatially distributed mode. The procedure, called PArameter Set Shuffling (PASS), provides good results for plain areas like Germany but, so far, it has not found application in other regions. This work presents the PASS approach in a more complex environment like the Alpine Region, by integrating snow information in the regionalization procedure.

2 Data and Methods

Data from the stream gauge network managed by the environmental protection agency of the Piemonte Region, in North-Western Italy, are used (Streamflow data, https://www.arpa.piemonte.it/rischinaturali/accesso-ai-dati/annali_meteoidrologici/annali-meteo-idro/banca-dati-idrologica.html, last access: 31 January 2023). The study area is characterized by a range of topographical elevations and a temperate continental climate, but conditions are not uniform in space. In the Alps, a tendency towards more humid conditions prevails, while the central area is characterized by lower annual precipitation and higher temperatures.

Figure 1 Stream gauges (red points), catchments boundaries (black lines) and pixels with meteorological inputs (blue cells) over Piemonte and Valle d'Aosta. The background represents mean elevation (m a.s.l.).

Precipitation [mm] and mean air temperature [°C] data are provided by a gridded dataset, obtained by interpolation of daily observations. These observations are taken from a dense network of meteorological stations and cover the period from 1957 up to the present day (Optimal Interpolation (OI) Dataset, http://www.arpa.piemonte.it/rischinaturali/tematismi/clima/confronti-storici/dati/dati.html, last access: 31 January 2023; Arpa Piemonte). The study domain consists in total of 305 pixels of size 0.125° × 0.125° (Fig. 1). Temperature data are used for the calculation of the potential evapotranspiration, by using the Blaney-Criddle equation based on the FAO procedure. A digital elevation model (DEM) at 90 m resolution (Digital elevation model, http://www.earthenv.org/, last access: 31 January 2023) and catchment boundaries for 197 catchments are also considered. Snow cover is mapped by combining MODIS products from Terra (MOD10A1) and Aqua (MYD10AQ) satellites, and in particular daily maps of the Normalized Difference Snow Index (NDSI) at 500 m spatial resolution. NDSI values larger than 40 are identified as snow. The combination is performed to reduce the effect of clouds; pixels which are classified as clouds or missing in Terra (NDSI > 100), are replaced by pixels from Aqua if these are classified as snow covered or snow free. The percentage of snow pixels in each grid cell of Fig. 1 is counted and the value is retained if the percentage of clouds is lower than 60 %. In this study, 79 catchment descriptors are used to infer regionalization rules, divided into 5 categories: climate (precipitation and temperature) statistics, calculated from the OI dataset; the 27 ETCCDI indices (http://etccdi.pacificclimate.org/list_27_indices.shtml, last access: 31 January 2023); morphological descriptors obtained from the digital elevation model (referred to elevation, slope and aspect); land-use classes derived from the Corine Land Cover dataset (2018) and a series of soil characteristics, obtained from the Harmonized World Soil Database.

Figure 2 (a) Non exceedance cumulative distribution function of model efficiencies (ME = 0.5 ⋅ KGE + 0.5 ⋅ SC) during the calibration period: “train_cal” shows MEs obtained by using locally calibrated lumped parameters for 76 training catchments, “test_cal” shows MEs obtained by using locally calibrated lumped parameters for 28 test catchments, “train_reg” shows MEs obtained by running the model with predicted distributed parameters for 76 training catchments, “test_reg” shows MEs obtained by running the model with predicted distributed parameters for 28 test catchments. (b) Non exceedance cumulative distribution function of model efficiencies (ME = 0.5 ⋅ KGE + 0.5 ⋅ SC) during the verification period. The continuous lines represent the mean model efficiency, while the dashed lines represent the 10th and 90th percentiles of model efficiencies for 30 parameter sets. (c) 30 normalized parameter sets found by local calibration for the Agogna river catchment at Momo, North-Western Italy. (d) 30 normalized parameter sets found by regional calibration for the Agogna river catchment at Momo, North-Western Italy.

The relationship between catchment descriptors and model parameters is established using decision trees (rpart R function). The PASS procedure is applied together with the TUWmodel (Appendix of Parajka et al., 2007), which follows the structure of the HBV model (Bergström, 1995). The algorithm used for local calibration is the differential evolution (DEoptim R package), which implements the global optimization of a real-valued function of a real-valued parameter vector (Mullen et al., 2011). We consider as model efficiency function the arithmetic average of two quantities:

• the Kling-Gupta Efficiency (KGE) for discharge:

  $$\mathrm{1}-\sqrt{{\left(r-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathrm{CV}}_{\mathrm{sim}}}{{\mathrm{CV}}_{\mathrm{obs}}}}-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\mu }}_{\mathrm{sim}}}{{\mathit{\mu }}_{\mathrm{obs}}}}-\mathrm{1}\right)}^{\mathrm{2}}}$$

  where r is the linear correlation coefficient between observed and simulated discharge, ${\mathrm{CV}}_{\mathrm{sim}}/{\mathrm{CV}}_{\mathrm{obs}}$ is the ratio between the coefficients of variation of simulated and observed discharges, ${\mathit{\mu }}_{\mathrm{sim}}/{\mathrm{\mu }}_{\mathrm{obs}}$ is the ratio between the means of simulated and observed discharges.

• a Snow Cover Efficiency, which is given by the following formula:

  $$\mathrm{SC}=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{n\cdot N}}\sum _{\mathrm{1}}^n\left(\sum _{\mathrm{1}}^N\left|{\mathrm{SC}}_{\mathrm{sim}}-{\mathrm{SC}}_{\mathrm{obs}}\right|\right)$$

  where N is the number of days, n is the number of cells, SC_sim is 1 if the simulated snow water equivalent (SWE) exceeds 1 mm and 0 otherwise, and SC_obs is 1 if the % of snow pixels within the cell is greater than 10 % and 0 otherwise.

The model efficiency is so defined:

$$\mathrm{ME}=\mathrm{0.5}\cdot \mathrm{KGE}+\mathrm{0.5}\cdot \mathrm{SC}$$

From the entire database, only stations having discharge data over the period 2000–2020 are considered for the analysis (104 stations). We split the timeseries into a calibration period (2010–2020) and a verification period (2000–2010). We perform the local calibration over the 104 catchments; from this set, only catchments with at least 5 parameter sets having local model efficiency (ME) > 0.75 and having a catchment area < 1000 km² are used for the training of PASS, following the assumptions of the method (Merz et al., 2020). These are the so called training catchments. For these sites, the best 30 parameter sets are used for running the regionalization. The remaining catchments, called test catchments, are instead used for validation.

3 Results and discussion

The results for the calibration period (2010–2020) and the verification period (2000–2010) are shown in Fig. 2a and b. Considering the calibration period, while the medians of the mean performances of training and test catchments in local calibration (“train cal” and “test cal”, respectively) are almost identical, applying the regionalization procedure, the median of MEs of regionally predicted distributed parameters sets for the training catchments (“train reg”) is 0.785, and the median regional efficiency for the test catchments (“test reg”) is 0.700.

Figure 3 Regional maps of distributed TUWmodel parameters obtained with PASS. The color indicates the mean value of the parameter among 30 regionalizations.

When moving to the verification period (2000–2010), the local model efficiencies have a substantial decrease. The median of regionally distributed MEs for training catchments decreases from 0.785 to 0.740, while the median of regionally distributed MEs for test catchments increases from 0.700 to 0.760. In general, we observe that the medians of regional MEs for training and test catchments are very close and the regional performance is not significantly degraded when moving from the calibration to the verification period. This confirms the robustness of the methodology in parameter estimation across spatial and temporal scales. In Fig. 2c and d, examples of normalized locally calibrated and regionally lumped calibrated parameter sets for the Agogna River at Momo are shown. Parameters found by regional calibration show a much smaller spread compared to locally calibrated parameters. As found in Merz et al. (2020), also for this Italian dataset regionalization reduces the effect of parameter equifinality.

Observing the maps of distributed model parameters, they appear to be highly consistent with the climatic and geomorphological features of the region; in particular snow parameters (SCF, DDF, Tr, Ts, Tm) seem to be in accordance with the three different geographical areas which characterize the region: the Alpine range on the West and North, the Apennines on the South and the Po River Valley (Fig. 3).

4 Conclusions

The results of this study indicate that the PASS procedure can be efficiently applied over the Alpine Region, with satisfactory results, also by taking advantage of snow cover information for the calibration and validation of the model. The model efficiencies for train catchments don't decrease significantly when moving from calibration to the verification period and actually for test catchments an increase is observed. This confirms the robustness of the procedure in the parameter estimation over different spatial and temporal domains. Overall, the results for this study area are in line with what has already been found in Merz et al. (2020). The spatial distribution of the parameters is consistent with the climatic and geomorphological characteristics of the region and this confirms that parameters are controlled by these variables. Following a process-based approach, the use of snow information can be crucial for a good model representation of the processes. This will allow to get “right results for the right reasons”, which is of particular importance in hydrology, especially for large-scale distributed modelling studies which require ungauged streamflow prediction.