Evaluation of a stochastic weather generator for long-term ensemble streamflow forecasts

ABSTRACT

Resampling historical time series remains one of the main approaches used to generate long-term probabilistic streamflow forecasts, while there is a need to develop more flexible approaches taking into account non-stationarities. One possible approach is to use a modelling chain consisting of a stochastic weather generator and a hydrological model. However, the ability of this modelling chain to generate adequate probabilistic streamflows must first be evaluated. The aim of this paper is to compare the performance of a stochastic weather generator against resampling historical meteorological time series in order to produce ensemble streamflow forecasts. The comparison framework is based on 30 years of forecasts for a single Canadian watershed. Forecasts resulting from the two methods are evaluated using the continuous ranked probability score (CRPS) and rank histograms. Results indicate that while there are differences between the methods, they nevertheless perform similarly, thus showing that weather generators can be used as substitutes for resampling the historical past.

KEYWORDS:

• stochastic weather generator
• resampling
• ensemble streamflow forecast
• long term

Editor A. Fiori Associate editor A. Requena

1 Introduction

Long-term streamflow forecasting (1–12 months ahead) is important to various water resource sectors, such as hydropower generation and optimal reservoir operation (Georgakakos 1989, Hamlet et al. 2002, Markoff and Cullen 2008). Consequently, providing reliable long-term streamflow forecasts has long been a concern for the hydrological science community.

A main issue with long-term streamflow forecasts involves properly incorporating uncertainties that grow larger as the forecast lead time increases (Georgakakos 1989). Deterministic forecasts have long been used, but have proven to be inadequate for long-term forecasts, and this has led researchers to focus on probabilistic approaches to properly encompass the large uncertainties associated with distant forecasts. An ensemble prediction system (ESP; Cooke 1906, Weather Bureau Brier 1944, Murphy and Winkler 1984, Day 1985) is a collection of two or more forecasts over the same future horizon, which represents the various probability states and uncertainties associated with the forecasting process (Sivillo et al. 1997). Various studies have demonstrated the benefit of using ESP in reservoir management and in the decision-making process (Faber and Stedinger 2001, Kim et al. 2007, Zhao et al. 2013). ESP is now widely used in research and operational forecasts around the world (Cloke and Pappenberger 2009).

The methods for making probabilistic forecasts can be classified into three main categories, namely numerical weather prediction systems (NWPs), resampling methods and stochastic weather generators (WGs). NWPs, obtained by assimilating current atmospheric and oceanic conditions into a coupled numerical mathematical model, are the most accurate tool for forecasting future weather conditions (Pielke 2013). However, forecasts from NWPs typically have skill up to a 6-day horizon for low precipitation (less than 2 mm) and for up to 4 days for high precipitation (between 2 and 10 mm) (Chen et al. 2010). NWPs are therefore not well suited for hydrological forecasts with longer time horizons as they cannot adequately represent long-term forecast uncertainty (Coulibaly 2003, Cloke and Pappenberger 2009). Accordingly, long-term forecasts typically favour resampling of the historical past as well as the use of WGs. Seasonal weather forecasts from coupled numerical models, both calibrated and uncalibrated, are currently receiving a lot of attention with the potential of outperforming ESP depending on season, region and lead time (Yuan 2016, Bennett et al. 2017, Arnal et al. 2018). However, resampling remains by far the most commonly used operational approach when it comes to generating long-term forecasts and estimating forecast uncertainties (Todorovic and Woolhiser 1975, Young 1994, Leander and Buishand 2009, King et al. 2014).

Resampling approaches were developed around 1990 to alleviate the problems inherent in parametric models (Lall and Sharma 1996), with one main such problem being the assumption of normality, which is not realized for most hydrological variables. This assumption requires data transformation, which is a primary cause of bias (Lall and Sharma 1996, Sharma et al. 1997, Prairie et al. 2006). In addition, parametric methods often have difficulties representing the complex nonlinear relationships observed between hydrometeorological variables. This weakness is what led to the development of nonparametric resampling methods, such as the index sequential method (ISM), the kernel-based approach and k-Nearest Neighbors (K-NN) bootstrapping methods (Kendall and Dracup 1991, Sharma and Sharma 1997, Salas and Lee 2009, Sharifazari and Araghinejad 2015, Li et al. 2017, Khaki et al. 2018). Although widely used thanks to their relative ease of implementation, most resampling approaches do, however, suffer from the same potential drawbacks.

Resampling methods are based on the historical record, and any deficiency in this record will affect the forecast quality. The methods are also unable to generate weather values larger (or smaller) than observed values (Rajagopalan and Lall 1999), which limits their ability to properly assess uncertainty related to rare events. The length of the historical record also brings limitations when it comes to appropriately sampling forecast uncertainty. Finally, non-stationarities due to internal climate variability and anthropogenic change call into question the validity of using the historical past for accurate long-term streamflow forecasts (Murphy and Winkler 1984, Clark et al. 2004a, 2004b). Despite all the different methods and techniques that have been proposed over the years, additional research is still needed to fully address these issues (Lall 1995, Prairie et al. 2006, 2007, Lee et al. 2010).

Stochastic WGs were developed in the 1980s with the aim of generating long-enough synthetic weather time series with specific statistical properties to be used in hydrologic models and risk assessment studies, and to extend and simulate weather time series at locations with no observed data (Wilks and Wilby 1999). The long time series of WGs have been used in many applications such as flood type classification (Turkington et al. 2016) and simulation of extreme floods (Okoli et al. 2019, Výleta and Valent 2019, Winter et al. 2019, Sikorska-Senoner et al. 2020) up to the probable maximum flood (Chen et al. 2015), as well as in many climate change studies (Semenov and Barrow 1997, Keller et al. 2017, Steinschneider et al. 2019, Ahn 2020). Many other studies have focused on improving WGs to represent higher statistical moments of both univariate and multivariate statistics of observed weather data (Hayhoe 2000, Hansen and Mavromatis 2001, Kyselý and Dubrovský 2005, Brissette et al. 2007, Chen et al. 2012a). Other studies have also looked into using WGs as constituting a downscaling tool in climate change impact studies because of the relative ease with which their parameters can be modified to represent climate variability (Semenov and Barrow 1997, Kilsby et al. 2007, Zhuang et al. 2016, Keller et al. 2017). Despite this advantage, only a few studies have looked at the potential of using stochastic WGs for long-term streamflow forecasting (Li et al. 2013, Breinl et al. 2015, Shield and Dai 2015, Breinl 2016), and an even smaller number implemented a WG into streamflow forecasts. Of note, Šípek and Daňhelka (2015) used a WG for ESP based on a limited number of observed years selected on the basis of large-scale climate indices, and Hwang et al. (2011) generated several precipitation scenarios to study the impact of uncertainty on streamflow simulations. Both studies pointed to the potential of WGs for improving ESP.

WGs have the potential to overcome the main drawbacks of resampling approaches and, notably, to take into account non-stationarities. Resampling approaches are much less flexible due to the limited number of years available. Before moving on to the evaluation of approaches dealing with non-stationarity, a necessary first step is to evaluate the ability of a stochastic WG for generating ensemble streamflow forecasts. Therefore, the goal of this paper is to compare the performance of ensemble streamflow forecasts built from resampling the historical record against that of a stochastic WG when coupled with a hydrological model to provide long-term probabilistic streamflow forecasts.

The remainder of this paper is divided into four sections. The watershed and data are first described, followed by the methodology employed. Results are then presented and discussed in the last two sections.

2 Description of watershed and data

The Lac Saint-Jean watershed was selected as a test case for this study. It is located in the eastern Canadian province of Quebec (Fig. 1). It has a surface area of 73 800 km² and consists of four main rivers flowing into a 1000 km² lake (Lac Saint-Jean), which is regulated at its outlet by a large hydropower dam. Rio Tinto, the world’s largest producer of aluminium, operates six power plants on this watershed, with an annual average capacity of 2000 megawatts, feeding a large aluminium plant. Economic development around Lac Saint-Jean is therefore closely linked to its hydropower potential (Dibike and Coulibaly 2005). Overall, the watershed is sparsely inhabited, and its land cover mostly consists of a homogeneous boreal forest.

Figure 1. Lac Saint-Jean watershed. W: West; N: North

Daily precipitation, daily maximum and minimum temperatures are provided in gridded format by the Natural Resources Canada (NRCAN) dataset. NRCAN is the only high-resolution Canada-wide daily-scale precipitation and temperature gridded dataset. It provides daily values with a 10-km spatial resolution (300 arcseconds). The NRCAN dataset is interpolated from the Canadian network of gauge surface stations using trivariate thin-plate smoothing splines, providing continuous functions as a function of latitude, longitude and ground-surface elevation (Hutchinson et al. 2009). The study period covered in this work spans the years 1950 to 2010. Daily scale naturalized streamflow data for Lac Saint-Jean was provided by Rio-Tinto, and covers the same time period. Figure 2 presents the annual cycle (and variability) of the hydrometeorological variables used in this study. Temperature and streamflow values are characterized by a strong seasonal pattern typical of a northern latitude climate, whereas the precipitation values show a weaker cycle with increased mean monthly totals in the summer and early fall.

Figure 2. Annual cycle of key hydrometeorological variables used in this study. (a) Minimum and maximum temperatures; (b) precipitation; (c) streamflow

3 Methodology

The aim of this project is to compare the performance of WGs against the benchmark method of equiprobable resampling for long-term streamflow forecasting. To this end, the framework for this comparison consists of three main steps. In the first step, for each forecast date, ensemble weather forecasts are generated from both approaches. In the second step, these weather forecasts are used with a hydrological model to generate ensemble streamflow forecasts. In the last step, the forecasts from both methods are compared against observed past conditions. Details are provided below.

Both methods are evaluated in hindcast mode over the 30-year, 1980–2009 period. The 1950–1979 time period is used as the initial historical record. Over the evaluation period, 1-year-ahead forecasts are made 12 times per year, on the first day of each month. The evaluation period therefore comprises 360 1-year forecasts (30 years × 12 forecasts). Similarly to what would be done in the real world, the recalibration of the hydrological model and WG is performed at the beginning of every new year, adding the preceding year to the historical record. Figure 3 presents the methodological framework.

Figure 3. Methodological framework

3.1 Ensemble weather forecast: approach 1 – resampling

In this study, the simple method of resampling the past observed climatology data is used as the benchmark approach. As discussed earlier, resampling past climatology is the most commonly used approach for long-term ensemble weather forecasting. This work uses equiprobable resampling with no reshuffling between precipitation and temperature years. The number of members in each of the ensemble weather forecasts is therefore equal to the number of years in the historical period.

3.2 Ensemble weather forecast: approach 2 – stochastic WG

The WG used in this study is WeaGETS, the uni-site multivariate WG (Chen et al. 2012a). WeaGETS is based on the original work of Richardson (1981) and Richardson and Wright (1984). It first generates precipitation occurrence using a first-, second- or third-order Markov chain. On wet days, precipitation is then generated using either an exponential or a gamma distribution. Minimum and maximum temperatures are generated conditionally on the wet/dry day status. Serial correlation and cross-correlations between precipitation and temperature are preserved using a first-order autoregressive process. Inter- and intra-annual variability is preserved using the spectral approach of Chen et al. (2010). All WeaGETS parameters are first computed on a monthly basis, and are then interpolated to the daily scale with Fourier harmonics, as proposed by Richardson (Richardson 1981) More details can be found in Chen et al. (2010, 2012a). For this study, a first-order Markov chain is used to produce the precipitation occurrence, and the two-parameter gamma distribution is used for estimating precipitation amounts. The gamma distribution has been used in many studies to model precipitation in a WG (Chen et al. 2016). It performs generally well but often fails at generating extreme precipitation that follows a distribution with a heavier tail, such as the generalized extreme value distribution (GEV) or Generalized Pareto distributions. Accordingly, several studies have compared various distribution functions to be used in WGs (Li et al. 2013, Breinl 2016, Breinl et al. 2017, Li and Shi 2019). These studies typically find that the more complex distribution functions, with three or more parameters, outperform the simpler ones, especially for larger quantiles of the observed distributions. However, for a study of long-term hydrological forecasts, the gamma distribution was deemed adequate. In addition, for large catchments such as the one under study, the impact of extreme convective precipitation on part of the catchment is lessened by the averaging of the streamflow-generating processes at the watershed scale.

Like most WGs, WeaGETS can generate an extremely large number of years (limited by computer storage/memory); 500 members were generated for each forecast. We used multiple randomly selected samples of N years to explore the impact of the number of members, but no significant difference was observed. Therefore, to allow for a fair comparison, for all forecasts the number of WG members is always identical to the number of years available in the historical past. Prior to using the WG, precipitation, minimum and maximum temperatures were averaged using all NRCAN grid points within the catchment boundaries. Despite the size of the catchment, a unisite WG (as opposed to a more complex multisite WG) was chosen because of the catchment’s very uniform topography and land cover (boreal forest). Despite the high spatial resolution of the NRCAN gridded dataset, the underlying station network over this high-latitude watershed is very sparse, and did not justify the use of a more complex multisite WG. In addition, Chen et al. (2016) showed that using a multisite WG on this same catchment did not result in any improvement in hydrological modelling, when compared to its simpler counterpart. In general, overparameterized methods may not always perform well, and the simpler methods are often preferred by operational agencies.

3.3 Hydrological model

The hydrological model used in this study is the lumped conceptual Service Hydrométéorologique Apports Modules Intermédiaires (HSAMI) model. HSAMI was developed by Hydro-Québec, and is used to make hourly and daily operational forecasts for more than 80 watersheds with surface areas ranging between 160 and 70,000 km². It has also been used in many Canadian hydrological studies, including climate change impact studies (Minville et al. 2008, Chen et al. 2011), multi-model simulations (Arsenault et al. 2015), model calibration experiments (Arsenault et al. 2013) and regionalization studies (Arsenault and Brissette 2014a, 2014b). HSAMI is a conceptual rainfall–runoff lumped model with 23 parameters. Two parameters account for evapotranspiration, 6 parameters for snowmelt simulation and 15 parameters for vertical and horizontal water movement. Vertical water flows are estimated by four interconnected linear reservoirs (snow on the ground; surface water; saturated and unsaturated zones). Horizontal routing is carried out by two unit hydrographs and one linear reservoir for low flows. The minimum required inputs are daily scale precipitation, maximum and minimum temperature, and flow discharge for the calibration period. If available, cloud cover fraction and snow water equivalent can also be used as inputs. The calibration of HSAMI was performed automatically using the covariance matrix adaptation evolution strategy (CMA-ES) algorithm (following the work of Arsenault et al. (2013)), with the Nash-Sutcliffe criterion used as an objective function. The calibration process followed the procedure outlined in Arsenault et al. (2018), and bypassed the validation step.

To avoid any bias resulting from the hydrological model, and to forgo the complex assimilation process, simulated streamflows were used instead of observations. Doing so results in a “perfect” hydrological model and eliminates all uncertainties resulting from the assimilation process and hydrological modelling biases. This ensures that all differences between the resampling and WG ensemble streamflow forecasts are entirely due to differences in the generation of the ensemble weather forecasts.

3.4 Evaluation procedure

The challenge for all probabilistic forecasts is to provide an unbiased ensemble with the proper dispersion. Despite the steadily growing body of literature on ensemble streamflow forecasting, the number of evaluation metrics remains relatively small. The main probabilistic score broadly used in the meteorological and hydrological community for ensemble verification remains the continuous ranked probability score (CRPS; Matheson and Winkler 1976). The CRPS generalizes the mean absolute error (MAE) to the case of probabilistic forecasts.

If F is the predictive cumulative distribution function (CDF) and y is the observation, the CRPS can be defined by:

(1) $\mathrm{C}\mathrm{R}\mathrm{P}\mathrm{S}\left(F,y\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\left[F\left(t\right)-H\left(t-y\right)\right]}^2\mathrm{d}\mathrm{t}$(1)

where $H\left(t-y\right)$ is the Heaviside function. It is 0 when $t<y$, and 1 otherwise (Hersbach 2000, Toth et al. 2003).

A perfect CRPS value is zero (all members of the ensemble predicting the same value, equal to the observation) and, therefore, the aim of a probabilistic forecast is to minimize the CRPS value. The CRPS allows for discrimination between different ensemble forecasts.

Another important part of ensemble verification is related to the dispersion of the forecast. To this end, the rank histogram (Talagrand diagram) has been widely used as a graphical evaluation method to visually assess the bias and spread of an ensemble forecast. The rank histogram measures how well the ensemble forecast represents the probabilistic distribution of observations (Hamill 2001). To plot the rank histogram, the forecast is typically separated into n + 1 bins (with n being the number of members in an ensemble), and a count is tabulated for the number of times observations fall within each bin. For a perfect ensemble forecast, the rank histogram will be flat since the probability of the observation falling within each bin would be equal. Alternatively, the number of bins can be related to the CDF of the ensemble. An asymmetrical rank histogram is the result of a biased forecast, whereas concave or concave shapes represent, respectively, an under- or over-dispersed forecast (Hersbach 2000, Hamill 2001, Wilks 2006, Zalachori et al. 2012). To measure the deviation from the flat and uniform rank histogram, the flatness coefficient (δ) described in Equation (2)(2) $\delta =\frac{1}{n+1}\sum _{z=1}^{z=n+1}\left|f\left(z\right)-y\right|$(2) is used:

(2) $\delta =\frac{1}{n+1}\sum _{z=1}^{z=n+1}\left|f\left(z\right)-y\right|$(2)

where $f\left(z\right)$ is the relative frequency in rank $z,y=\frac{1}{n+1}$ is the theoretical relative frequency, and $n$ is the number of ensemble members.

In a perfectly flat rank histogram, the value of δ is equal to 0 (Velázquez et al. 2010).

4 Results

Prior to using the modelled streamflows as pseudo-observations, it is nonetheless important to ensure that HSAMI performs accurately over the Lac Saint-Jean watershed. Figure 4 shows the mean hydrograph derived from observation and modelled streamflows over the 1950–1980 period. The performance of HSAMI is excellent (Nash-Sutcliffe of mean daily hydrograph is 0.9885) and justifies the use of a conceptual lumped model, despite the large size of the Lac Saint-Jean watershed. This performance also justifies the use of a single-site WG, rather than its more complex multi-site version (Brissette et al. 2007, Chen et al. 2014).

Figure 4. Mean annual hydrograph of simulated and observed inflows into Lac Saint-Jean

As was shown in the methodology framework, the ESP forecasts are made by two methods at each forecast date and over the entire forecast period. CRPS values are calculated daily for each forecast. Figure 5 shows the ensemble streamflow forecasts made on 1 January 1980, and the associated CRPS values for both forecast approaches. The larger apparent spread observed for the WG is because all 500 generated members were used to plot this figure. For the comparison of all forecasts, the first N members (between 30 and 59 depending on the hindcast year) of the generated ensemble were chosen, as discussed earlier. Otherwise, the two approaches provide forecasts with similar CRPS values.

Figure 5. One-year ensemble streamflow forecasts made on 1 January 1980 using resampling (left) and weather generator (right) approaches. The black line represents observations. There are 30 members for resampling (1950–1979) and 500 members for the weather generator (first row). The CRPS values are shown in the second row

Figure 6 presents mean monthly CRPS values aggregated across all forecasts. Twelve 1-year forecasts are made for every year, one on the first of each month. The mean monthly CRPS value is computed for each year by taking the mean of each daily CRPS (30 days × 12 forecasts). The box plots presented in Fig. 6 are made with the mean monthly CRPS values computed for each forecasted year (1980–2009). Values are then used to construct the box plots. The solid box plot rectangles represent the interquartile spread (25th and 75th quantiles), with the median near the middle, while the whiskers show the extent of the values. Outliers are identified with the + sign.

Figure 6. Box plot results of monthly CRPS based on ESPs forecasted by resampling and the weather generator method on 1 January, 1 April, 1 July and 1 October

These results indicate that the two forecasting approaches perform similarly in all cases. The largest CRPS values are observed during the spring flood in April, May and June, when uncertainty and streamflow values are at their maximum, whereas the lowest values are related to the winter low flows, where uncertainty and streamflow values are relatively small. The impact of the forecast date on CRPS values is quite clear when looking at the 1 April forecast, which results in lower uncertainty in the spring flood. To better outline differences between the two forecasting approaches, Fig. 7 presents the 30-year average of all mean monthly CRPS values for all 12 forecast dates.

Figure 7. Mean monthly CRPS of ensemble streamflow forecasts using the resampling (top) and weather generator (bottom) approaches

As was the case for Fig. 6, the mean monthly values are obtained by averaging daily CRPS values within each month, for all 30 years of the hindcasting period and all 12 forecast dates. Therefore, each mean monthly value represents the average for 360 forecasted days. The vertical axis presents the date of the forecast, and the horizontal axis shows the month of the forecast. Accordingly, the diagonal is distinctly darker since it corresponds to the first month of each forecasting date, when uncertainty is lowest. The general pattern of the mean monthly CRPS is very similar for both approaches. The lowest mean monthly CRPS values are consistently found in the winter, when the streamflows are consistently small. Mean monthly CRPS values are highest during the snowmelt period (April to June), when streamflow values and variability are highest. To better outline differences between the two approaches, Fig. 8 presents the difference between the estimates of mean monthly CRPS values.

Figure 8. Difference between the mean monthly continuous ranked probability score (CRPS) values from the resampling and weather generator approaches. Red indicates that the resampling approach has a lower mean monthly CRPS, while blue indicates the reverse

The difference between the two approaches is minimal from January to March. The resampling forecast ensemble performs generally better in April and May, especially for longer lead times (above the diagonal). The WG approach generally gives lower mean monthly CRPS values in the summer, but consistently larger ones in November and December. Nonetheless, differences remain small, with the exception of May forecasts issued 7–12 months in advance.

To further compare the two approaches, rank histograms are presented in Figs. 9 and 10. Each histogram is composed of 10 800 forecasted days (30 days per month × 12 1-year forecasts × 30 years). The histogram flatness value (Equation 3) is presented within each histogram figure.

Figure 9. Plotted rank histograms based on monthly ESP forecasts using the resampling method

Figure 10. Plotted rank histograms based on monthly ESP forecasts using the weather generator method

A constant number of 30 ensemble members were used for all forecasts. The WG randomly generated 30 years for each forecasting date. For resampling, when the number of historical years was larger than 30, the ensemble members were chosen randomly from amongst all available years.

Overall, all histograms are relatively flat. All forecasts for both approaches are relatively well dispersed. There is, however, one main difference between the two approaches. Despite having similar mean CRPS values (Figs. 6 and 7), the rank histograms for the WG approach display a small but consistent negative bias for the winter months (November to March), with flatness values consistently higher than for the resampling values. The two approaches otherwise perform similarly during all other months, although resampling results in a larger positive bias during the post-flood period, from May to November.

Finally, Fig. 11 presents the rank histogram of all 1-year monthly forecasts made from 1 January 1980 to 31 December 2009. The resampling rank histogram is flat, with a very slight positive bias. Taken globally, the results indicate that the WG’s performance is not nearly as good as resampling, as it is negatively biased and slightly under-dispersed.

Figure 11. Global rank histogram pooled from all 1-year monthly forecasts for the resampling (left) and weather generator (right) approaches

5 Discussion

The results presented are based on 360 1-year forecasts issued on the first of each month, from January 1980 to December 2009. Overall, the results, presented in terms of CRPS and rank histograms, show that the two approaches perform similarly. In terms of the CRPS metric, the two approaches are globally equivalent. The WG method generates slightly better forecasts just before the flood (January to April), whereas resampling is generally slightly better for the forecasts issued after the flood. The performance is also similar when looking at bias and dispersion from rank histograms, although the WG forecasts are globally slightly under-dispersed and negatively biased. This shows the importance of using more than one metric to evaluate forecast quality (Hamill 2001, Renner et al. 2009, Demargne et al. 2010). The differences between the two approaches remain small, and become clearer after all forecasts are pooled together (Fig. 11). The reasons for the differences are difficult to explain. In principle, stochastic WGs are built to reproduce uni- and multivariate statistics of observed time series. The WG used in this work has been extensively tested for its ability to reproduce characteristics of precipitation and temperature series (Caron et al. 2008, Li et al. 2010, Chen et al. 2014), and especially for its ability to correct for the underestimation of intra- and inter-annual variance (Chen et al. 2010). To better understand the reasons behind this behaviour, we plotted rank histograms of accumulated precipitation at 3, 6, 9 and 12 months, as shown in Fig. 12.

Figure 12. Rank histograms of accumulated precipitation at 3, 6, 9 and 12 months from the weather generator (left) and resampling (right) approaches

This shows that the WG precipitation is under-dispersed and slightly negatively biased, which tracks well with the results observed for streamflow as displayed in Fig. 11. This appears to result from an underestimation of the temporal autocorrelation of precipitation by the WG. The structure of WeaGETS does not specifically prescribe the temporal autocorrelation of precipitation, and all other requirements (e.g. preserve wet/day series, cross-correlations with temperature) are not sufficient to induce a temporal autocorrelation. The variance correction scheme of Chen et al. (2010) preserves inter-annual autocorrelation but does not address it at the intra-annual scale, which is most important for streamflow forecasts. Figure 13 presents CRPS of monthly precipitation for both the WG and resampling.

Figure 13. Mean monthly CRPS box plots for the 1-year, 1 January forecast over the entire 30-year hindcast period, 1980–2009

The results show that the WG CRPS values are slightly better than those from resampling, which is consistent with the results of Fig. 6 for streamflows. A similar exercise was performed for temperature (results not shown) and no noticeable difference was found between resampling and WG outputs. In the case of temperature, the WG uses an autoregressive scheme which specifically mandates temporal autocorrelation of the generated time series. It therefore appears that most of the observed differences can be tracked to the underestimation of intra-annual precipitation autocorrelation by the WG.

The performance of WG outputs in terms of driving environmental models (such as a hydrology model) has been less studied in the literature, even though this is one of the main reasons why WGs were developed in the first place. Such studies performed with hydrological models indicate that WG series provide realistic streamflow series for a wide variety of conditions. The small differences observed between WGs and original observed series however can be increased due to the impact of hydrological model (Khalili et al. 2011, Li et al. 2013). The results presented in this work are consistent with these findings, in that WG streamflow forecasts appear to perform similarly to those from resampling, but with some relatively minor differences that are difficult to track back to the WG precipitation and temperature time series. The hydrological model acts as a complex non-linear filter of precipitation and temperature, and small changes in one or both series may result in somewhat larger differences. If the WG time series were statistically identical to the observed series, one would expect the two streamflow ensembles to perform identically. Any differences can therefore be tracked back to minor inadequacies in the WG series characteristics.

The presence of biases in forecasts is typically rooted in uncertainties due to a combination of various input, output, model structure and parameter uncertainties (Krzysztofowicz 2002, Han et al. 2007, Li et al. 2009). Such biases have been noted in resampled time series in other studies (Seo et al. 2006, Schaake et al. 2007, Bogner and Kalas 2008). In particular, biases in long-term forecasts may be due to the non-stationarities present in historical time series (Ceola et al. 2014, Meng et al. 2019) due to internal climate variability and long-term trends linked to anthropogenic forcing.

Non-stationarity in hydro climatic time series is now well documented in many regions of the world (Westra et al. 2013). Non-stationarity raises important questions about equiprobable historical resampling. Using non-equiprobable resampling is a complex endeavour, for which there is currently no easy solution. It is, however, relatively easy to incorporate non-stationarity into a stochastic WG. A few studies have provided frameworks that could be applied to the problem of generating long-term ensemble streamflow forecasts (Semenov and Barrow 1997, Chen et al. 2012b, Jones et al. 2016, Keller et al. 2017, Li Liu et al. 2017). For example, the variability structure of historical meteorological time series at seasonal, monthly and daily scales could be preserved, while mean annual and seasonal values could be conditioned on large-scale teleconnection indices or by weighting years differently, with the more recent ones being more important. Comparatively, resampling approaches are much less flexible due to the limited number of years available. If the historical record is long enough, non-equiprobable resampling offers the possibility of selecting a subset of years judged more relevant based on criteria such as similarity to current/recent conditions (e.g. wet, warm) or based on a limited number of climate indices (e.g. El Niño Southern Oscillation). Such selection is, however, far from simple, and the “less relevant” years left out of the process may still contain important information on internal climate variability that would be useful to long-term forecasts. A WG could be used to generate stochastic variability around a small subset of years based on the similarity to present climate indices, as was done by Šípek and Daňhelka (2015), for example. In practice, however, correlations between climate indices are typically small in many regions of the world.

In any event, approaches dealing with climate non-stationarity for long-term streamflow forecasting must recognize the two main components of climate change: internal variability and anthropogenic forcing. Internal or unforced variability comes from the chaotic behaviour of the climate system, and impacts the regional climate at seasonal to multidecadal scales. Internal variability is most often represented by patterns of precipitation and temperature anomalies at local and regional scales. Anthropogenic forcing results from the increase in greenhouse gas emissions) and results in slow but steady changes in the climate. For temperature, anthropogenic forcing already displays a clear signature, with warming exceeding 2°C annually at high latitudes, and has exceeded that due to internal variability in most parts of the world (Lehner et al. 2017), whereas for precipitation, internal variability is expected to remain the dominant variability component until later this century (Martel et al. 2018). Incorporating temperature and precipitation non-stationarity into long-term streamflow forecasts should therefore follow a different process for each variable.

6 Conclusion

This paper presents a comparison of two methods for generating long-term ensemble streamflow forecasts. Both methods make use of a hydrological model to transform precipitation and temperature time series into streamflows. The first approach resamples historical time series, while the second makes use of a stochastic WG calibrated using the same historical series. The main goal of this paper is therefore to validate the use of a stochastic WG for ensemble streamflow forecasting. To this end, forecasts made from the combination of a stochastic WG and a hydrological model are compared to those obtained from equiprobable resampling of past meteorological variables. The evaluation period consisted of 1-year-long forecasts issued on the first day of each month over a 30-year hindcast period (1980–2009). Forecasts resulting from the two methods are evaluated with respect to CRPS and rank histograms.

Results indicate that while there are differences between the two methods, they largely perform similarly, which thus indicates that WGs can be used as substitutes for resampling the historical past. Potential approaches to modify WGs to take into account non-stationarities are discussed.

Acknowledgements

This research project was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through a Collaborative Research and Development (CRD) Grant, and by Hydro-Québec, Rio Tinto and the Ouranos Consortium on Regional Climatology and Adaptation to Climate Change. The authors acknowledge Korbinian Breinl (Vienna University of Technology) and an anonymous reviewer for their constructive comments on this article.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) [CRNSG # RDCPJ 428468 – 11].