A geostatistical data-assimilation technique for enhancing macro-scale rainfall–runoff simulations

Our study develops and tests a geostatistical technique for locally enhancing macro-scale rainfall–runoff simulations on the basis of observed streamflow data that were not used in calibration. We consider Tyrol (Austria and Italy) and two different types of daily streamflow data: macro-scale rainfall–runoff simulations at 11 prediction nodes and observations at 46 gauged catchments. The technique consists of three main steps: (1) period-of-record flow–duration curves (FDCs) are geostatistically predicted at target ungauged basins, for which macro-scale model runs are available; (2) residuals between geostatistically predicted FDCs and FDCs constructed from simulated streamflow series are computed; (3) the relationship between duration and residuals is used for enhancing simulated time series at target basins. We apply the technique in cross-validation to 11 gauged catchments, for which simulated and observed streamflow series are available over the period 1980–2010. Our results show that (1) the procedure can significantly enhance macro-scale simulations (regional LNSE increases from nearly zero to ≈0.7) and (2) improvements are significant for low gauging network densities (i.e. 1 gauge per 2000 km²).

1 Introduction

The steady increase in computational capabilities together with the expanding accessibility of regional and global datasets (e.g. soil properties, land-cover, morphology, climate characteristics, satellite-based gridded precipitation) trigger the development of regional- to continental-scale and global-scale hydrological models (Archfield et al., 2015), hereafter referred to as macro-scale models.

During the last decade, several of these macro-scale models have become operational and thus continuously provide data automatically for decision-making. For instance, the distributed rainfall–runoff-routing model LISFLOOD (De Roo et al., 2000) provides daily forecast for operational warning services through the systems of EFAS (Pappenberger et al., 2013) and GLOFAS (Alfieri et al., 2013); the LARSIM models (Haag and Luce, 2008) are used operationally for simulating streamflow at large areas in southern Germany, Luxembourg, Austria, Switzerland, and the eastern part of France; the WATFLOOD, developed at the University of Waterloo, is used operationally in Canada (Kouwen et al., 1993); and the S-HYPE model (Strömqvist et al., 2012) is running operationally for flood or drought forecasting and water quality assessments for the Swedish landmass, providing high-resolution information to authorities and citizens (Hjerdt et al., 2011).

Other macro-scale models are used for off-line water assessments and research purposes. For instance, the global WaterGAP Global Hydrological model (Alcamo et al., 2003) assists in water accounting; the SAFRAN-ISBA-MODCOU model (Habets et al., 2008) has been applied over the entire French territory to combine a meteorological analysis system, a land surface model, and a hydrogeological model; the PGB-IPH model (Pontes et al., 2017) has been applied to many South American basins; and the SWIM model (Krysanova et al., 1998) couples water balance simulations with water quality for small to mid-size watersheds, i.e. regional meso-scale.

The macro-scale hydrological models are getting more and more popular due to three main reasons: (1) they can provide users with a large-scale representation of hydrological behaviour, which is fundamental information for effectively addressing several water resources planning and management problems (e.g. surface water availability assessment, instream water quality studies, ecohydrological studies); (2) they can be used to compute a variety of hydrological signatures everywhere along the stream network at the resolution of the model; (3) model outputs in some cases are open-access and freely distributed, so that their regional runs represent a wealth of information for addressing the problem of hydrological predictions in data-scarce regions of the world (Beck et al., 2016; Donnelly et al., 2016; Pechlivanidis and Arheimer, 2015). Accurate regional hydrological simulations undoubtedly foster and support the implementation of improved large-scale and trans-boundary policies for water resources system management and flood-risk mitigation or climate change adaptation (Arheimer et al., 2017; Falter et al., 2016; Sampson et al., 2015; de Paiva et al., 2013).

However, improved accuracy in terms of average regional performance does not necessarily imply homogeneous improvements in local performance. In fact, due to the difficulties to perform local calibrations and validations of macro-scale models over the entire modelled regions, local performance can be rather diverse (Donnelly et al., 2016; de Paiva et al., 2013). Factors controlling the heterogeneity of local performance may be various, for instance the quality of macro-scale input data, local water management, representativeness of model structure chosen, the influence of local geophysical and micro-climatic factors.

There is a recognized and noteworthy value in readily available and easy accessible simulated daily streamflow series for scarcely gauged, or ungauged, areas of the world to enhance awareness and decision-making (Arheimer et al., 2011; Hjerdt et al., 2011). Nevertheless, the harmonization and enhancement of local performances of macro-scale models is still a scientific challenge that is worth addressing in operational hydrology, and which raises different research questions, such as the following:

• How could we deal with locally biased simulations?

• Can we assimilate additional data to improve model performance without re-calibration?

• Is there a minimum gauging network density that makes the post-modelling data assimilation viable and effective?

Recent literature shows the significant potential of kriging-based techniques for performing regional prediction of streamflow indices in ungauged locations (Castiglioni et al., 2011; Pugliese et al., 2014; Skøien et al., 2006). Among such techniques, topological kriging, or top-kriging (Skøien et al., 2006), has shown high prediction accuracy and excellent adaptability to a variety of water-related applications, such as prediction of low-flow indices (Castiglioni et al., 2009), interpolation of river temperatures (Laaha et al., 2013), estimation of flood quantiles (Archfield et al., 2013), regionalization of flow–duration curves (Castellarin, 2014; Castellarin et al., 2018; Pugliese et al., 2014, 2016), estimation of daily runoff in ungauged basins (Parajka et al., 2015) and reconstruction of historical daily streamflow series (Farmer, 2016).

Our study aims to develop and test a geostatistical data-assimilation procedure for better agreement between locally observed streamflow and model results from macro-scale rainfall–runoff models. The procedure employs top-kriging for geostatistically interpolating empirical period-of-record flow–duration curves (FDCs) along the stream network available at gauged basins. Interpolated FDCs are assimilated into simulated daily streamflow series at ungauged stream-network nodes, enhancing the local accuracy of simulated daily streamflow series. We test our method by improving E-HYPE model simulations (Donnelly et al., 2016; Hundecha et al., 2016), which provides approximately 30 years of simulated daily streamflows freely and openly accessible for 35 408 prediction nodes in Europe (mean catchment size of 215 km², see also http://hypeweb.smhi.se/europehype/time-series/, last access: 30 July 2018). We address the Tyrolean region, as this area gives particularly poor simulations in the E-HYPE version 3 and thus would benefit from statistical enhancement of results. For the geostatistical interpolation, we use a group of 46 gauged catchments obtained from Austrian and Italian water services, and not used when setting up E-HYPE. With the observed streamflows we construct and interpolate empirical FDCs. Then, we assess the value and potential of assimilating this streamflow information into E-HYPE simulated series. In particular, (1) we cross-validate the proposed data-assimilation procedure for 11 E-HYPE prediction nodes located nearby an existing stream gauge, and (2) we assess the enhancement of simulated series resulting from the geostatistical data assimilation under different hypotheses on the spatial density of the stream-gauging network.

The paper is structured as follows: first, methods and procedures are presented in a general way, then we illustrate the case study and application. In particular, Sect. 2 presents the proposed procedure, while Sect. 3 details the system of cross-validations and sensitivity analyses we adopted for assessing the procedure. Section 4 illustrates the study area and E-HYPE simulation data. The last three sections report results, discussion and conclusions, respectively.

2 A new geostatistical streamflow-data assimilation method

2.1 Geostatistical interpolation of empirical flow–duration curves (TNDTK)

Top-kriging is a powerful geostatistical procedure developed by Skøien et al. (2006) for the prediction of hydrological variables. Like all kriging approaches, top-kriging produces predictions of hydrological variables at ungauged sites with a linear combination of the empirical information collected at neighbouring gauging stations. Through this method, the unknown value of the streamflow index of interest at prediction location x₀, Z(x₀) can be estimated as a weighted average of the regionalized variable, measured within the neighbourhood:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }Z\left(x_{\mathrm{0}}\right)=\sum _{j=\mathrm{1}}^n{\mathit{\lambda }}_{j}Z\left(x_j\right),\end{array}$$

where λ_j is the kriging weight for the empirical value Z(x_j) at location x_j, and n is the number of neighbouring stations used for interpolation. Kriging weights λ_j can be found by solving the typical ordinary kriging linear system (Eq. 2a) with the constraint of unbiased estimation (Eq. 2b):

$$\begin{array}{}\text{(2a)}& {\displaystyle }& {\displaystyle }\sum _{j=\mathrm{1}}^n{\mathit{\gamma }}_{i,j}{\mathit{\lambda }}_j+\mathit{\theta }={\mathit{\gamma }}_{\mathrm{0},i}i=\mathrm{1},\mathrm{\dots },n,\text{(2b)}& {\displaystyle }& {\displaystyle }\sum _{j=\mathrm{1}}^n{\mathit{\lambda }}_j=\mathrm{1},\end{array}$$

where θ is the Lagrange parameter and γ_i,j is the semi-variance between catchment i and j (Isaaks and Srivastava, 1990). The semi-variance, or variogram, represents the spatial variability of the regionalized variable Z. Unique from any other method of kriging, top-kriging considers the variable defined over a non-zero support S, the catchment drainage area (Cressie, 1993; Skøien et al., 2006). The kriging system of Eqs. (2a) and (2b) remains the same, but the semi-variances between the measurements need to be obtained by regularization, i.e. smoothing the point variogram over the support area.

Figure 1 A sketch of the total negative deviation (TND).

The point variogram can then be back-calculated by fitting aggregated variogram values to the sample variogram (Skøien et al., 2006). Pugliese et al. (2014) proposed a method for using top-kriging to predict FDCs at ungauged locations that they termed total negative deviation top-kriging (TNDTK). The authors reduce the dimensionality of the problem by seeking a unique index of site-specific FDCs. Unlike other regional approaches (Castellarin et al., 2013), the kriging-based method interpolates the entire curve, therefore ensuring its monotonicity (Castellarin, 2014; Pugliese et al., 2014). This is accomplished by first standardizing the empirical FDCs at site x, Ψ(x, d), for some reference value, Q^*(x), to yield a dimensionless FDC:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\mathit{\psi }(x,d)={\displaystyle \frac{\mathrm{\Psi }(x,d)}{Q^{*}\left(x\right)}},\end{array}$$

where d denotes a specific duration. Pugliese et al. (2014) identified an overall point index that effectively summarizes the entire curve. This index, which the authors termed total negative deviation (TND), is derived by integrating the area between the lower limb of the FDC and the reference streamflow value Q^* (see Fig. 1). Empirical TND values are computed as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathrm{TND}\left(x\right)=\sum _{i=\mathrm{1}}^m\left|q_i\left(x\right)-\mathrm{1}\right|{\mathit{\delta }}_i,\end{array}$$

where $q_i=\frac{Q_i}{Q^{*}}$ represents the ith empirical dimensionless quantile standardized for the selected reference value Q^*, δ_i is half of the frequency interval between the (i+1)th and (i−1)th quantile and the summation involves only the m standardized quantiles lower than 1. The equality between a given streamflow value and the reference value Q^* is represented by a horizontal dashed line in Fig. 1, i.e. the threshold given by the equation $\frac{Q}{Q^{*}}=\mathrm{1}$. The range of the summation, m, in Eq. (4) is a function of the maximum duration d_max, which is itself a function of that sample with minimum length across gauged sites in the study region. Having calculated empirical TNDs, Pugliese et al. (2014) propose using the TNDs as a regionalized variable to develop site-specific weighting schemes. The same weights, derived through the solution of the linear kriging system (Eqs. 2a and 2b), are used for a batch prediction of the continuous, dimensionless FDC for the ungauged site x₀:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\widehat{\mathit{\psi }}\left(x_{\mathrm{0}},d\right)=\sum _{j=\mathrm{1}}^n{\mathit{\lambda }}_j\mathit{\psi }\left(x_j,d\right)\forall d\in (\mathrm{0},\mathrm{1}),\end{array}$$

where λ_j, with j=1, …, n, is the weights resulting from the kriging interpolation of TNDs, ψ(x_j, d) is the dimensionless empirical FDC at the donor site x_j, and $\widehat{\mathit{\psi }}(x_{\mathrm{0}}$, d) is the predicted dimensionless FDC. It is worth highlighting that the computation of the linear kriging system (Eqs. 2a and 2b) depends on n, the number of neighbouring sites on which to base the spatial interpolation, a fact that will be explored below.

If a reliable model for predicting Q^* at the ungauged site x₀ can be developed, the prediction of the dimensional FDC, $\widehat{\mathrm{\Psi }}(x_{\mathrm{0}}$, d), is obtained as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\widehat{\mathrm{\Psi }}\left(x_{\mathrm{0}},d\right)={\widehat{Q}}^{*}\left(x_{\mathrm{0}}\right)\widehat{\mathit{\psi }}\left(x_{\mathrm{0}},d\right)\forall d\in (\mathrm{0},\mathrm{1}),\end{array}$$

where ${\widehat{Q}}^{*}\left(x_{\mathrm{0}}\right)$ is the prediction of Q^* at the ungauged site x₀ and $\widehat{\mathit{\psi }}(x_{\mathrm{0}}$, d) has the same meaning as in Eq. (5).

The main objective of our study is to improve rainfall–runoff model simulations, whereas an assessment of the reliability of the regional metric TND is out of the scope of this paper. Nevertheless, further details on TNDTK (i.e. cross-validations in different geomorphological and climatic regions, sensitivity analyses, comprehensive assessments and comparisons to state-of-the-art models for predicting FDCs in ungauged basins) are reported in recent studies to which an interested reader is referred (Kim et al., 2017; Pugliese et al., 2014, 2016).

2.2 Algorithm for assimilation of local streamflow data

Following the approach proposed by Smakhtin and Masse (2000), we present a novel procedure for predicting the model residuals that may be associated with macro-scale rainfall–runoff model simulations (e.g. LISFLOOD, HYPE, PGB-IPH; see Introduction). This method relies on a regional prediction of the long-term FDC in the same site where these simulations are available.

For instance, let Ψ(x₀, d) be the “true” unknown FDC for a given catchment x₀ and ${\widehat{\mathrm{\Psi }}}_{\mathrm{SIM}}(x_{\mathrm{0}}$, d) be its prediction constructed on the basis of the daily streamflows simulated through the macro-scale model. We can assume that a general relationship between the two curves exists and reads as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\mathrm{\Psi }\left(x_{\mathrm{0}},d\right)={\widehat{\mathrm{\Psi }}}_{\mathrm{SIM}}\left(x_{\mathrm{0}},d\right)+\mathit{\epsilon }\left(x_{\mathrm{0}},d\right)\forall d\in (\mathrm{0},\mathrm{1}),\end{array}$$

where ε(x₀, d) are the model residuals defined over the duration domain d, which we may term the residual–duration curve (εDC). Evidently, although the “true” residual–duration curve is unknown at ungauged basins, one can nevertheless estimate such a curve on the basis of geostatistically interpolated flow–duration curves ${\widehat{\mathrm{\Psi }}}_{\mathrm{TNDTK}}(x_{\mathrm{0}}$, d) introduced in Sect. 2.1,

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }\widehat{\mathit{\epsilon }}\left(x_{\mathrm{0}},d\right)={\widehat{\mathrm{\Psi }}}_{\mathrm{TNDTK}}\left(x_{\mathrm{0}},d\right)-{\widehat{\mathrm{\Psi }}}_{\mathrm{SIM}}\left(x_{\mathrm{0}},d\right)\forall d\in (\mathrm{0},\mathrm{1}).\end{array}$$

The estimated residual–duration curve obtained from the regional prediction of the long-term flow–duration curve can then be used for assimilating local streamflow information into the simulated daily streamflow series. The procedure is sketched in Fig. 2: (1) given a simulated streamflow series (red line in the top-right), select a specific day t and the corresponding discharge Q(t); (2) retrieve the duration d associated with Q(t) from the flow–duration curve constructed from simulated data (red line in the top-left quadrant); (3) read the estimated residual $\widehat{\mathit{\epsilon }}\left(t\right)$ off of the predicted residual–duration curve (blue line in the bottom-left quadrant); and (4) assimilate the residual into the simulated series as $Q\left(t\right)+\widehat{\mathit{\epsilon }}\left(t\right)$. The iteration of the algorithm through all time steps leads to an enhanced simulated series (blue line in the top-right quadrant).

This new assimilation procedure shares some analogies with a technique called “quantile mapping”, used in the context of bias corrections for global climate model predictions (Komma et al., 2007). The procedure we propose in this context, though, is a rather general tool that can be applied to, for example, any macro-scale rainfall–runoff model for locally enhancing long simulated streamflow series without the need to re-run computationally intensive simulations, provided that the model itself is behavioural and validated on the basis of streamflow data that were not available for model calibration. The performance of the assimilation procedure depends on a variety of drivers, e.g. the quality of the simulated streamflows, which can be severely impacted by the local quality of input data even for a behavioural and well-calibrated model, the quality of streamflow data and the density of the stream-gauging network (see Sect. 3.1.3), the accuracy of the chosen regional model for predicting FDCs (Castellarin, 2014; Castellarin et al., 2013). Regarding the latter element, indeed, the proposed procedure deeply relies on accurate regional FDC predictions by means of an unbiased regional model, which necessarily must be validated beforehand for the area of interest. Otherwise, detriments are likely to be expected.

Figure 2 Illustration of the proposed data-assimilation procedure for a given simulated time series. Top-right panel: real streamflow series (unknown, since the basin is ungauged, black dashed line); macro-scale model simulation (red solid line); geostatistically enhanced streamflow series (blue solid line). Top-left panel: FDC predictions obtained from simulated streamflows (red solid line) and via geostatistical interpolation (black solid line), with real (unknown) FDC (black dashed line). Bottom-left panel: estimated residual–duration curve (blue solid line) computed as the difference between the two predicted FDCs in the top-left panel. Bottom-right panel: time series of residuals (blue dashed line).

3 Testing the proposed algorithm: cross-validation procedures and sensitivity analysis

3.1 Structure of the analysis

The rationale of the analyses implemented in this study drives a sequence of operations, which can be summarized as follows:

1. Streamflow simulations from a given model, e.g. a rainfall–runoff macro-scale model, are available in a well-defined study area.

2. We suppose that recorded streamflow data are made available for a reasonable number of stream gauge stations within the study area, and we apply a suitable model for the regionalization of FDCs.

3. We validate the regional model with respect to available streamflow observations. In this case we adopted a leave-one-out cross-validation (Kroll and Song, 2013; Salinas et al., 2013; Srinivas et al., 2008; Wan Jaafar et al., 2011), even though different validation schemes might be preferred in other regions (Castellarin et al., 2018; Pugliese et al., 2016).

4. We validate the assimilation procedure by sequentially (a) neglecting all the streamflow information at a given nodes of the river network, (b) predicting FDCs, and (c) applying the assimilation method illustrated in Sect. 2.2.

5. We evaluate the sensitivity of the assimilation procedure to the stream gauge density of the study area.

The methodologies adopted for addressing points 3–5 are illustrated in Sects. 3.1.1, 3.1.2, and 3.1.3, in this order; the accuracy of predictions (e.g. regional FDCs, assimilated and simulated streamflow series) relative to their empirical or observed counterparts is quantified through performance indices described in Sect. 3.2.

3.1.1 Cross-validation of the FDC geostatistical interpolator (point 3 above)

We proposed to assess the accuracy of the geostatistical predictor of FDCs (i.e. TNDTK, see Sect. 2.1) in cross-validation with respect to available streamflow observations from a sufficient number of gauged catchments. We chose the mean annual flow (MAF), computed as the average flow of recorded historical streamflow series, as the reference value Q^* (see details in Sect. 2.1).

TNDTK operates by first applying top-kriging to empirical TND values (see Sect. 2.1), which we performed by calculating a binned sample variogram first, and then by modelling binned empirical data with a five-parameter “modified” exponential theoretical variogram (Skøien et al., 2006). The fitted theoretical point variogram and its five parameters were obtained through the weighted least squares (WLS) regression method from Cressie (1993) by simultaneously fitting all regularized binned variograms that were computed for various area classes (Skøien, 2018). Recent applications of TNDTK indicate n=6 as an optimal number of neighbouring donor stations, and thus we chose the same value for this case study as well (Pugliese et al., 2014, 2016). Then, TNDTK uses the kriging weights obtained for predicting TND values for interpolating the dimensionless FDCs at the location of interest as the weighted average of dimensionless empirical FDCs constructed from the n=6 neighbouring gauged sites (see Eq. 5, in which λ_j, with j=1, …, 6 and n=6, is the kriging weights). While the computation of TNDs does not require any specific resampling scheme of the FDCs, the prediction of the curves in ungauged locations is carried out using a fixed number of quantiles that should be selected to thoroughly represent the variability from high to low flows. Thus, observed FDCs are resampled to 20 equally spaced points, in the normal space, leading to the widest range of durations compatible with the shorter observed streamflow series in the dataset. We adopted a LOOCV procedure (Pugliese et al., 2014, 2016) to test the accuracy and uncertainty associated with FDC predictions. This simulates the ungauged conditions at each and every gauged site in the study area by (1) removing it in turn from the dataset and (2) referring to the n=6 neighbouring gauges for predicting its dimensionless FDC. Given that the geostatistical assimilation procedure uses dimensional FDCs, we also tested the suitability of standard top-kriging for predicting MAF at ungauged locations in the study area, still through a LOOCV procedure (Blöschl et al., 2013). For MAF interpolation, we adopted the same settings used for predicting TND values (i.e. a five-parameter “modified” exponential theoretical variogram and n=6 neighbouring sites). We then used cross-validated dimensionless FDCs and MAF predictions at each and every gauging station in the study area to obtain cross-validated predictions of dimensional FDCs for each measuring node through Eq. (6).

3.1.2 Cross-validation of the geostatistical assimilation procedure (point 4 above)

We applied the proposed assimilation procedure as outlined in Sect. 2.2, by, firstly, assessing the efficiency of the procedure through a leave-one-out cross-validation. We predicted the FDC associated with each simulation node by using TNDTK and by, also, neglecting the hydrological information coming from the closest (i.e. immediately upstream or downstream) gauged catchment, therefore assuming that no streamflow information is available near the simulation node. The workflow of the validation algorithm is as follows:

• a.

  We select one pair among n_p possible pairs, let us term it pair i–j, where i stands for simulation node and j stands for the corresponding stream gauge.

• b.

  We drop the daily streamflow series observed at stream gauge j from the set of observed series.

• c.

  We interpolate FDC at simulation node i through TNDTK (as illustrated in Sect. 2.1) using the remaining n_g−1 gauged sites, where n_g is total number of stream gauges in the study region.

• d.

  We apply the assimilation procedure outlined in Sect. 2.2 and depicted in Fig. 2 to the streamflow series simulated for the simulation node i.

• e.

  We compare the original simulated daily streamflow series and the geostatistically enhanced one at prediction node i with the daily streamflow series observed at stream gauge j times the corresponding area ratio A_i∕A_j (i.e. A_i is the drainage area of simulation node i, A_j is the drainage area of stream gauge j; see also Sect. 4.1).

• f.

  We repeat all previous steps for each one of the remaining n_p pairs.

For the sake of consistency, we anticipate here that we will refer to the procedure presented above as GAE-HYPE (i.e. geostatistically assimilated E-HYPE) in the remainder of the paper. The acronym clearly recalls the rainfall–runoff model used in this study (see Sect. 4.2); however, the procedure disregards a specific rainfall–runoff model.

Finally, it is worth highlighting here that FDCs obtained from either the geostatistical model or the rainfall–runoff simulation model are resampled to 20 equally spaced points across the normally transformed duration intervals (Pugliese et al., 2014). Thus, as a result, the produced εDCs reflect the same sampling scheme of the curves. Nevertheless, the procedure does not foresee any restriction to the resolution of the resampled curve, allowing for a finer resampling scheme in other analyses.

3.1.3 Stream-gauging network density and effectiveness of geostatistical data assimilation (point 5 above)

Since the proposed geostatistical data-assimilation procedure (GAE-HYPE) relies upon the local availability of stream gauges records, understanding to what extent the performances of the assimilation method are driven by gauging network density is a fundamental of paramount importance. Therefore, we performed a sensitivity analysis and assessed the degree of enhancement of simulated daily streamflow sequences associated with different scenarios of streamflow data availability, repeating for each scenario the procedure described in Sect. 3.1.2. Thus, we randomly discarded some of the gauges available over the study area and varied the total number of available stations continuously from the lowest to the highest gauge density. At each density scenario, we performed exactly the same kriging settings and same LOOCV illustrated in detail in Sect. 3.1.1 and 3.1.2, respectively.

Figure 3 Study area: Tyrol. Catchment boundaries for 11 E-HYPE prediction nodes (red) and 46 stream gauges (black).

Table 1 Study catchments: streamflow properties standardized by drainage area (m³ s⁻¹ km⁻²) for either gauged catchments or E-HYPE catchments, mean annual flow (q_MAF), and 50 % and 95 % streamflow quantiles (q₅₀ and q₉₅, respectively).

3.2 Performance indices

We assessed the performances for predicting regional FDCs by means of Nash–Sutcliffe efficiency (Nash and Sutcliffe, 1970) computed for log-transformed streamflows (LNSEs). These indices are computed as follows:

$$\begin{array}{ll}{\displaystyle }{\mathrm{LNSE}}_{\mathrm{FDC},j}=& {\displaystyle }\mathrm{1}-{\displaystyle \frac{\sum _{k=\mathrm{1}}^{n_{\mathrm{d}}}{\left(\ln \mathrm{\Psi }\left(x_j,d_k\right)-\ln \widehat{\mathrm{\Psi }}\left(x_j,d_k\right)\right)}^{\mathrm{2}}}{\sum _{k=\mathrm{1}}^{n_{\mathrm{d}}}{\left(\ln \mathrm{\Psi }\left(x_j,d_k\right)-{\mathit{\mu }}_j\right)}^{\mathrm{2}}}},\\ \text{(9)}& {\displaystyle }& {\displaystyle }j=\mathrm{1},\mathrm{\dots },n_{\mathrm{g}},\end{array}$$

where Ψ(x_j, d_k) and $\widehat{\mathrm{\Psi }}(x_j$, d_k) are the empirical and the predicted kth streamflow quantiles at site x_j, respectively, μ_j is the mean of empirical log-transformed streamflow quantiles at site x_j, n_d is the number of discretization points throughout duration range, and n_g is the number of stream gauges.

Another useful metric of performance for the assessment of FDC predictions is the overall absolute curve error (Ganora et al., 2009), which reads as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_{\mathrm{FDC},j}=\sum _{k=\mathrm{1}}^{n_{\mathrm{d}}}\left|\mathrm{\Psi }\left(x_j,d_k\right)-\widehat{\mathrm{\Psi }}\left(x_j,d_k\right)\right|j=\mathrm{1},\mathrm{\dots },n_{\mathrm{g}},\end{array}$$

where Ψ(x_j, d_k) and $\widehat{\mathrm{\Psi }}(x_j$, d_k) have the same meaning as in Eq. (9).

Similarly, concerning streamflow time series, the assessment of modelled streamflows is carried out with LNSE, but, in this case, it reads as follows:

$$\begin{array}{ll}{\displaystyle }{\mathrm{LNSE}}_{\mathrm{mod},j}=& {\displaystyle }\mathrm{1}-{\displaystyle \frac{\sum _{t=\mathrm{1}}^{t_{\mathrm{s}}}\left(\ln Q_{\mathrm{emp},j}\left(t\right)-\ln Q_{\mathrm{mod},j}\left(t\right)\right)}{\sum _{t=\mathrm{1}}^{t_{\mathrm{s}}}{\left(\ln Q_{\mathrm{emp},j}\left(t\right)-{\mathit{\omega }}_j\right)}^{\mathrm{2}}}},\\ \text{(11)}& {\displaystyle }& {\displaystyle }j=\mathrm{1},\mathrm{\dots },n_{\mathrm{p}},\end{array}$$

where Q_emp,j(t) and Q_mod,j(t) are the empirical and predicted streamflow at site x_j and time t, respectively, ω_j is the mean of empirical log-transformed streamflow at site x_j, t_s and n_p is the number of selected simulation nodes.

Furthermore, we assessed the efficiency of the data-assimilation procedure through the following metric:

$$\begin{array}{ll}{\displaystyle }{\mathrm{LNSE}}_{\mathrm{ratio},j}=& {\displaystyle \frac{{\mathrm{LNSE}}_{\mathrm{mod}\mathrm{2},j}-{\mathrm{LNSE}}_{\mathrm{mod}\mathrm{1},j}}{\mathrm{1}-{\mathrm{LNSE}}_{\mathrm{mod}\mathrm{1},j}}},\\ \text{(12)}& {\displaystyle }& {\displaystyle }j=\mathrm{1},\mathrm{\dots },n_{\mathrm{p}}.\end{array}$$

LNSE_ratio quantifies the degree of enhancement of “model 2” relative to “model 1” standardized by the maximum possible improvement (i.e. 1−LNSE_mod1). An LNSE_ratio close to zero means no significant enhancement (detriment of original sequences in case of negative values), whereas an LNSE_ratio close to 1 indicates that no further enhancement is possible. Such an index is derived from the reciprocal root-mean-squared error ratio between the two models.

Finally, in order to verify whether or not the proposed assimilation procedure GAE-HYPE outperforms rainfall–runoff simulations, for different gauge density scenarios (see Sect. 3.1.3), we used the Wilcoxon signed-rank test with the null hypothesis that simulation model LNSEs are greater than GAE-HYPE ones at 5 % significance level (Hollander and Wolfe, 1999).

4 Data

4.1 Study region

We focus on a large alpine region located in Tyrol (Italy, Austria and, for small portion only, Switzerland). Our analyses consider two types of data, observed daily streamflows and E-HYPE simulated daily streamflows (see Sect. 4.2), representing different sets of catchments (Fig. 3). In this study, E-HYPE represents the rainfall–runoff model selected to evaluate the procedure presented in Sect. 2.2, which has shown significantly poor results in this region. Indeed, this alpine area is particularly suitable for hydro-power generation, and therefore the presence of dams along the stream network could likely alter the streamflow regime downstream, producing a significant alteration of the natural flow conditions. E-HYPE only simulates the dams present in the global database of GranD (Lehner et al., 2011), which might not be representative for hydropower production at the local scale (Arheimer et al., 2017). Thus, we removed from the initial group of gauged catchments all basins for which the streamflow regime is highly or significantly altered by upstream dams. Table 1 reports the main characteristics of streamflow regimes for 46 gauged basins and 11 selected E-HYPE prediction nodes. Among all E-HYPE prediction nodes available in Tyrol we selected only those whose catchments were the closest to gauged ones, i.e. difference in terms of drainage areas < 14 % and distance between catchment centroids < 15 km. These criteria resulted in the selection of 11 E-HYPE prediction nodes that are evenly distributed in the study region (see red lines in Fig. 3). We addressed the limited differences existing between drainage areas of E-HYPE and gauged basins by adopting the drainage-area ratio technique (Farmer and Vogel, 2013), that is by rescaling daily streamflows according to drainage areas of the corresponding catchment. Such a method assumes the same unit daily streamflow for any pair of hydrologically similar catchments i and j, which reads as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle \frac{Q_i\left(t\right)}{A_i}}={\displaystyle \frac{Q_j\left(t\right)}{A_j}},\end{array}$$

where Q_i(t) represents the daily streamflow at day t for catchment i, and Q_j(t) is the daily streamflow at day t for catchment j. In our application, i and j could correspond to any given pair (stream gauge, E-HYPE prediction node), and A_i and A_j the corresponding drainage areas. Finally, it is worth pointing out that neither the observed nor the E-HYPE series present zero values in their recording (or simulation) periods for each of the 11 selected nodes.

Figure 4 Schematic concept of the HYPE model (all equations are available at http://hypecode.smhi.se/, last access: 30 July 2018).

4.2 Pan-European rainfall–runoff simulation: E-HYPE

The HYdrological Predictions for the Environment (HYPE) model is a hydrological model for small-scale and large-scale assessments of water resources and water quality, developed at the Swedish Meteorological and Hydrological Institute (SMHI) during 2005–2007 (Lindström et al., 2010). The European application, E-HYPE, has been proved to be a powerful tool for water resources managers and practitioners, addressing nutrient concentration in river flow as well as water forecasts on short or seasonal timescale. It is also widely used to estimate snow storage and accumulated TWh (terawatt hour) of water inflow to hydropower dams and in climate change impact analysis (Donnelly et al., 2017). The website Hypeweb (http://hypeweb.smhi.se, last access: 30 July 2018) provides visualization and free downloading of 30 years of continuous and consistent daily streamflow simulations across the European river network at rather fine scale (i.e. the average size of elementary catchments is equal to 215 km²) as well as forecasts and climate change impact analysis.

The HYPE model is open-access and can be downloaded with documentation and model set-up guidelines from the model website (http://hypecode.smhi.se/, last access: 30 July 2018). It simulates water flow and substances on their way from precipitation through soil, river and lakes to the river outlet (Lindström et al., 2010). River basins are divided into sub-basins, which in turn are divided into classes (the finest calculation units) depending on land use, soil type and elevation (Fig. 4).

The soil is modelled as several layers, which may have different thicknesses for each class. In E-HYPE, each sub-basin can have up to some 40 soil and land-use classes, which are lumped within the sub-basins, while the watercourses are routed through the river network. The model parameters can be associated with land use (e.g. evapotranspiration) or soil type (e.g. water content in soil), or be common for the whole catchment or a region with geophysical similarities (Hundecha et al., 2016). This way of coupling the parameters with geographic information makes the model better suited for simulations in ungauged catchments.

5 Results

The application of the geostatistical method TNDTK through an LOOCV procedure reveals an agreement between empirical values and predictions as shown in Figs. 5 and 6. Specifically, Fig. 5 reports empirical (x axis) against geostatistically predicted (y axis) MAF values as well as LNSEs obtained in cross-validation (i.e. 0.96). Figure 6a shows a scatter diagram between observed (x axis) and predicted streamflows (y axis) from FDCs. This agreement is also confirmed by the distributions of on-site LNSE values (see box plots in Fig. 6b); the median LNSE is equal to 0.97, while mean LNSE is ca. 0.90. The performance obtained in cross-validation legitimizes the use of TNDTK for predicting FDCs in the study area at the 11 E-HYPE prediction nodes of interest, for which TNDTK delivers high prediction capability, with LNSE values above 0.97 (see the spatial distribution of efficiency values in Fig. 10a).

Figure 5 Top-kriging predictions of mean annual flow (MAF) in cross-validation mode.

Table 2 Nash–Sutcliffe Efficiencies computed on log-transformed daily streamflows for E-HYPE and GAE-HYPE: median values for the 11 prediction nodes considered in the study; smallest enhancement (IDs 3675 – 9001070), largest enhancement (IDs 201236 – 9608296).

Figure 6 (a) Scatter diagrams of empirical (x axis) vs. predicted (y axis) streamflows. (b) Box-plot representation of on-site LNSE values, summarizing the first, second (median) and third quartiles along with whiskers extending to the most extreme non-outlying data point (outliers are highlighted as circles and lay at more than 1.5 times the interquartile-range from the nearest quartile); the average on-site LNSE value is reported in (a) and illustrated as a dashed line in (b).

Figure 7 11 E-HYPE prediction nodes: (a) comparison between TNDTK (x axis) and E-HYPE (y axis) in terms of distances δs between empirical and predicted FDCs; the 1:1 line represents equivalent performance for TNDTK and E-HYPE; (b) standardized residual–duration curves computed as illustrated in Eq. (8) (TNDTK streamflow quantiles are used for standardization).

Figure 7 reports the results obtained by applying the aforementioned cross-validation algorithm. Circles in Fig. 7a represent the cumulative absolute error δ (see Eq. 10) computed for each catchment pair i–j, between empirical FDCs and predicted FDCs for either E-HYPE (δ_EHYPE on the y axis) or TNDTK (δ_TNDTK on the x axis) predictions. This figure clearly shows that for 9 out of 11 target sites the geostatistical method TNDTK outperforms E-HYPE in predicting FDCs. Moreover, one of the two sites for which E-HYPE outperforms TNDTK shows nearly the same performance as TNDTK (i.e. the circle is very close to the 1:1 line), while the other one (i.e. site 3675-9001070, highlighted with a black dot in the figure) is associated with the worst performance of TNDTK relative to E-HYPE (see also Sect. 6 on this). Figure 7b reports estimated residual–duration curves ($\widehat{\mathit{\epsilon }}$DCs) for the selected sites. For the sake of representation, we report standardized residuals in the y axis, i.e. residuals divided by the corresponding streamflow quantiles predicted via TNDTK; we referred to TNDTK quantiles for standardization since the real empirical FDC is supposed to be unknown (see cross-validation algorithm illustrated in Sect. 3.1.2). Overall, $\widehat{\mathit{\epsilon }}$DCs show negative values for lower durations and positive values for higher durations (see also Fig. 8). This means that, in Tyrol, E-HYPE tends to overestimate streamflow in wet periods as well as to underestimate streamflows in drier ones relative to the geostatistically predicted FDCs (i.e. TNDTK, see the left panels in Fig. 8). We eventually used the $\widehat{\mathit{\epsilon }}$DC curves, which are estimates of E-HYPE residuals, to assimilate locally available streamflow data into E-HYPE simulated series as illustrated in Fig. 2, obtaining what we termed GAE-HYPE simulations (see right panels of Fig. 8).

Representativeness of simulations (i.e. E-HYPE and GAE-HYPE simulated daily streamflows) is assessed through a Nash–Sutcliffe efficiency of log-flows (LNSE) computed by referring to a recorded streamflow time period of the paired stream gauge (see pairing method adopted in Sect. 4.1). Improvements are obtained with the proposed data-assimilation procedure relative to E-HYPE (Fig. 9). Indeed, we obtained an enhancement of LNSE values of GAE-HYPE simulations relative to the original E-HYPE ones for all the 11 selected sites, which show LNSE increments from −0.462 to 0.594 in the best case (catchment pair IDs: 201236–9608296) and from 0.527 to 0.690 in the worst case (catchment pair IDs: 3675–9001070, see also Table 2). The median on-site LNSE value increases from 0.045 to 0.685, which ultimately underlines the benefits introduced with the proposed method. Figure 9, also, illustrates the impact of geostatistical data assimilation for the two E-HYPE prediction nodes mentioned above (i.e. the one characterized by the best improvements in terms of overall LNSE value, and the one associated with the most limited improvement).

Moreover, looking at the spatial distribution of LNSE values across the 11 selected prediction nodes within the study area depicted in Fig. 10, it is clear how the proposed enhancement strategy benefits from the unbiased estimations of FDCs. In fact, TNDTK shows homogeneous and rather high performance for predicting FDCs (Fig. 10a); also, Fig. 10b and c reveal that the enhancement capabilities of the assimilation procedure are lower for those catchments where E-HYPE performs better (see elementary catchments filled in yellow to green in Fig. 10), whereas the assimilation procedure proves to be powerful when E-HYPE performs worse (see elementary catchments coloured in orange to red), bringing efficiencies from negative to positive values in all cases (from green to blue).

Figure 8 Examples of comparison between observed streamflow series (black dashed lines) and simulated daily streamflows via E-HYPE (red solid lines) and GAE-HYPE (blue solid lines) for two representative sites and a given year, showing two cases for which the geostatistical assimilation procedure resulted in sizeable (a) and limited (b) improvements, respectively.

The assessment of gauge density impacts on the proposed procedure reveals that enhancements are obtained even with the lowest gauge density scenario (i.e. seven gauging stations). Figure 11 displays a clear pattern, showing an improvement in the degree of enhancement associated with an increasing gauge density. Moreover, Fig. 11 shows how the degree of enhancement flattens out in cases in which there are more than approximately 25 gauges available per 10 000 km². Finally, the p values resulting from the Wilcoxon signed-rank tests (see Sect. 3.2) highlight that the assimilation procedure outperforms E-HYPE: the null hypothesis is rejected, with p values always lower than 0.04 %, regardless of the particular stream gauge density scenario.

Figure 9 Scatter diagrams of empirical vs. simulated daily streamflows for either E-HYPE (red dots) and GAE-HYPE (blue dots) for two representative sites, showing the cases in which the data-assimilation procedure respectively produced the largest (a) and smallest (b) degree of enhancement for the study area, respectively.

6 Discussion

This new geostatistical procedure enables practitioners and water resources managers and planners to profit from the wealth of hydrological information, by adjusting open-data products with local observations. We enhanced the streamflow series simulated by macro- and continental-scale rainfall–runoff models at ungauged prediction nodes by assimilating streamflow observations, which are locally available in the region of interest, without having to redo the original hydrological model calculations. This is a recurrent condition since local streamflow data are released under different license terms and policies: some of them could be public and open-access, while some others might not be openly and freely accessible by the broad public. The E-HYPE model obviously lacks storage capacity in the Tyrol region and the proposed approach to enhance the results should be seen as temporary until a new model version accounting for this is available. We do not propose the procedure as a general fix for structurally unsuitable (Beven and Binley, 1992) models, which have been proved to be unfit for either the area of interest or the water problem at hand. For a more sustainable solution, we suggest using another model structure or re-calibrating the model, instead of post-processing the output. However, this procedure makes sense for making a first assessment of water issues in regions where information is otherwise missing, but only macro-scale models are readily available.

Figure 10 Spatial distribution of Nash–Sutcliffe efficiency computed for log-transformed streamflows (LNSEs) at the 11 E-HYPE prediction nodes considered in the study: geostatistically predicted flow–duration curves (FDC TNDTK, a); predicted daily streamflow time series (E-HYPE, b, and GAE-HYPE, c, respectively); the locations of the two sites considered in Figs. 9 and 8 are highlighted with black triangles.

Figure 11 LNSE ratio (see Eq. 12) as a function of stream gauge availability: black dots represent the average of 11 LNSE ratio values, while crosses indicate their 95 % confidence interval.

Our study shows for Tyrol that it is possible to significantly enhance rainfall–runoff simulations resulting from macro-scale, regional or continental hydrological models by geostatistically assimilating (geographically sparse) streamflow observations (see e.g. Figs. 8, 9 and 10); provided that available streamflow series are long enough to obtain a good empirical approximation of the long-term FDC for the site of interest (Castellarin et al., 2013). Indeed, series length of the observed streamflow dataset controls the magnitude of duration extremes (i.e. duration boundary interval), which, in turn, might affect the adopted resampling scheme needed for predicting FDCs at simulation nodes. Nonetheless, in some specific application, very high or very low durations might be of particular interest (e.g. studies focused on flood or drought only); therefore a preliminary investigation on the resampling scheme (e.g. duration extremes, duration intervals, resolution of the points to represent correctly the whole curve) should always be taken into account.

One of the main advantages of the proposed method is that the end user can get enhanced streamflow simulations without any further model calibration or refinement. Even though one could argue that when additional streamflow data become available at neighbouring gauges it should be used for improving the performance of the model at the site of interest, calibrating and validating macro-scale and regional model could be a time-consuming and computationally demanding task. The proposed procedure, instead, is neither computational nor data-intensive, and is implemented only using observed streamflow data and a GIS vector layer with catchment boundaries (see e.g. Fig. 3). The application requires the identification of a suitable regional model for predicting FDC in ungauged basis (see e.g. Fig. 6). However, it has advantages, such as (a) a regional model can be a very informative and useful tool for water resources managers and planners, and (b) the subsequent advantages obtained from the data-assimilation procedure is transferred downstream in the entire regional river network (see Fig. 11).

One important limitation of the proposed method is that, once a target prediction node is considered, any given simulated streamflow value is associated with a single duration, which corresponds to a particular estimated residual, which will be used in turn for correcting the streamflow value itself (see Figs. 2 and 8). Essentially, this means that the volumes from E-HYPE are discarded while the sequencing of E-HYPE simulations is retained. Moreover, this algorithm cannot possibly account for seasonal (or interannual) modifications in the hydrological behaviour of the catchment. Indeed, as shown in the time series comparison in Fig. 8, when the geostatistical prediction of FDCs is unreliable, the assimilation procedure reflects such low accuracy (i.e. the procedure fails to correctly capture high-flow and low-flow regimes; see e.g. the resulting FDC from GAE-HYPE simulations at the catchment pair 3675–9001070 in Fig. 8), propagating this bias throughout the whole simulated series. Finally, designing a theoretical framework that combines statistical data-driven approaches with deterministic process-driven ones is seen by many as the correct way for tackling the “prediction in ungauged basins” (PUB) problem and further advancing the scientific research in this area (Di Prinzio et al., 2011). We believe that our geostatistical data-assimilation procedure for macro-scale hydrological models is one example in this direction. Future analyses will focus on the relaxation of the main limitation of the approach (i.e. the incorporation of seasonal patterns in a data-assimilation procedure) and on the extension of its applicability to anthropologically altered streamflow regimes.

7 Conclusions

This research work focuses on the development of an innovative method for enhancing streamflow series simulated by macro-, continental-, and global-scale rainfall–runoff models by means of a geostatistical prediction of model residuals. We focus on Tyrol as study region and E-HYPE (European-HYdrological Predictions for the Environment, from the Swedish Meteorological and Hydrological Institute, SMHI) as a macro-scale hydrological model; nevertheless, the geostatistical data-assimilation procedure is general and can be applied to simulated streamflow series coming from other macro-scale rainfall–runoff models. The proposed data-assimilation procedure utilizes streamflow data that are locally available for the area of interest, which were not considered in the implementation of the macro-scale hydrological model; it (1) adopts top-kriging for regionally interpolating empirical period-of-record flow–duration curves (FDCs) that can be constructed from locally available streamflow data; (2) constructs residual–duration relationships at any prediction node in the study region where simulated streamflow series are available, by comparing FDCs resulting from geostatistical interpolation (top-kriging) and rainfall–runoff simulation (E-HYPE); and (3) uses the residual–duration curve to enhance macro-scale simulated streamflows.

The cross-validation tests of the proposed approach with different scenarios of streamflow data availability show the significant advantages of geostatistical data assimilation even for very low stream-gauging network densities (i.e. ca. 1 gauge per 2000 km²). It can become a stand-alone numerical tool to be used for enhancing results from macro-scale models anywhere along the stream network of a given region. Potential applications are envisaged for a variety of water resources management and planning problems that require accurate streamflow series (e.g. regional assessment of hydropower potential, habitat suitability studies, surface water allocation, civil protection management strategies, climate change trends, safety of river structures). Future analyses will address the main limitation of the proposed geostatistical data-assimilation procedure, aiming to incorporate observed seasonal and inter-annual variations of the hydrological behaviour of the study region as well as other hydrological features (e.g. baseflow index or peak-flow data) into the geostatistical regionalization of model residuals.