Articles | Volume 22, issue 6
Research article
29 Jun 2018
Research article |  | 29 Jun 2018

Seasonal streamflow forecasting in the upper Indus Basin of Pakistan: an assessment of methods

Stephen P. Charles, Quan J. Wang, Mobin-ud-Din Ahmad, Danial Hashmi, Andrew Schepen, Geoff Podger, and David E. Robertson

Timely and skilful seasonal streamflow forecasts are used by water managers in many regions of the world for seasonal water allocation outlooks for irrigators, reservoir operations, environmental flow management, water markets and drought response strategies. In Australia, the Bayesian joint probability (BJP) statistical approach has been deployed by the Australian Bureau of Meteorology to provide seasonal streamflow forecasts across the country since 2010. Here we assess the BJP approach, using antecedent conditions and climate indices as predictors, to produce Kharif season (April–September) streamflow forecasts for inflow to Pakistan's two largest upper Indus Basin (UIB) water supply dams, Tarbela (on the Indus) and Mangla (on the Jhelum). For Mangla, we compare these BJP forecasts to (i) ensemble streamflow predictions (ESPs) from the snowmelt runoff model (SRM) and (ii) a hybrid approach using the BJP with SRM–ESP forecast means as an additional predictor. For Tarbela, we only assess BJP forecasts using antecedent and climate predictors as we did not have access to SRM for this location. Cross validation of the streamflow forecasts shows that the BJP approach using two predictors (March flow and an El Niño Southern Oscillation, ENSO, climate index) provides skilful probabilistic forecasts that are reliable in uncertainty spread for both Mangla and Tarbela. For Mangla, the SRM approach leads to forecasts that exhibit some bias and are unreliable in uncertainty spread, and the hybrid approach does not result in better forecast skill. Skill levels for Kharif (April–September), early Kharif (April–June) and late Kharif (July–September) BJP forecasts vary between the two locations. Forecasts for Mangla show high skill for early Kharif and moderate skill for all Kharif and late Kharif, whereas forecasts for Tarbela also show moderate skill for all Kharif and late Kharif, but low skill for early Kharif. The BJP approach is simple to apply, with small input data requirements and automated calibration and forecast generation. It offers a tool for rapid deployment at many locations across the UIB to provide probabilistic seasonal streamflow forecasts that can inform Pakistan's basin water management.

1 Introduction

The Asian Development Bank rates water security in Pakistan as “hazardous” (the lowest of five classes), ranking it 46th out of 48 countries in the Asia–Pacific region, with only Kiribati and Afghanistan ranked lower (Asian Development Bank, 2016). Other studies confirm Pakistan's relatively high levels of exploitation of river flows and groundwater, associated water stress and resultant exposure to climate change (Döll et al., 2009; Wada et al., 2011; Schlosser et al., 2014; Kirby et al., 2017). Given the high demands on the main water source, the Indus River, its year-to-year flow variability has a significant impact on security of supply in the Indus Basin Irrigation System (IBIS) of Pakistan. Better management outcomes could be achieved if a reliable understanding of Kharif (April–September, “all Kharif”) water availability at the beginning of the season were available. This would improve IBIS water allocation planning, a critical need given the highly seasonal flows ( 80 % annual flow occurs in the Kharif season), limited storage capacity (10 % of inflows) and increasing water demand for agriculture and energy production. Thus we assess methods for providing seasonal streamflow forecasts for the two largest water supply dams, Tarbela (on the Indus) and Mangla (on the Jhelum), in the upper Indus Basin (UIB) of Pakistan.

Seasonal streamflow forecasts can be a valuable source of information for water resource managers (Chiew et al., 2003; Anghileri et al., 2016), with both statistical and dynamical forecasting approaches developed and implemented internationally (Yuan et al., 2015). Sources of seasonal streamflow predictability come from initial hydrological or antecedent conditions (e.g. water held in storage in a catchment, in the soil, as ground water, in surface stores or as snow/ice) and also from the skill of seasonal climate forecasts (Bennett et al., 2016; Doblas-Reyes et al., 2013; Li et al., 2009; Shukla and Lettenmaier, 2011; Shukla et al., 2013; van Dijk et al., 2013; Koster et al., 2010; Wood et al., 2015; Yossef et al., 2013). Statistical approaches relate antecedent catchment conditions and/or climate indices to streamflow using techniques such as multiple linear regression (Maurer and Lettenmaier, 2003). Statistical approaches require predictor–predictand records of sufficient length to determine robust relationships, stationarity in the relationships and rigorous cross-validation to avoid over-fitting or an inflated skill assessment (Robertson and Wang, 2012; Schepen et al., 2012). Dynamical approaches use hydrological models initialised with observed inputs up to the beginning of the forecast season (to account for antecedent conditions) that can be driven either by historical or modelled climate inputs (Yuan et al., 2015; Zheng et al., 2013). For example, in ensemble streamflow prediction (ESP), hydrological models are driven by each historical season's precipitation and temperature series to produce an ensemble of flow forecasts, with this ensemble providing a distribution of plausible flows for the forecast period (Wood and Lettenmaier, 2008). ESP forecasts can also be used as an input predictor to statistical techniques (Robertson et al., 2013). Dynamical (i.e. climate model) forecasts of precipitation and temperature are often not sufficiently skilful in this region. For example, Kim et al. (2012) assessed retrospective seasonal forecasts of the Asian summer monsoon from ECMWF System 4 (Molteni et al., 2011) and NCEP CFSv2 (Saha et al., 2014), finding low skill for precipitation prediction and poor simulations of the Indian summer monsoon circulation. Cash et al. (2017) assessed monthly North American Multi-Model Ensemble (Kirtman et al., 2013) hindcasts initialised on 1 May for May to November for South Asia, including the mountainous areas of Afghanistan and Pakistan. They concluded that the multi-model ensemble-mean temperature and precipitation forecasts, while generally exceeding the skill of any individual model, provided little benefit over climatology. Given these findings, we have not investigated the use of dynamical forecasts in this assessment. Alternatively, statistically based forecast methods using robust relationships between climate drivers, antecedent catchment conditions and resultant streamflow can be valuable research and management tools when properly implemented (Plummer et al., 2009; Schepen et al., 2016). Thus in this assessment, we assess three statistically based forecast options for their practical feasibility in developing seasonal streamflow forecasting models for the study region:

  1. a statistical approach using the Bayesian joint probability (BJP) model with predictors accounting for antecedent basin conditions and climate drivers;

  2. an ESP approach using the snowmelt runoff model (SRM);

  3. a hybrid approach using option (1) with an additional predictor – the mean ESP forecasts from (2).

The study is reported as follows: Sect. 2 outlines the study area, details of the case studies and data used and climate influences; Sect. 3 presents the BJP statistical approach, the SRM–ESP approach and the verification metrics used to assess forecast skill, bias, reliability and robustness; Sect. 4 presents the results of the BJP and SRM skill scores and performance diagnostics. Section 5 discusses the performance of the forecast approaches and Sect. 6 concludes with the main findings and recommendations.

Figure 1Map of UIB showing sub-basins and location of major dams.


Figure 2Annual cycle of mean precipitation, PET and inflow for (a) Jhelum at Mangla and (b) Indus at Tarbela (note different volume scales).


2 Case study and data

2.1 Upper Indus Basin

Pakistan's water supply, crucial for its extensive irrigated agriculture industry, hydropower generation and industrial and municipal water supply, is predominantly sourced from Indus river flow, with groundwater a secondary although important contributor to most demands (with the exception of hydropower). The glaciated and snow covered sub-basins of the UIB, encompassing glaciated headwater catchments within the northern Hindu Kush, Karakoram and western Himalayan mountain ranges, dominate water generation within the Indus Basin (Alford et al., 2014). The UIB's tributaries include the Indus at Kharmong: Shigar, Shyok and Astore in the Karakoram Himalaya; the Jhelum, Chenab, Ravi and Sutlej in the western Himalaya; the Hunza, Gilgit, Kabul, Swat and Chitral in the Hindu Kush mountains (Fig. 1). These basins can be classified as having a flow regime that is either glacier-melt dominated (Hunza, Shigar and Shyok) or snowmelt dominated (Jhelum, Kabul, Gilgit, Astore and Swat; Hasson et al., 2014). The predominant source of flow in the UIB is snowmelt, with glacier melt a secondary source, with 80 % of flow occurring during the June–September summer period. Interannual flow variability is thus controlled by two processes, snow accumulation as determined by winter precipitation and temperature and meltwater generation as determined by summer temperatures. Hence snowmelt-generated flow is a function of winter precipitation and temperature and also summer temperature, whereas glacier melt is primarily a function of summer temperature, although glacier melt is also influenced by snow cover (Charles, 2016).

Inflows to two major reservoirs, Tarbela Dam on the Upper Indus and Mangla Dam on the Jhelum River, a major tributary to the Indus system, are investigated (Fig. 1). Daily inflow data from 1975 to 2015 were obtained from Pakistan Water and Power Development Authority (WAPDA). Figure 2 presents the seasonal hydroclimatic cycle for these two basins, showing double-peaked (winter and summer) precipitation with inflow peaking in May for Mangla and July for Tarbela. The Tarbela Dam on the main stem of the Indus is one of the largest individual storages in the UIB, crucial for hydropower generation and irrigation supply (Ahmed and Sanchez, 2011). Annual inflows to Tarbela constitute 70 % melt water, of which snowmelt contributes 44 % and glacial melt contribute 26 % (Mukhopadhyay and Khan, 2015). The Mangla Dam on the Jhelum River is (since enlargement) a similarly sized storage as Tarbela and one of the most important resources in Pakistan for electricity generation and water supply for irrigation (Mahmood et al., 2015). For the Jhelum, the area upstream of Mangla is reported as 33 500 km2 with an elevation range from 300 to 6285 m and mean of nearly 2400 m, the relatively low altitude ensures that there is only 0.7 % coverage by glaciers or perennial snow according to GLIMS glacier database as cited by Bogacki and Ismail (2016). In contrast, the Indus upstream of Tarbela is over 5 times larger (173 345 km2) with higher elevation (to 8238 m, as reported by Immerzeel et al., 2009) and 11.5 % is covered by perennial glaciers (Ismail and Bogacki, 2018), such that (as noted above and evident in Fig. 2, showing inflow exceeding precipitation) glacier ice melt makes a significant contribution to annual flow.

2.2 Climate influences

Useful climate indices should be related to the weather prior to the forecast season, providing an indication of snow accumulation, and also to the weather within the forecast season, influencing temperature and hence snow and glacier melt rates. A literature review identified indices related to the North Atlantic Oscillation (NAO) and El Niño Southern Oscillation (ENSO) as the most likely to provide skill for the UIB (Charles, 2016). These both influence the direction of prevailing winds bringing moisture into the region and thus determine precipitation and temperature conditions influencing the depth and areal extent of snow accumulation in the winter and early spring preceding the Kharif (April–September) high-flow season.

The NAO is a measure of the strength of the pressure gradient between the subtropics and polar regions in the North Atlantic, representing a dominant source of variability in circulation and winds influencing the region (Hurrell, 1995; Bierkens and van Beek, 2009). It has a direct influence on the interannual variability of the westerly winds (westerly disturbances) and their water content traversing Europe, the Mediterranean and the Middle East region into the mountains of the UIB (Yadav et al., 2009a; Syed et al., 2010; Filippi et al., 2014). Indices of the NAO have been linked with UIB weather station winter precipitation (Archer and Fowler, 2004; Afzal et al., 2013; Filippi et al., 2014), western Indus Basin's winter snow cover and weather station precipitation (Hasson et al., 2014), Pakistan weather station temperature (del Río et al., 2013) and winter precipitation in north-west India (Kar and Rana, 2014).

ENSO is a dominant pattern of multi-year variability driven by ocean–atmosphere interactions in the tropical Pacific (Wolter and Timlin, 2011), influencing climate globally including the variability of both western disturbances and monsoon processes experienced by the region. The commonly used SOI (Southern Oscillation Index) has been linked with winter Hindu Kush Himalayan region precipitation (Afzal et al., 2013), Indian Summer Monsoon precipitation (Ashok et al., 2004; Ashok and Saji, 2007), central south-west Asian winter precipitation (Syed et al., 2006), Pakistan weather station temperature (del Río et al., 2013) and north-west India winter precipitation (Kar and Rana, 2014). Stronger links have been reported between ENSO, western disturbances and interannual winter precipitation variability in recent decades (Yadav et al., 2009a, b).

3 Methods

3.1 BJP forecasting models

The statistical seasonal forecasting model used is the Bayesian joint probability (BJP) approach of Wang et al. (2009). The BJP has state-of-the-art capabilities in developing seasonal forecast models that optimally utilise information available on antecedent catchment conditions, large-scale climate forcing (through climate indices) and flow forecast scenarios from hydrological models (Robertson et al., 2013; Robertson and Wang, 2012; Schepen et al., 2012; Wang and Robertson, 2011). The BJP models simulate predictor–predictand relationships using conditional multivariate normal distributions, with predictor and predictand data transformed to normal using either a log-sinh (Wang et al., 2012b) or Yeo–Johnson (Yeo and Johnson, 2000) transformation. BJP parameters are inferred using Markov chain Monte Carlo (MCMC) methods to account for parameter uncertainty, which can be due to factors such as short data records. Probabilistic (ensemble) forecasts are produced by generating samples from the estimated conditional multivariate normal distributions. When predictor–predictand relationships are weak, the BJP produces reliable forecasts that approximate climatology. The full technical details of the BJP modelling approach are presented in Wang et al. (2009) and Wang and Robertson (2011).

3.2 SRM forecasting model

The Snowmelt Runoff Model (SRM) of Martinec et al. (2008) has been used in several studies in the basin (Butt and Bilal, 2011; Romshoo et al., 2015; Tahir et al., 2011; Bogacki and Ismail, 2016; Ismail and Bogacki, 2018). WAPDA has procured a version of SRM implemented in Excel©, “ExcelSRM”, and for their 2012 case study for Jhelum inflow into the Mangla Reservoir (NESPAK et al., 2012), ExcelSRM was calibrated using data for 2003 to 2010 and subsequently validated against inflows for 2000 to 2002 and 2007 and 2011 (Bogacki and Ismail, 2016). In contrast to the probabilistic forecasts produced by the BJP, the SRM is a deterministic model and so produces a single forecast for a given set of inputs.

Given the inadequacy of seasonal meteorological forecasts for the region (Bogacki and Ismail, 2016; Cash et al., 2017; Ismail and Bogacki, 2018), an ESP approach is used to forecast a range of possible Kharif season inflows. That is, SRM is initiated with snow cover observed at the end of March, and then run to produce 6-month Kharif season scenarios using the P and T inputs from each year in the available historical record, together with the modified depletion curve approach (Rango and Martinec, 1982) to simulate snow cover depletion as a function of that scenario year's degree-days series. This approach results in an ensemble of simulated inflows and, as well as assessing the SRM–ESP forecasts themselves (from a research version of ExcelSRM we obtained in 2015; i.e. not the current version used operationally), this study has used the mean of the ensemble members for each year as an additional predictor series for input to the BJP (option 3 as outlined in the Introduction).

3.3 Verification

BJP forecast model performance is verified using leave-one-out cross-validated results (Wang and Robertson, 2011). That is, to avoid artificially inflating the skill, for each Kharif season the calibration and assessment does not use that season's data for BJP parameter estimation. The cross-validated BJP forecast performance was assessed for 1975–2015 (41 seasons), with the BJP models calibrated on a seasonal basis (i.e. 40 data points) using 1000 MCMC samples for each of the leave-one-out calibrations.

As noted, the BJP can also use hydrological model simulated flow as an additional predictor (Robertson et al., 2013) and in this case we have SRM simulations for the Mangla inflow for a subset of the investigation period (2001–2015). The SRM is an exception to the leave-one-out cross validation as it is calibrated using all data for the 2003 to 2010 period, with parameters manually tuned “… in order to keep parameters at smooth values and to maintain a reasonable trend in time” (NESPAK et al., 2012). The mean flow simulations obtained from driving each year's SRM with all years' available precipitation and temperature are therefore not independent forecasts, and so it is not surprising that for the 2003 to 2010 period the SRM forecasts are closer to the observed flows than the median cross-validated BJP estimates. When used as a predictor to the BJP, SRM forecasts are applied with leave-one-out cross validation i.e. for each year in the 2001–2015 period, simulations from all of the other 14 years are used to produce an ESP with the resulting mean used as a predictor for that year. Note the BJP is able to extract skill from biased dynamic hydrological model forecasts, as long as the hydrological model simulation bias is systematic and stationary (i.e. not random or with a trend).

Verification assesses the overall skill and the bias, reliability and robustness of the forecasts. This includes assessing whether the bias and reliability of the forecasts vary for different periods of the record (temporal stability) or for different event sizes, e.g. whether there is a limitation in forecasting high- or low-flow seasons. Skill scores, quantifying the skill of the forecasts, allow the direct comparison of the performance of forecasting models that use different sets of predictors. Two common skill scores used here are the root mean squared error (RMSE), which assesses the forecast median, and the continuous ranked probability score (CRPS), which assesses the reduction in error of the whole forecast probability distribution (Robertson and Wang, 2013). The skill scores are reported as percentage reductions in error scores of the forecasts relative to the observed historical (climatological) median, for RMSE, and relative to the full distribution of the observed historical (climatological) events, for CRPS. The “sharpness” of a probabilistic forecast distribution (i.e. a narrower peaked distribution rather than a wide, flat distribution) is also a characteristic relevant to forecast skill (Gneiting et al., 2007). Sharp forecasts with narrow forecast intervals reduce the range of possible outcomes that are anticipated, increasing their usefulness for decision makers (Li et al., 2016). This skill can be quantified, for example, as the percentage reduction in the inter-quartile range (IQR) between the forecast's distribution and the observed historical (climatological) distribution (Crochemore et al., 2017). RMSE, CRPS and IQR skill scores are interpreted as “How are the skill score categories defined?” from (last access: 27 June 2018):

  • 0 is considered to be a forecast with no skill (equivalent skill to predicting using historical averages or historical reference);

  • less than 5 is considered to be a forecast with very low skill;

  • 5–15 is considered to be a forecast with low skill;

  • 15–30 is considered to be a forecast with moderate skill;

  • greater than 30 is considered to be a forecast with high skill.

Reliability refers to the statistical similarity between the forecast probabilities and the relative frequencies of events in the observations, which can be verified using probability integral transforms (PITs). The PIT represents the non-exceedance probability of observed streamflow obtained from the cumulative distribution function of the ensemble forecast. If the forecast ensemble spread is appropriate and free of bias then observations will be contained within the forecast ensemble spread, with reliable forecasts having PIT values that follow a uniform distribution between 0 and 1 (Laio and Tamea, 2007). Thus PIT plots are an efficient diagnostic to visually evaluate whether the forecast probability distributions are too wide or too narrow or are biased (under or over estimating) in their prediction of the observed distribution (Wang and Robertson, 2011). As outlined by Thyer et al. (2009), PIT plot points falling on the 1 : 1 line indicate that the predicted distribution is a perfect match to the observed; observed PIT values of 0.0 or 1.0 indicate the corresponding observed data falls outside the predicted range, hence the predictive uncertainty is significantly underestimated; PIT values clustered around the mid-range (i.e. a low slope in the 0.4–0.6 uniform variate range) indicate the predictive uncertainty is overestimated; PIT values clustered around the tails (i.e. a high slope in the 0.4–0.6 uniform variate range) indicate the predictive uncertainty is underestimated; if PIT values at the theoretical median are higher than those of the uniform variate the predictions have an underprediction bias; if they are lower than the uniform variate, then the predictions have an overprediction bias.

Table 1Cross-validated skill scores of BJP forecasts for Kharif season for Jhelum at Mangla. Bootstrap 10th to 90th percentile resampling ranges are shown in brackets.

a 1976–2015; b 1975–2015; c 2001–2015 SRM–ESP mean, with L1OCV; d 2001–2015 SRM–ESP with L1OCV.

Download Print Version | Download XLSX

4 Results

4.1 Skill scores

BJP models were trialled with combinations of predictors accounting for antecedent flow (flow immediately preceding the forecast season, i.e. March flow for the Kharif forecast) and NAO- or ENSO-based climate indices identified from the literature, as introduced in Sect. 2.2 (Charles, 2016). MODIS (Hall et al., 2010) snow-cover area and GLDAS-2.1 (Rodell et al., 2004) snow-water equivalent (SWE), additional measures of antecedent conditions, were also assessed as potential predictors. A significant limitation is the shorter record lengths for MODIS and GLDAS, as available data for both start in 2000. This is of particular concern for the BJP's leave-one-out cross-validation, as using short records to identify dynamical mechanisms is susceptible to spurious skill. Correlation analysis, cognisant of the short 2000–2015 period, found that these snow products have a similar or lower correlation with Kharif flow (QKharif) compared to March flow (QMarch), and are relatively highly correlated with QMarch. Thus the limitation of short record length, lack of higher correlation with flow than that of the QMarch predictor, and relatively high cross-correlation with QMarch, leads us to conclude that they would not be expected to provide additional skill as a predictor for the BJP.

Table 2Cross-validated skill scores of BJP forecasts for Kharif season for Indus at Tarbela. Bootstrapped 10th to 90th percentile resampling ranges are shown in brackets.

a 1976–2015; b 1975–2015; c no SRM for Indus; d denotes lag to calendar year before flow season.

Download Print Version | Download XLSX

SRM–ESP scenario-mean forecasts were an additional predictor trialled for Jhelum at Mangla. Higher skill was generally obtained for predictor combinations using a flow predictor (March flow) together with either the Multivariate ENSO index (MEI;, last access: 27 June 2018; Wolter and Timlin, 1998) or the Southern Oscillation signal index (SSI;, last access: 27 June 2018; Trenberth, 1984) as a climate predictor. The seasons of the trialled climate predictors (i.e. their time-lag preceding the forecast season) were selected based on their highest linear correlations with flow (not shown). Table 1 presents cross-validated BJP forecast skill scores using the trialled combinations of antecedent flow and climate predictors for the Kharif season for Jhelum at Mangla, together with bootstrapped 10th to 90th percentile ranges to assess model uncertainty. These ranges were obtained by resampling 1000 random sequences of years of the same length as the observed record, i.e. with replacement, and calculating skill scores for each sample. Combinations including the SRM forecasts (ESP mean) as a predictor are included; however, because the SRM results are only available for the 15-year period 2001–2015, they are only providing skill during the 2001–2015 period when included as a BJP predictor for the full 41-year period (1975–2015). These results show

  • The antecedent predictor (March flow, QMarch) provides greater skill than any of the individual climate predictors used.

  • Two-predictor models using QMarch and the SRMKharif predictor give poorer skill scores compared to using QMarch alone.

  • Two-predictor models using QMarch and one climate predictor slightly improve (in most cases) the skill scores compared to using QMarch alone.

  • Given the large uncertainty in skill scores, we do not aim to select a “best” model. However, as there are many models with positive skill (i.e. better than climatology), using skilful models is plausible. Ideas on how to do this are discussed in Sect. 5.

Addition of the SRMKharif predictor to the two-predictor models using QMarch and one climate predictor does not improve skill scores. Table 2 presents the skill scores for the Kharif season forecasts for the Indus at Tarbela, also using the antecedent flow and climate predictors but without SRM forecasts in this case (as SRM was not available for Indus at Tarbela for this study). Similarly to the results for the Jhelum at Mangla, for the Indus at Tarbela BJP forecasts:

  • The antecedent predictor (March flow, QMarch) provides greater skill than any of the climate predictors used.

  • On the whole, a single climate predictor produces low skill compared to that obtained using QMarch, with a notable exception that the MEIMayJun (i.e. the year before) predictor produces skill scores comparable to those obtained using QMarch. The selection of MEIMayJun as a predictor is discussed further in Sect. 5.

  • Two-predictor models using QMarch and one climate predictor improve the skill scores compared to using QMarch alone.

  • Given there is large uncertainty in skill scores we do not aim to select a “best” model. However, as there are many models with positive skill (i.e. better than climatology) then using skilful models is plausible. Ideas on how to do this are discussed in Sect. 5.

In addition to calibration for the full Kharif season, BJP calibrations were also undertaken for the early Kharif (April–June) and the late Kharif (July–September) using the relevant flow and ENSO-based predictors (e.g. for late Kharif the June flow was used as an antecedent predictor). A comparison of the resulting skill scores are shown in Fig. 3 for the BJP models that gave the highest skill gain relative to climatology for the Kharif, early Kharif and late Kharif periods. It is interesting to contrast the performance for the two locations, with Kharif and late Kharif giving similar results across the two locations, whereas for early Kharif a marked difference is seen, with high skill for Jhelum at Mangla contrasting the low skill for Indus at Tarbela. Determining the underlying physical reasons for this contrast would require further investigation, with possible causes discussed in Sect. 5.

4.2 Performance diagnostics

Here we assess the cross-validated performance of forecasts from BJP models using an antecedent and climate predictor combination (option 1) selected on the basis of skill scores and, for Mangla, we compare these with the SRM–ESP forecasts (option 2). We do not compare results for the BJP models using the SRM–ESP mean as an additional predictor (option 3), as the addition of this predictor added little or no skill to BJP forecasts (Table 1).

Figure 3BJP cross-validated skill scores, % skill gain relative to climatology, for CRPS, RMSE and IQR. Less than 5 is considered to be a forecast with very low skill. Between 5 and 15 is considered low skill. Between 15 and 30 is considered moderate skill, and higher than 30 is considered to be a forecast with high skill.


We use PIT plots for verification of the reliability and robustness of the forecast probability distributions, to assess whether there are biases in the forecasts, or whether the forecast probability distributions are too wide or too narrow (Laio and Tamea, 2007). For reliable forecasts, the PIT values should follow a uniform distribution and hence follow a 1 : 1 line when plotted against a standard uniform variate. For the BJP forecasts for the full 1975–2015 period, both Jhelum at Mangla (Fig. 4a) and Indus at Tarbela (Fig. 5a) are reliable (i.e. forecast probability distributions are unbiased and of appropriate spread), evidenced by the forecast's PIT values plotting close to the 1 : 1 lines and within the Kolmogorov 5 % significance band. Comparison of Jhelum at Mangla BJP and SRM–ESP forecasts, for the shorter 2001–2015 period for which SRM results are available, show a contrast between reliable BJP forecasts (Fig. 6a) and biased SRM forecasts, including five values of 0 or 1 indicating that the SRM forecast distribution is too narrow (Fig. 7a).

A feature of robust forecasts is their stability across the full period of record and range of flow magnitudes. Figures 4b and 5b show a uniform spread of PIT values and hence stability across the full period of record, for Jhelum at Mangla and Indus at Tarbela BJP forecasts respectively. Similarly, forecast stability across the range of flow magnitudes is verified by the uniformity of PIT values against forecast median, as shown in Figs. 4c and 5c for Jhelum at Mangla and Indus at Tarbela respectively. For the Jhelum at Mangla BJP and SRM–ESP comparison, stability across time and flow magnitude are harder to assess given the short 15-year sample size. Figure 6b and c show reasonable stability of the BJP forecasts, although there is a trend over time for this part of the record. The equivalent SRM plots (Fig. 7b and c) are not as robust, with (as noted previously) five values at 0 or 1 (from the chronological plot: 2001 is at 0, and 2010, 2012, 2014 and 2015 are at 1).

Robustness is also assessed by plotting forecast quantile ranges and observed flows against the forecast median (Figs. 4d and 5d) and chronologically (Figs. 4e and 5e, for Jhelum at Mangla and Indus at Tarbela respectively). These show that BJP forecasts reasonably account for the range of observed variability for both locations. The relatively less robust SRM–ESP forecasts are shown in Fig. 7d and e, again highlighting the overly narrow forecast distribution range for this version of SRM.

Overall, the performance shown in these figures highlight the reliability and robustness of the BJP forecasts for the Kharif season for Jhelum at Mangla and Indus at Tarbela. For the Jhelum at Mangla BJP and SRM–ESP forecasts for the shorter 2001–2015 period, the results contrast with the reliable and robust BJP against some limitations from the SRM, which will be discussed further in the next section.

5 Discussion

The SRM forecasts are examples of the commonly applied ESP approach (Shi et al., 2008; Shukla and Lettenmaier, 2011; Wood et al., 2005). As such, Wood and Schaake (2008) note

One strength of the ESP approach is that it accounts for uncertainty in future climate, which in some seasons is the major component of forecast uncertainty, by assuming that historical climate variability is a good estimate of current climate uncertainty. A weakness of the approach, however, is that when the uncertainty of the current (“initial”) hydrologic state is a significant component of the overall forecast uncertainty …, the deterministic estimate of the forecast ensemble's initial hydrologic state leads to an overconfident forecast – that is, one having a spread that is narrower than the total forecast uncertainties warrant.

This can be seen in the poor verification performance of SRM–ESP shown in Sect. 4.2, with the SRM–ESP from the 15-year sample unable to account for the observed range of flows, i.e. the ESP range is too narrow, even though in terms of overall RMSE and mean absolute percentage error (MAPE) statistics the SRM–ESP mean and BJP median are comparable (18.7 % RMSE and 14.3 % MAPE for BJP median; 18.4 % RMSE and 12.7 % MAPE for SRM–ESP mean). Given that the observed climate of each individual year would, in most cases, be within the range of the ensemble of climate inputs used to produce the ESP, this indicates SRM formulation could be too strongly reliant on antecedent conditions at the beginning of the forecast season. An additional source of bias, as evidenced by the years of poorest performance being outside the years used for parameter estimation, could be over-fitting with the model parameters tied too closely to the range of observed predictor–predictand relationships in the 2003–2010 calibration period.

Figure 4BJP cross-validated forecasts for Jhelum at Mangla Kharif for 1975–2015; (a) PIT uniform probability plot (1 : 1 black line, theoretical uniform distribution; grey lines, Kolmogorov 5 % significance bands; blue points, PIT values of forecast streamflow); (b) chronological PIT plot; (c) median PIT plot; (d) forecast quantiles and observed plotted according to forecast median (1 : 1 line, forecast median; dark vertical line, forecast [0.25, 0.75] quantile range; light and dark vertical line, forecast [0.10, 0.90] quantile range; dots, observed inflow); (e) chronological forecast quantile range and observations (dark blue, forecast [0.25, 0.75]; light and dark blue, forecast [0.10, 0.90]; crosses, forecast median; dots, observed).


Figure 5As in Fig. 4 for BJP cross-validated forecasts for Indus at Tarbela Kharif for 1975–2015.


Figure 6As in Fig. 4 for BJP cross-validated forecasts for Jhelum at Mangla Kharif for SRM period of 2001–2015, except for (e) see 2001–2015 period of Fig. 4e.


Figure 7As in Fig. 4 for SRM ESP forecasts for Jhelum at Mangla Kharif for SRM period of 2001–2015.


The poorer BJP performance for Indus at Tarbela, as seen in the skill scores relative to Jhelum at Mangla (Fig. 3), could be related to the differences in flow generation mechanisms. As the predictors are the source of skill in the statistical BJP approach, examination of correlations between the predictors and flow for the individual months within the Kharif season is insightful. Table 3 shows the (intuitively) expected pattern for Jhelum at Mangla of QMarch having a maximum correlation with April flow (0.84) and then maintaining a relatively high correlation until July (0.62) before dropping off for August (0.12). A different process appears to be influencing Indus at Tarbela, as the initial highest QMarch correlation with April flow (0.66) drops immediately to 0.18 for May before oscillating between 0.25 (August) and 0.55 (September) for subsequent Kharif months. Similarly, the climate predictor's correlations with the individual month's flows show more of a gradual reduction for Jhelum at Mangla (high for the first 4 months), whereas for Indus at Tarbela again an oscillatory relationship is seen. These higher correlations in the late Kharif (relative to the early Kharif) for Indus at Tarbela could be related to the correspondingly higher relative skill scores shown in Fig. 3 for late Kharif, corresponding to late-season glacier melt processes that are a significant component of the inflow to Tarbela but not Mangla (Mukhopadhyay and Khan, 2015). Future research could investigate whether dynamical seasonal forecasts of temperature have skill of relevance to forecasting glacier melt; however, as noted above such skill has not been determined to date (e.g. Cash et al., 2017), and is beyond the scope of this assessment given our focus of developing practical and easily implementable forecast tools using readily available inputs.

Table 3Correlation between the flow of the individual months of the Kharif season and the predictors used by BJP (values in bold significant at p< 0.05).

* year before.

Download Print Version | Download XLSX

It is also interesting to reflect on the relative performance of the NAO climate predictor, which does not provide any skill for inflow to Mangla (Table 1) but offers comparable skill to several of the ENSO indices trialled for Tarbela (Table 2). This indicates NAO may have some skill with regard to late season glacier melt. Overall, these results concur with investigations showing a stronger relationship between ENSO and precipitation and weaker relationship between NAO and precipitation in recent decades (Yadav et al., 2009a, b) resulting in the prevalence of ENSO as the better predictor of winter snowpack magnitude.

For Mangla, the predictor combination that gave the best Kharif season cross-validated skill scores included an ENSO-based predictor (SSIMarch) immediately before the season (Table 1), which makes sense intuitively as it represents a climate driver of both the snow accumulation before the Kharif season and precipitation conditions during. In contrast, for Tarbela a much earlier ENSO-based predictor (MEIMayJun, i.e. May–June the year before) provides higher skill scores than the equivalent predictor immediately before the season (MEIFebMar; Table 2). To try to understand the dynamical mechanism by which MEIMayJun is providing skill in forecasting QKharif, we compared MEI correlations with GLDAS SWEMarch, QMarch and QKharif. Results were inconclusive, and perhaps impeded by the short record lengths given SWE is only available from 2000, as while MEIMayJun has a higher correlation with QKharif than MEIFebMar (0.76 versus 0.63 respectively), it has a slightly lower correlation with SWEMarch (0.48 versus 0.52 respectively). Hence MEIMayJun does not appear to be a long-lead predictor of snow accumulation, and so the differences in skill scores may be due to spurious correlations. Therefore we recommend comparing and assessing both this model and the QMarch and MEIFebMar model for future events. More generally, the skill score uncertainty ranges presented in Tables 1 and 2 highlight that no “best” forecast model can be selected for either basin. Attempting to select a best model would ignore model uncertainty and thus not make best use of forecast skill across the range of models trialled. To address this, probabilistic forecasts from multiple BJP models can be combined using Bayesian model averaging to produce combined forecasts with higher skill than that obtainable from any individual model (Wang et al., 2012a). Thus trialling a Bayesian model averaging (BMA) approach is recommended, although it is beyond the scope of this current work.

6 Conclusion

This study has assessed the performance and practical feasibility of three options for producing Kharif (April–September) seasonal streamflow forecasts for the Jhelum River inflows to the Mangla Dam in the UIB of Pakistan: option 1, the BJP statistical forecasting technique; option 2, the SRM physically based model run in ESP mode; option 3, a hybrid of option 1 with the mean ESP forecasts from option 2 used as an additional predictor for input to the BJP. The option 1 BJP forecast model used antecedent catchment and climatic predictors, with the predictors selected based on BJP skill score performance. The selected predictors represent hydrological conditions immediately preceding the forecast season (i.e. flow of the preceding month – March in this case) and ENSO-based climate indices related to drivers of winter snow accumulation. For an additional comparison, the option 1 BJP approach was also undertaken for the Indus River inflows to the Tarbela Dam.

Overall findings were as follows:

  • The best performing BJP models for Tarbela and Mangla inflows are consistent in that both used an antecedent flow predictor and a climate predictor representing ENSO.

  • For Tarbela the QMarch and MEIMayJun model gave the best skill; however, because we could not determine the dynamical mechanism(s) by which the relatively long lag between MEIMayJun influences snowpack accumulation and flow, we cannot rule out the possibility that the skill is due to spurious correlation. Therefore we recommend comparing and assessing both this model and the QMarch and MEIFebMar model for future events and, more generally, that BMA be trialled in future research to combine the skill of multiple BJP models as, for example, undertaken in Australia in Pokhrel et al. (2013).

  • There are pragmatic benefits to selecting a BJP model using only antecedent and climate predictors, rather than including SRM mean ESP as an additional predictor even in cases when SRM does provide skill, given that flow and climate predictors are readily available and thus BJP forecasts are easily and quickly produced. The SRM, being a deterministic model, is a much more technically involved and data-intensive approach to forecast generation (Bogacki and Ismail, 2016; Ismail and Bogacki, 2018).

  • Cross-validated performance of the BJP seasonal forecasts for the 1975 to 2015 Kharif seasons, as shown in the diagnostic and verification statistics presented, highlight that the BJP produces forecasts that are statistically unbiased, robust and reliable. In contrast, the SRM–ESP forecasts show bias particularly for the most recent years outside the SRM calibration period, potentially indicating limitations with the SRM due to lack of cross-validated calibration and resultant over-fitting. Thus SRM–ESP forecasts are overly confident, underestimating the full uncertainty that is captured by the BJP approach.

  • High skill was obtained for BJP forecasts of early Kharif flow for Jhelum at Mangla. Moderate skill was obtained for the full and late Kharif season forecasts for both Jhelum at Mangla and Indus at Tarbela. Lower skill was seen for early Kharif for Indus at Tarbela.

In future research, BJP forecast models could readily be developed and assessed for other tributaries, e.g. the Chenab and Kabul, subject to availability of flow data. This would allow an overall assessment of UIB flow forecasting for the major contributing basins. The present method used by the Indus River System Authority (IRSA) to forecast UIB Kharif streamflow is based on historical analogues. The IRSA use their database of the previous 60 years of flow to select years where the historic March flows are within 5 % of the current March flow and use the corresponding historical Kharif flows (within 5 %) for their forecast scenario. The selection of the historical scenario is also informed by forecasts from the Pakistan Meteorological Department, forecasts provided by WAPDA (e.g. SRM forecasts) and present snow conditions in the catchment. The forecasts are continuously revised as the season progresses.

Sufficiently skilful BJP forecasts could also inform scenario selection, providing for the first time a probabilistic approach to forecasts in contrast to a single forecast as currently used. However probabilistic forecasts (such as the BJP) can be misinterpreted if the water management professionals using them to inform decisions are unfamiliar with them (Pagano et al., 2002; Ramos et al., 2013; Rayner et al., 2005; Whateley et al., 2015). Hence the successful transfer of BJP forecast tools to operational use within Pakistan would require guidance for building BJP models and generating forecasts, test cases with example results, face-to-face training and on-going support.

Data availability

Data sets, with the exception of inflow data, are available upon request by contacting the correspondence author.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Sub-seasonal to seasonal hydrological forecasting”. It is not associated with a conference.


The Pakistan Mission of the Australian Department of Foreign Affairs and Trade provided funds for this study through the Sustainable Development Investment Portfolio Indus project. CSIRO provided additional support through its Land and Water portfolio. We thank the Pakistan Meteorological Department (PMD) for sharing datasets and the Pakistan Water and Power Development Authority Glacier Monitoring Research Centre (WAPDA-GMRC) for sharing datasets and the SRM used in this collaborative work. Thanks also to Luis Neumann and Tony Zhao, CSIRO Land and Water, for help and advice. We express our gratitude to the two anonymous reviewers and the editor, Andy Wood, for their input, which helped improve the paper. We also acknowledge helpful comments received from Muhammad Fraz Ismail regarding the SRM.

Edited by: Andy Wood
Reviewed by: two anonymous referees


Afzal, M., Haroon, M. A., Rana, A. S., and Imran, A.: Influence of North Atlantic Oscillations and Southern Oscillations on winter precipitation of northern Pakistan, Pakistan J. Meteorol., 9, 1–8, 2013. 

Ahmed, K. and Sanchez, M.: A study of the factors and processes involved in the sedimentation of Tarbela reservoir, Pakistan, Environ. Earth Sci., 62, 927–933,, 2011. 

Alford, D., Archer, D., Bookhagen, B., Grabs, W., Halvorson, S., Hewitt, K., Immerzeel, W., Kamp, U., and Krumwiede, B.: Monitoring of glaciers, climate and runoff in the Hindu Kush-Himalaya mountains, South Asia Water Initiative Report No. 67668-SAS, World Bank, Washington, D.C., 2014. 

Anghileri, D., Voisin, N., Castelletti, A., Pianosi, F., Nijssen, B., and Lettenmaier, D. P.: Value of long-term streamflow forecasts to reservoir operations for water supply in snow-dominated river catchments, Water Resour. Res., 52, 4209–4225,, 2016. 

Archer, D. R. and Fowler, H. J.: Spatial and temporal variations in precipitation in the Upper Indus Basin, global teleconnections and hydrological implications, Hydrol. Earth Syst. Sci., 8, 47–61,, 2004. 

Ashok, K. and Saji, N. H.: On the impacts of ENSO and Indian Ocean dipole events on sub-regional Indian summer monsoon rainfall, Nat. Hazards, 42, 273–285,, 2007. 

Ashok, K., Guan, Z., Saji, N. H., and Yamagata, T.: Individual and combined influences of ENSO and the Indian Ocean Dipole on the Indian Summer Monsoon, J. Climate, 17, 3141–3155,<3141:IACIOE>2.0.CO;2, 2004. 

Asian Development Bank: Asian water development outlook 2016: Strengthening water security in Asia and the Pacific, Asian Development Bank, Mandaluyong City, Philippines, 2016. 

Bennett, J. C., Wang, Q. J., Li, M., Robertson, D. E., and Schepen, A.: Reliable long-range ensemble streamflow forecasts: Combining calibrated climate forecasts with a conceptual runoff model and a staged error model, Water Resour. Res., 52, 8238–8259,, 2016. 

Bierkens, M. F. P. and van Beek, L. P. H.: Seasonal Predictability of European Discharge: NAO and Hydrological Response Time, J. Hydrometeorol., 10, 953–968,, 2009. 

Bogacki, W. and Ismail, M. F.: Seasonal forecast of Kharif flows from Upper Jhelum catchment, Proc. IAHS, 374, 137–142,, 2016. 

Butt, M. J. and Bilal, M.: Application of snowmelt runoff model for water resource management, Hydrol. Process., 25, 3735–3747,, 2011. 

Cash, B. A., Manganello, J. V., and Kinter, J. L.: Evaluation of NMME temperature and precipitation bias and forecast skill for South Asia, Clim. Dynam.,, in press, 2017. 

Charles, S. P.: Hydroclimate of the Indus – synthesis of the literature relevant to Indus basin hydroclimate processes, trends, seasonal forecasting and climate change, CSIRO Sustainable Development Investment Portfolio project, CSIRO Land and Water, Australia, 48 pp., 2016. 

Chiew, F. H. S., Zhou, S. L., and McMahon, T. A.: Use of seasonal streamflow forecasts in water resources management, J. Hydrol., 270, 135–144,, 2003. 

Crochemore, L., Ramos, M. H., Pappenberger, F., and Perrin, C.: Seasonal streamflow forecasting by conditioning climatology with precipitation indices, Hydrol. Earth Syst. Sci., 21, 1573–1591,, 2017. 

del Río, S., Anjum Iqbal, M., Cano-Ortiz, A., Herrero, L., Hassan, A., and Penas, A.: Recent mean temperature trends in Pakistan and links with teleconnection patterns, Int. J. Climatol., 33, 277–290,, 2013. 

Doblas-Reyes, F. J., García-Serrano, J., Lienert, F., Biescas, A. P., and Rodrigues, L. R. L.: Seasonal climate predictability and forecasting: status and prospects, Wiley Interdisciplin. Rev.: Clim. Change, 4, 245–268,, 2013. 

Döll, P., Fiedler, K., and Zhang, J.: Global-scale analysis of river flow alterations due to water withdrawals and reservoirs, Hydrol. Earth Syst. Sci., 13, 2413–2432,, 2009. 

Filippi, L., Palazzi, E., von Hardenberg, J., and Provenzale, A.: Multidecadal Variations in the Relationship between the NAO and Winter Precipitation in the Hindu Kush–Karakoram, J. Climate, 27, 7890–7902,, 2014. 

Gneiting, T., Balabdaoui, F., and Raftery, A. E.: Probabilistic forecasts, calibration and sharpness, J. Roy. Stat. Soc. B, 69, 243–268,, 2007. 

Hall, D. K., Riggs, G. A., Foster, J. L., and Kumar, S. V.: Development and evaluation of a cloud-gap-filled MODIS daily snow-cover product, Remote Sens. Environ., 114, 496–503,, 2010. 

Hasson, S., Lucarini, V., Khan, M. R., Petitta, M., Bolch, T., and Gioli, G.: Early 21st century snow cover state over the western river basins of the Indus River system, Hydrol. Earth Syst. Sci., 18, 4077–4100,, 2014. 

Hurrell, J. W.: Decadal Trends in the North Atlantic Oscillation: Regional Temperatures and Precipitation, Science, 269, 676–679, 1995. 

Immerzeel, W. W., Droogers, P., de Jong, S. M., and Bierkens, M. F. P.: Large-scale monitoring of snow cover and runoff simulation in Himalayan river basins using remote sensing, Remote Sens. Environ., 113, 40–49,, 2009. 

Ismail, M. F. and Bogacki, W.: Scenario approach for the seasonal forecast of Kharif flows from the Upper Indus Basin, Hydrol. Earth Syst. Sci., 22, 1391–1409,, 2018. 

Kar, S. and Rana, S.: Interannual variability of winter precipitation over northwest India and adjoining region: impact of global forcings, Theor. Appl. Climatol., 116, 609–623,, 2014. 

Kim, H.-M., Webster, P., Curry, J., and Toma, V.: Asian summer monsoon prediction in ECMWF System 4 and NCEP CFSv2 retrospective seasonal forecasts, Clim. Dynam., 39, 2975–2991,, 2012. 

Kirby, M., Ahmad, M.-u.-D., Mainuddin, M., Khaliq, T., and Cheema, M. J. M.: Agricultural production, water use and food availability in Pakistan: Historical trends, and projections to 2050, Agr. Water Manage., 179, 34–46,, 2017. 

Kirtman, B. P., Min, D., Infanti, J. M., Kinter, J. L., Paolino, D. A., Zhang, Q., van den Dool, H., Saha, S., Mendez, M. P., Becker, E., Peng, P., Tripp, P., Huang, J., DeWitt, D. G., Tippett, M. K., Barnston, A. G., Li, S., Rosati, A., Schubert, S. D., Rienecker, M., Suarez, M., Li, Z. E., Marshak, J., Lim, Y.-K., Tribbia, J., Pegion, K., Merryfield, W. J., Denis, B., and Wood, E. F.: The North American Multimodel Ensemble: Phase-1 Seasonal-to-Interannual Prediction; P hase-2 toward Developing Intraseasonal Prediction, B. Am. Meteorol. Soc., 95, 585–601,, 2013. 

Koster, R. D., Mahanama, S. P. P., Livneh, B., Lettenmaier, D. P., and Reichle, R. H.: Skill in streamflow forecasts derived from large-scale estimates of soil moisture and snow, Nat. Geosci., 3, 613–616, 2010. 

Laio, F. and Tamea, S.: Verification tools for probabilistic forecasts of continuous hydrological variables, Hydrol. Earth Syst. Sci., 11, 1267–1277,, 2007. 

Li, H., Luo, L., Wood, E. F., and Schaake, J.: The role of initial conditions and forcing uncertainties in seasonal hydrologic forecasting, J. Geophys. Res.-Atmos., 114, D04114,, 2009. 

Li, M., Wang, Q. J., Bennett, J. C., and Robertson, D. E.: Error reduction and representation in stages (ERRIS) in hydrological modelling for ensemble streamflow forecasting, Hydrol. Earth Syst. Sci., 20, 3561–3579,, 2016. 

Mahmood, R., Babel, M. S., and Jia, S.: Assessment of temporal and spatial changes of future climate in the Jhelum river basin, Pakistan and India, Weather Clim. Extr., 10, 40–55,, 2015. 

Martinec, J., Rango, A., and Roberts, R.: Snowmelt Runoff Model (SRM) User's Manual, New Mexico State University, Las Cruces, New Mexico, 178 pp, (last access: 27 June 2018), 2008. 

Maurer, E. P. and Lettenmaier, D. P.: Predictability of seasonal runoff in the Mississippi River basin, J. Geophys. Res.-Atmos., 108, 8607,, 2003. 

Molteni, F., Stockdale, T. M., Balmaseda, A., Balsamo, G., Buizza, R., Ferranti, L., Magnusson, L., Mogensen, K., Palmer, T., and Vitart, F.: The new ECMWF seasonal forecast system (System 4), ECMWF Technical Memorandum 656, ECMWF – European Centre for Medium Range Weather Forecasts, Shinfield Park, Reading, Berkshire, England, 2011. 

Mukhopadhyay, B. and Khan, A.: A reevaluation of the snowmelt and glacial melt in river flows within Upper Indus Basin and its significance in a changing climate, J. Hydrol., 527, 119–132,, 2015. 

NESPAK, AHT, and DELTARES: Hydrological Flow Forecast Model for Mangla Catchment, Upgrading of Tools, Water Resources Database, Management Systems and Models Under Sub Component “B1” of WCAP, Final Report, Lahore, Pakistan, 2012. 

Pagano, T. C., Hartmann, H. C., and Sorooshian, S.: Factors affecting seasonal forecast use in Arizona water management: a case study of the 1997–98 El Niño, Clim. Res., 21, 259–269, 2002. 

Plummer, N., Tuteja, N., Wang, Q. J., Wang, E., Robertson, D., Zhou, S., Schepen, A., Alves, O., Timbal, B., and Puri, K.: A seasonal water availability prediction service: opportunities and challenges, in: 18th World IMACS/MODSIM Congress, 13–17 July 2009, Cairns, Australia, 80–94, 2009. 

Pokhrel, P., Wang, Q. J., and Robertson, D. E.: The value of model averaging and dynamical climate model predictions for improving statistical seasonal streamflow forecasts over Australia, Water Resour. Res., 49, 6671–6687,, 2013. 

Ramos, M. H., van Andel, S. J., and Pappenberger, F.: Do probabilistic forecasts lead to better decisions?, Hydrol. Earth Syst. Sci., 17, 2219–2232,, 2013. 

Rango, A. and Martinec, J.: Snow accumulation derived from modified depletion curves of snow coverage, Hydrological Aspects of Alpine and High Mountain Areas, in: Proceedings of the Exeter Symposium, July 1982, Exeter, 83–90, 1982. 

Rayner, S., Lach, D., and Ingram, H.: Weather Forecasts are for Wimps: Why Water Resource Managers Do Not Use Climate Forecasts, Climatic Change, 69, 197–227,, 2005. 

Robertson, D. E. and Wang, Q. J.: Seasonal Forecasts of Unregulated Inflows into the Murray River, Australia, Water Resour. Manage., 27, 2747–2769,, 2013. 

Robertson, D. E. and Wang, Q. J.: A Bayesian approach to predictor selection for seasonal streamflow forecasting, J. Hydrometeorol., 13, 155–171,, 2012. 

Robertson, D. E., Pokhrel, P., and Wang, Q. J.: Improving statistical forecasts of seasonal streamflows using hydrological model output, Hydrol. Earth Syst. Sci., 17, 579–593,, 2013. 

Rodell, M., Houser, P. R., Jambor, U., Gottschalck, J., Mitchell, K., Meng, C. J., Arsenault, K., Cosgrove, B., Radakovich, J., Bosilovich, M., Entin, J. K., Walker, J. P., Lohmann, D., and Toll, D.: The Global Land Data Assimilation System, B. Am. Meteorol. Soc., 85, 381–394,, 2004. 

Romshoo, S. A., Dar, R. A., Rashid, I., Marazi, A., Ali, N., and Zaz, S. N.: Implications of Shrinking Cryosphere Under Changing Climate on the Streamflows in the Lidder Catchment in the Upper Indus Basin, India, Arct. Antarct. Alp. Res., 47, 627–644,, 2015. 

Saha, S., Moorthi, S., Wu, X., Wang, J., Nadiga, S., Tripp, P., Behringer, D., Hou, Y.-T., Chuang, H.-Y., Iredell, M., Ek, M., Meng, J., Yang, R., Mendez, M. P., van den Dool, H., Zhang, Q., Wang, W., Chen, M., and Becker, E.: The NCEP Climate Forecast System Version 2, J. Climate, 27, 2185–2208,, 2014. 

Schepen, A., Wang, Q. J., and Robertson, D.: Evidence for Using Lagged Climate Indices to Forecast Australian Seasonal Rainfall, J. Climate, 25, 1230–1246,, 2012. 

Schepen, A., Zhao, T., Wang, Q. J., Zhou, S., and Feikema, P.: Optimising seasonal streamflow forecast lead time for operational decision making in Australia, Hydrol. Earth Syst. Sci., 20, 4117–4128,, 2016. 

Schlosser, C. A., Strzepek, K., Gao, X., Fant, C., Blanc, É., Paltsev, S., Jacoby, H., Reilly, J., and Gueneau, A.: The future of global water stress: An integrated assessment, Earth's Future, 2, 341–361,, 2014. 

Shi, X., Wood, A. W., and Lettenmaier, D. P.: How essential is hydrologic model calibration to seasonal streamflow forecasting?, J. Hydrometeorol., 9, 1350–1363,, 2008. 

Shukla, S. and Lettenmaier, D. P.: Seasonal hydrologic prediction in the United States: understanding the role of initial hydrologic conditions and seasonal climate forecast skill, Hydrol. Earth Syst. Sci., 15, 3529–3538,, 2011. 

Shukla, S., Sheffield, J., Wood, E. F., and Lettenmaier, D. P.: On the sources of global land surface hydrologic predictability, Hydrol. Earth Syst. Sci., 17, 2781–2796,, 2013. 

Syed, F. S., Giorgi, F., Pal, J. S., and King, M. P.: Effect of remote forcings on the winter precipitation of central southwest Asia part 1: observations, Theor. Appl. Climatol., 86, 147–160,, 2006. 

Syed, F. S., Giorgi, F., Pal, J. S., and Keay, K.: Regional climate model simulation of winter climate over Central–Southwest Asia, with emphasis on NAO and ENSO effects, Int. J. Climatol., 30, 220–235,, 2010. 

Tahir, A. A., Chevallier, P., Arnaud, Y., Neppel, L., and Ahmad, B.: Modeling snowmelt-runoff under climate scenarios in the Hunza River basin, Karakoram Range, Northern Pakistan, Journal of Hydrology, 409, 104-117,, 2011. 

Thyer, M., Renard, B., Kavetski, D., Kuczera, G., Franks, S. W., and Srikanthan, S.: Critical evaluation of parameter consistency and predictive uncertainty in hydrological modeling: A case study using Bayesian total error analysis, Water Resour. Res., 45, W00B14,, 2009. 

Trenberth, K. E.: Signal Versus Noise in the Southern Oscillation, Mon. Weather Rev., 112, 326–332,<0326:SVNITS>2.0.CO;2, 1984. 

van Dijk, A. I. J. M., Peña-Arancibia, J. L., Wood, E. F., Sheffield, J., and Beck, H. E.: Global analysis of seasonal streamflow predictability using an ensemble prediction system and observations from 6192 small catchments worldwide, Water Resour. Res., 49, 2729–2746,, 2013. 

Wada, Y., van Beek, L. P. H., and Bierkens, M. F. P.: Modelling global water stress of the recent past: on the relative importance of trends in water demand and climate variability, Hydrol. Earth Syst. Sci., 15, 3785–3808,, 2011. 

Wang, Q. J. and Robertson, D. E.: Multisite probabilistic forecasting of seasonal flows for streams with zero value occurrences, Water Resour. Res., 47, W02546,, 2011. 

Wang, Q. J., Robertson, D. E., and Chiew, F. H. S.: A Bayesian joint probability modeling approach for seasonal forecasting of streamflows at multiple sites, Water Resour. Res., 45, W05407,, 2009. 

Wang, Q. J., Schepen, A., and Robertson, D. E.: Merging Seasonal Rainfall Forecasts from Multiple Statistical Models through Bayesian Model Averaging, J. Climate, 25, 5524–5537,, 2012a. 

Wang, Q. J., Shrestha, D. L., Robertson, D. E., and Pokhrel, P.: A log-sinh transformation for data normalization and variance stabilization, Water Resour. Res., 48, W05514,, 2012b. 

Whateley, S., Palmer, R. N., and Brown, C.: Seasonal Hydroclimatic Forecasts as Innovations and the Challenges of Adoption by Water Managers, J. Water Resour. Pl. Manage., 141, 04014071,, 2015. 

Wolter, K., and Timlin, M. S.: Measuring the strength of ENSO events: How does 1997/98 rank?, Weather, 53, 315-324,, 1998. 

Wolter, K. and Timlin, M. S.: El Niño/Southern Oscillation behaviour since 1871 as diagnosed in an extended multivariate ENSO index (MEI.ext), Int. J. Climatol., 31, 1074–1087,, 2011. 

Wood, A. W. and Lettenmaier, D. P.: An ensemble approach for attribution of hydrologic prediction uncertainty, Geophys. Res. Lett., 35, L14401,, 2008. 

Wood, A. W. and Schaake, J. C.: Correcting Errors in Streamflow Forecast Ensemble Mean and Spread, J. Hydrometeorol., 9, 132–148,, 2008. 

Wood, A. W., Kumar, A., and Lettenmaier, D. P.: A retrospective assessment of National Centers for Environmental Prediction climate model-based ensemble hydrologic forecasting in the western United States, J. Geophys. Res.-Atmos., 110, D04105,, 2005. 

Wood, A. W., Hopson, T., Newman, A., Brekke, L., Arnold, J., and Clark, M.: Quantifying Streamflow Forecast Skill Elasticity to Initial Condition and Climate Prediction Skill, J. Hydrometeorol., 17, 651–668,, 2015. 

Yadav, R. K., Rupa Kumar, K., and Rajeevan, M.: Increasing influence of ENSO and decreasing influence of AO/NAO in the recent decades over northwest India winter precipitation, J. Geophys. Res.-Atmos., 114, D12112,, 2009a.  

Yadav, R. K., Yoo, J. H., Kucharski, F., and Abid, M. A.: Why is ENSO influencing northwest India winter precipitation in recent decades?, J. Climate, 23, 1979–1993,, 2009b. 

Yeo, I.-K. and Johnson, R. A.: A New Family of Power Transformations to Improve Normality or Symmetry, Biometrika, 87, 954–959, 2000. 

Yossef, N. C., Winsemius, H., Weerts, A., van Beek, R., and Bierkens, M. F. P.: Skill of a global seasonal streamflow forecasting system, relative roles of initial conditions and meteorological forcing, Water Resour. Res., 49, 4687–4699,, 2013. 

Yuan, X., Wood, E. F., and Ma, Z.: A review on climate-model-based seasonal hydrologic forecasting: physical understanding and system development, Wiley Interdisciplin. Rev.: Water, 2, 523–536,, 2015. 

Zheng, H. X., Wang, Q. J., Aynul, K., Shao, Q. X., Shin, D., and Tuteja, N.: Evaluation of Downscaled POAMA M24 for Monthly and 3-Monthly Streamflow Forecasts, in: 20th International Congress on Modelling and Simulation (Modsim2013), 1–6 December 2013, Adelaide, Australia, 2799–2805, 2013. 

Short summary
Predictions of irrigation-season water availability are important for water-limited Pakistan. We assess a Bayesian joint probability approach, using flow and climate indices as predictors, to produce streamflow forecasts for inflow to Pakistan's two largest dams. The approach produces skilful and reliable forecasts. As it is simple and quick to apply, it can be used to provide probabilistic seasonal streamflow forecasts that can inform Pakistan's water management.