Statistical forecast of seasonal discharge in Central Asia for water resources management: development of a generic linear modelling tool for operational use

The semi-arid regions of Central Asia crucially depend on the water resources supplied by the mountainous areas 15 of the Tien Shan, Pamir and Altai mountains. During the summer months the snow and glacier melt dominated river discharge originating in the mountains provides the main water resource available for agricultural production, but also for storage in reservoirs for energy generation during the winter months. Thus a reliable seasonal forecast of the water resources is crucial for a sustainable management and planning of water resources. In fact, seasonal forecasts are mandatory tasks of all national hydro-meteorological services in the region. In order to support the operational seasonal forecast procedures of hydro20 meteorological services, this study aims at the development of a generic tool for deriving statistical forecast models of seasonal river discharge. The generic model is kept as simple as possible in order to be driven by available meteorological and hydrological data, and be applicable for all catchments in the region. As snowmelt dominates summer runoff, the main meteorological predictors for the forecast models are monthly values of winter precipitation and temperature, satellite based snow cover data and antecedent discharge. This basic predictor set was further extended by multi-monthly means of the 25 individual predictors, as well as composites of the predictors. Forecast models are derived based on these predictors as linear combinations of up to 3 or 4 predictors. A user selectable number of best models is extracted automatically by the developed model fitting algorithm, which includes a test for robustness by a leave-one-out cross validation. Based on the cross validation the predictive uncertainty was quantified for every prediction model. Forecasts of the mean seasonal discharge of the period April to September are derived every month starting from January until June. The application of the model for several 30 catchments in Central Asia ranging from small to the largest rivers (240 km2 to 290,000 km2 catchment area) – for the period

2000-2015 provided skilful forecasts for most catchments already in January with adjusted R 2 values of the best model in the range of 0.3 -0.8. The skill of the prediction increased every following month, i.e. with reduced lead time, with adjusted R 2 values usually in the range 0.8 -0.9 for the best and 0.7 -0.8 for the ensemble mean in April just before the prediction period.
The later forecasts in May and June improve further due to the high predictive power of the discharge in the first 2 months of the snow melt period. The improved skill of the model ensemble with decreasing lead time resulted in very narrow predictive 5 uncertainty bands at the beginning of the snow melt period. In summary, the proposed generic automatic forecast model development tool provides robust predictions for seasonal water availability in Central Asia, which will be tested against the official forecasts in the upcoming years, with the vision of operational implementation.

Introduction
The Central Asian region encompassing the five countries Kazakhstan, Kyrgyzstan, Tajikistan, Turkmenistan and Uzbekistan as well as northern parts of Afghanistan and north-western regions of China is characterized by the presence of two major mountain systems Tien Shan and Pamir drained by a number of endorheic river systems such as Amudarya, Syrdarya, Ili, Tarim and a few smaller ones. The Central Asian river basins are characterized by the semi-arid climate with strong seasonal 15 variation of precipitation. Most precipitation falls as snow during winter and spring months in Western and Northern Tien Shan (Aizen et al., 1995(Aizen et al., , 1996Sorg et al., 2012). In contrast, Central and Eastern Tien Shan receive their largest precipitation input during the summer months. The Pamir mountains receive the highest portion of precipitation during winter and spring months with minimum in summer (Schiemann et al., 2008;Sorg et al., 2012).
Precipitation also exhibits a high spatial variation, with e.g. less than 50 mm/year in the desert areas of Tarim and around 100 20 mm/year on leeward slopes of Central Pamir to more than 1000 mm/year in the mountain regions exposed to the westerly air flows being a major moisture source in the region (Aizen et al., 1996;Bothe et al., 2012;Hagg et al., 2013;Schiemann et al., 2008). This makes Tien Shan and Pamir Mountains the regional 'water towers', with snow melt to be the dominant water source during spring and early summer months. During summer, glacier melt and liquid precipitation gain importance depending on the basin location and degree of glacierisation (Aizen et al., 1996). The Tien Shan and Pamir mountains exhibit 25 particularly high relative water yield compared to the lowland parts of these catchments (Viviroli et al., 2007). Related to the economic water demands in the lowland plains primarily for irrigated agriculture, the Tien Shan and Pamir mountains are among the most important contributors of stream water worldwide (Viviroli et al., 2007). They are also among those river basins with the highest share of glacier melt water in summer, particularly in drought years (Pritchard, 2017). Within the Aral Sea basin, to which the Amudarya and Syrdarya rivers drain, the irrigated area amounts to approximately 8.2-8.4 million ha 30 (Conrad et al., 2016;FAO, 2013) Additionally, considerable irrigation areas are located in the Aksu/Tarim basin, where agricultural land doubled in the period 1989-2011 and land use for cotton production increased even 6-fold (Feike et al., 2015). Irrigated agriculture in Central Asia (CA) is mainly fed by the stream water diversion with only small portion of groundwater withdrawal (FAO, 2013;Siebert et al., 2010). Hence, reliable prediction of seasonal runoff during vegetation period (April -September) is crucial for agricultural planning and yield estimation in the low lying countries in the Aral Sea basin, as well as for the management of reservoir capacities including dam safety operations in the upper parts of the catchments. Seasonal forecasts are one of the major responsibilities of the state hydro-meteorological (hydromet) services of the Central Asian 5 countries and are regularly released starting from January till June with the primary forecast issued end of March -beginning of April for the upcoming 6-months period. In some post-soviet countries, these forecasts are typically developed based on the empirical relationships for individual basins relating precipitation, temperature and snow depth/SWE records to seasonal discharge, partly available only in analogue form as look-up tables or graphs (Hydromet Services, unpublished questionnaire survey undertaken within the CAWa project). Particularly, point measurements of snow depth and/or snow water equivalent 10 (SWE), which have been carried out by helicopter flights or footpath surveys in mountain regions in the past decades, are costly or not feasible due to access problems nowadays. Other Hydromet Services apply the hydrological forecast model AISHF (Agaltseva et al., 1997) developed at the Uzbek Hydro-meteorological Service (Uzhydromet), which computes discharge hydrographs by considering temperature, snow accumulation and melt. Snow pack is accumulated in winter and temperature and precipitation are taken from an analogous year to drive the model in the forecast mode. Hydrometeorological 15 services rely on the available meteorological and hydrological data acquired by the network of climate and discharge stations, which, however, strongly diminished during the 1990s  and slowly recovers nowadays, partly with substantial international support (e.g. Schöne et al. (2013); CAHMP Programme by World Bank; previous programmes by SDC and USAID). Hence, to fulfill their task, hydrometeorological services need the timely to near real-time data and simple methodologies capable of utilizing available information. 20 Schär et al. (2004) showed the potential of the ERA-15 precipitation data from December-April period to explain about 85% of the seasonal runoff variability in May-September in the large-scale Syrdarya river basin. The explained variance for the Amudarya River amounted, however, to only about 25%, presumably due to poor precipitation modelling in the ERA dataset, strong influence of glacier melt and water abstraction for irrigation purposes. Similarly, Barlow and Tippett (2008) explored the predictive power of NCEP-NCAR cold-season (November-March) precipitation for warm-season (April-August) discharge 25 forecast using canonical correlation analysis. Though for some of the 24 Central Asian gauges, no skillful prediction could be achieved, for a few catchments 20 to 50% explained variance could be attained. Archer and Fowler (2008) utilized temperature and discharge records additionally to precipitation for spring and summer seasonal flow forecast on the southern slopes of Himalaya in northern Pakistan using multiple linear regression models. Despite good predictions of spring and early summer flows, late summer discharges were poorly forecasted due to the strong influence of summer monsoon. Recently, Dixon and 30 Wilby (2015) demonstrated the skill of a linear regression model for the Naryn basin, Kyrgyzstan, based on TRMM precipitation from October-March to explain 65% of the seasonal flow variance in the vegetation period. The authors selected specific TRMM pixels in the catchments showing the highest correlation to seasonal discharge. They also explored the predictive skills of multiple linear regression models additionally including temperature and antecedent discharge and testing different lead times from one to three months. They showed that forecasts based on multiple linear regression models are always superior to zero order forecasts, i.e. the mean flow.
The fact that substantial snow accumulation in Central Asian mountain regions during the winter and spring months significantly governs runoff in the vegetation period can be effectively utilized for seasonal forecasts. For a similar climatic setting, Pal et al. (2013) included the measurements of snow water equivalent at point locations into multiple linear regression 5 models along with precipitation, antecedent discharge and temperature-based predictors. Linear models with multiple predictor combinations achieved skilful forecasts of the spring (March-June/April-June) seasonal flow in northern India on the southern Himalaya slopes. Point snow measurements are, however, rarely available and remotely sensed snow cover extent can provide a viable alternative. Based on the monitored snow cover extent, e.g. using optical satellite imagery, and additionally considering temperature and precipitation to implicitly approximate snow water equivalent (SWE) a solid basis for seasonal 10 discharge forecast can be formed. The MODIS snow cover product was shown to deliver high accuracy for the Central Asian region . Methodologies to remove cloud obstruction of optical imagery have matured over the past decade (Gafurov and Bárdossy, 2009;Gafurov et al., 2016) and tools for the automated image acquisition and processing reached the operational level . MODIS snow cover data was used for runoff forecast in the Argentinian Andes in the high-water season (September-April), though no cloud elimination algorithms were applied (Delbart et al., 2015). Snow cover 15 in September-October could explain about 60% of the high-water season discharge variance. However, no skilful forecast with lead times greater than zero were possible. Rosenberg et al. (2011) proposed a hybrid (statistical -hydrological model) framework for seasonal flow prediction in Californian catchments using accumulated precipitation in antecedent period and SWE modelled by a distributed hydrological model. These two predictors were linked to seasonal discharge by principal component and Z-score regression (Rosenberg et al., 2011). The hybrid approach was found comparable and in some cases 20 superior to a purely statistical approach, however, at the cost of effort for hydrological simulation of the SWE dynamics.
Based on the finding of the studies listed above, we propose a simple methodology for the operational forecast of seasonal runoff for the vegetation period (April-September) for all Central Asian catchments, which areas range over three orders of magnitude. The method is based on multiple linear regression models with automatic predictor selection, whereas the predictors are based on the readily available precipitation, temperature and discharge gauge records and additionally leverage 25 by the operationally processed cloud-free MODIS snow cover product. It is argued, that in linear modelling the use of meteorological data from a single gauging station for a large catchment is justified, as long as the variability of the station records are representative for the whole catchment. We demonstrate the model predictive skill and robustness in a crossvalidation and discuss the relative importance of the automatically selected predictors depending on the prediction lead time.

Study sites and data
For the testing of the forecast models 13 catchments were selected. They cover a wide range of geographical regions, ranging from catchments along the western slopes of the Altai mountains in Eastern Kazakhstan (Uba, Ulba), over catchments at the western and northern rim of the Tien-Shan (Chirchik, Talas, Ala-Archa, Chu, Chilik, Charyn) and central Tien-Shan (Karadarya, Naryn) mountains (cf. Aizen et al., 2007), to the northern and central Pamir (Amudarya) and the northern 5 Hindukush (Murgap). The size of the catchments varies over three orders of magnitude from 239 km 2 to 288,000 km 2 . Figure   1 provides an overview of the location and size of the catchments, while Table 1 additionally lists the discharge and meteorological gauging stations used for the seasonal flow forecast. Note that the Naryn catchments are nested. The Upper Naryn represents a high alpine catchment and is the headwater catchment of the larger Naryn catchment draining into the Toktogul reservoir. This separation was undertaken in order to test the proposed method also for a high alpine catchment with 10 a comparably high degree of glacialisation. The wide range of catchment locations, climatic conditions and sizes enable a testing of the proposed forecast models under different boundary conditions, and thus provides an indication of the applicability, robustness and transferability of the approach.
The catchment boundaries are derived to map the catchment area draining to the selected discharge stations. For the meteorological data (temperature and precipitation) meteorological stations run by the individual Hydromet Services were 15 selected. Ideally those are located in the catchment area and have sufficient data coverage of at least 16 years (starting in 2000 in order to be consistent with the MODIS temporal coverage). However, in some catchments meteorological stations fulfilling these criteria were not available. For those catchments stations nearby were selected for the prediction.  For both discharge and meteorological data monthly values were obtained for the stations listed in Table 1, i.e. monthly mean discharges, monthly mean temperatures and monthly precipitation sums. For the presented study meteorological station data was used, because of the operational availability to the CA Hydromet Services. Gridded re-analysis products like ERA-Interim 5 typically have a latency of weeks to months, and thus cannot be used for operational forecasts to fulfil the mandatory regulations while station temperature and precipitation data are likely not representative for basin average values, it is assumed that the variability of the catchment averages and the station data is similar. This, in turn, enables the use of the station data in the statistical forecast using multiple linear regressions.
In addition to the station data, mean monthly snow coverages for the individual catchments were calculated using daily snow 10 cover data derived by the MODSNOW-Tool (Gafurov and Bárdossy, 2009;Gafurov et al., 2016). MODSNOW uses the MODIS satellite snow cover product and applies a sophisticated cloud elimination algorithm (Gafurov and Bárdossy, 2009;Gafurov et al., 2016) to obtain cloud free daily snow cover images. The MODSNOW-Tool runs operationally in most of the CA Hydromet Services, thus enabling the use of snow cover information for operational forecasts.
Due to the use of MODIS snow cover, which is available since March 2000, the time series of the data used for the construction 15 of the forecast models had to be limited to post-2000. The time period for the model development and testing was thus set to 2000 -2015, for which mostly complete continuous time series for all data and stations were available.
The seasonal discharge, i.e. the predictand of the forecasts, is calculated as the mean monthly discharge for the period April to September. Figure 2 shows the seasonal discharges for all catchment considered in this study. The top panel highlights the differences in the magnitude of the seasonal discharge, spanning almost three orders of magnitude (cf. also Table 1). Discontinuous lines indicate data gaps. In order to illustrate differences in the inter-annual variability of the seasonal discharge the lower panel of Figure 2 plots the seasonal discharges normalized to zero mean and standard deviation of 1. This plot indicates similar but also 5 different inter-annual variability patterns of the different catchments. In order to distinguish between similar and different inter-annual variabilities cross-correlations of the seasonal discharges are calculated and hierarchically clustered (Figure 3).

Seasonal discharge variability
Cluster memberships were established using the Ward algorithm clustering the catchments based on the dissimilarities of the correlation between the seasonal discharge time series of the different catchments. The correlation matrix in Figure 3 shows that the seasonal discharges mainly cluster according to their geographical location. The variability of the seasonal discharge 10 of the two catchments in the Altai region (Uba, Ulba) is distinctively different to all the others. Also the two most southern catchments (Amudarya and Murgap) form a distinct cluster that is joined by the most western catchment of the northern Tien Shan, Chirchik. However, Chirchik is also well correlated to the largest group, the catchments in the Tien Shan, which all show similar inter-annual variability of the seasonal discharge. An exception to this is the smallest catchment in the study, Ala-Archa, which is not correlated to any of the other catchments, presumably due to the strong influence of local meteorology 15 and glacier-melt dominated discharge formation in the summer months.
The analysis of the inter-annual variability thus maps the geographical and climatic differences of the catchments considered in this study to a large extent. These differences in variability, but also in the magnitude of the discharges and catchment size imply that the forecast methods can be tested against a wide range of boundary conditions and seasonal variabilities. If skilful forecasts are obtained for all catchments, it can be argued that the approach delivers robust forecasts that are not obtained by 20 chance or due to similar variabilities in all catchments. If successful, it could also be inferred that the approach can be transferred to other regions with similar streamflow generation characteristics.

Method
As mentioned in the introduction, the seasonal discharge during the vegetation period of April to September in CA is dominated by snow melt in the mountains. Therefore a good estimation of the snow accumulation and snow water equivalent in the catchments during the winter months may provide reliable forecasts of the discharge during the vegetation period. However, 10 data about the depth and snow water equivalent are not regularly acquired except for some dedicated research sites. Thus alternative data containing proxy information about the snow depth and water equivalent must be used. Therefore predictors for the forecast models were derived from mean monthly temperature records, monthly sums of precipitation and monthly mean snow coverage of the catchments. It is argued that the combination of these factors is able to serve as proxy data for snow depth and water equivalent. While the precipitation directly contains information about the snow fall amount and thus 15 accumulation, temperature may contain information on the wetness of the snow pack. In combination with snow coverage, temperature and precipitation may thus provide information about the snow volume and water content. In addition to the climate data monthly antecedent discharge can serve as an indicator about the magnitude of the snow melt process and groundwater storage state and release, and is used as predictor, too.
For some regions, particular the Altai catchments, early summer (May -July) precipitation plays a larger role for the seasonal discharge. This precipitation is partly considered as observations in the late forecasts presented here. However, reliable 5 information about the spring precipitation in advance could possibly improve the early forecasts. But due to the low predictability of the typically convection type summer precipitation (Gerlitz et al., 2016) this is not considered in the predictor set.
Evaporative losses in the presented mountain catchments are considered low due to the low summer temperatures and fast catchment response and high water flow velocities in the rivers. Higher losses can occur in reservoir lakes, but with the 10 exception of the large Amudarya basin there are no reservoirs present in the selected catchments. In the Amudary catchment the Nurek reservoir lake at the Vakhsh river exist. However, evaporative losses from the lake surface area of 98 km 2 can be considered negligible in comparison to the large catchment size. Therefore evaporation is not directly considered as predictor for the forecasts.
Moreover, the catchment area of the Vaksh river at the conjunction with the Panj river amounts to 31,415 km², equivalent to 15 about 11% of the Amudarya catchment at Kerky considered here. Assuming further that the reservoir can manage only a fraction of the total discharge of the Vakhsh river, and that the effects of the water retention are further buffered by the seasonal mean discharge spanning six months, it can be assumed that the regulating effect of the Nurek dam on the overall seasonal discharge of the Amudarya at Kerky is rather low. Additionally, the dam is operational since 1980, therefore a discontinuity in the time series 2000-2015 can be ruled out. We thus argue that the anthropogenic influence of the seasonal discharge time 20 series of the Amudarya is negligible for the presented study.

Generation of the predictor set
The core set of predictors consists of the monthly values preceding the prediction date. According to the operational forecast schemes of the CA Hydromet Services a series of different prediction dates were defined. The first prediction of the seasonal mean discharge for the vegetation period (April to September) is issued on January 1 st , followed by predictions on February 25 1 st , March 1 st , April 1 st , May 1 st , and June 1 st . The predictions January to March are preliminary forecasts, while the prediction on April 1 st is the most important for the water resource planning in the CA states. The following forecasts serve as corrections of the April forecast. They are actually partial hindcasts, as the predictors already cover a part of the prediction season. This means that the predictors for the late forecasts are not fully independent from the predictand. This decision was on the first hand made due to the foreseen implementation of the method in the Central Asian Hydromet Services. The presented procedure 30 is in line with the official forecast procedures in the Central Asian Hydromet Services. In order to obtain acceptance of the proposed method in the services and their use in the official forecast procedures it is advisable to follow the prescribed procedures. It is required from the Hydromet Services to issue updated (corrected) forecasts, which include the entire vegetation period (April-September), The water regulation procedures and e.g. agricultural yield estimation are traditionally based on bulk numbers for the entire period. If these procedures are not followed, the obtained results, which are better than the forecasts issued with the existing procedures, might not be implemented and come into practise, and thus a chance would be missed to bring research results into application. Furthermore, shortened seasonal discharge time series, where the predictors and predictand are independent, are highly correlated (cf. Annex 1). This means that results obtained with seasonal discharges 5 calculated only for months following the prediction data will obtain very similar results to the ones using the full discharge time series.
For the prediction up to the 1 st of April the monthly values over the whole winter period, i.e. from October onwards are used.
For later predictions this was limited to data of the prediction year, i.e. from January onwards, in order to keep the number of predictor combinations in reasonable limits. The monthly predictor values were accompanied by multi-monthly means, 10 spanning over two and three months prior to the prediction date, and mean values for the whole predictor period defined above, i.e. either from October to the prediction month, or from January to the prediction month, respectively. Furthermore, composites were calculated from the climatological data in order to extend the predictor set. They are introduced in order to explore their potential to reflect snow wetness. It is argued that composites can improve the prediction by linear models, as some non-linear interactions might be reflected better by composites compared to the raw data (as shown in e.g. 15 Hall et al., 2017). Analogously to the original data, monthly and multi-monthly composites were derived. For the composites, products of "temperature and precipitation", "temperature and snow coverage", "precipitation, snow coverage and temperature", "precipitation and snow coverage" were used. Antecedent discharge was not included in the composites, because this should not influence the snow cover characteristic.

Statistical modelling 20
For the development of the statistical forecast models standard multiple linear regression (MLR) was applied. It is argued that the discharge generation from snow melt over whole catchments and on a seasonal time scale can be approximated by linear models. In fact, this was shown by a large number of studies using hydrological models based on linear concepts like linear storage, e.g. Duethmann et al. (2014) and Duethmann et al. (2015) in Central Asia. Additionally a number of studies have shown that linear regression is a valid approach for seasonal forecasts (e.g. Delbart et al., 2015;Dixon and Wilby, 2015;Seibert 25 et al., 2017). A linear modelling approach is thus seen as a valid approach for seasonal forecasts in the study region from a general point of view. However, in order to statistically support the assumption that the runoff generating processes can be approximated by linear models, the formal assumptions of MLR were also tested: the assumption of normal distribution of the residuals was tested by the Shapiro-Wilk test, the independence of the residuals was tested by calculating the autocorrelation with lag 1, and the heteroscedasticity of the residuals was tested by the Breusch-Pagan test. 30 In the model selection procedure all possible predictor combinations, which are different for every prediction month as described in 3.1, are used in the MLR for the construction of forecast models. However, some restrictions were put on the predictor combinations in order to avoid overfitting and thus spurious regression results: 1. The predictors are grouped into 8 groups: snow cover, temperature, precipitation, antecedent discharge, and the four composite types.
2. The maximum number of predictors in a regression is limited to four.
3. Only one predictor from each group of predictors can be used in an individual regression model. This resulted in 7,728 predictor combinations, i.e. multiple linear models to be tested in January, and increased to 155,690 5 possible models in April. A complete list of the predictors for the different prediction months is provided in the Annex. The coefficients for all these linear models were automatically fitted during the MLR by the least squares method. Only models with all predictors statistically significant at p = 0.1 and with an overall model significance of at least p = 0.1 were retained.
The best models were selected based on the lowest Predicted Residual Error Mean of Squares (PREMS) value obtained by a Leave-One-Out Cross Validation (LOOCV). In the LOOCV one year of the seasonal discharge time series is removed from 10 the data set for fitting the MLR. The missing data point is then estimated by the model fitted to the remaining data. The PREMS value is the mean of squared errors of all seasonal discharges left out and the associated predicted LOOCV values. PREMS is thus defined as: with e(i) being the residuals of the LOOCV: 15 where � ( ) is the regression estimate of based on a regression equation computed leaving out the i th observation of the overall number of n observations. The PREMS was used in this study instead of the usual PRESS (Predictive Residual Error Sum of Squares) in order to avoid biases possibly introduced by missing predictor or predictand data. Using the sum of squares could favour models with missing data compared to models providing predictions for all 16 years. Using the mean of the squares can 20 avoid this to a large extent.
The LOOCV is testing the MLR for robustness and can avoid overfitting and incidental good MLR results valid for the whole data set only. In order to avoid an over-estimation of the forecast skill the seasonal discharge time series were tested for autocorrelation, which could lead to spurious estimation of model robustness.
Model skill was evaluated in a number of measure: adjusted R2, root mean square error RMSE, and mean absolute error MAE. 25 The robustness of the model ensemble was quantified as the ration of the adjusted R2 based on the LOOCV residuals to the adjusted R2 of the complete model residuals. The reliability of the model was analysed by PIT diagrams (e.g. Crochemore et al., 2017) and quantified as PIT-scores (Renard et al., 2010).
In the presented study not only the single-best model according to PREMS of the LOOCV-MLR was selected as prediction model, but rather the best 20 models, if more than 20 models pass the significance tests. This selection aims at the analysis of 30 the differences between the best models in terms of performance and predictors, but also serves as a model ensemble for the forecast of the seasonal discharge. The distribution of the residuals of the best 20 forecast models was evaluated to provide 80% predictive uncertainty intervals for every forecast. However, it has to be noted that the choice to use the best 20 models is subjective, and this number can be increased or reduced. Moreover, a different set of specific models from the best models can be selected according to their performance measures and temporal dynamics by experts knowledgeable of the individual catchments. Sufficient amount of freedom was left for the selection of the number of best models to be retained, in order to 5 enable an expert selection of models by the forecasters of the Central Asian Hydromet Services. The forecasters check every model retained for their performances (quantitatively and qualitatively), and select the models accordingly.

Predictor importance
The predictors of the selected best models were analysed for their importance, i.e. their share of the overall explained variance 10 (R 2 ) of the individual models. This was achieved by the lmg algorithm implemented in the R-package relaimpo (Grömping, 2006). lmg is based on sequential R 2 s, but explicitly eliminating the dependence on predictor orderings by averaging over orderings using simple unweighted averages. In sequential R 2 calculations, the model is re-run with a single predictor only and the explained variance is calculated. Then the next predictor is added and the gain in explained variance is calculated. By this procedure the variance explained by individual predictors can be quantified. However, in this procedure the sequence of 15 predictors added influences the share of explained variance associated to the individual predictors. Therefore the lmg algorithm tests all possible predictor sequences and calculates the mean importance of every sequence in order to overcome the problem of predictor ordering in sequential R 2 s. The predictor importance calculation yields information about the importance of the individual predictors at different forecast points in time for the catchments under study, which in turn can be used for a discussion of the factors responsible for the winter snow accumulation and snow water content in the catchments. 20 However, such a discussion is complicated by the use of the composite predictors. Therefore the importance of composite predictors is divided into equal proportions to the components of the composites. If more than one composite is used in a model, the proportions associated to the component factors (snow cover, precipitation, temperature) are summed up and displayed as parts of the composite importance in the figures presented in 4.2. This analysis is not meant to provide a quantitative estimation for the component importance of the composite predictors, but rather to enhance the discussion and 25 interpretation of the predictors of the selected forecast models.
In addition to the importance for an individual model (here the best LOOCV model), the mean importance over the best 20 LOOCV models is calculated. This is achieved by calculating the fractions of the sum of importance of an individual predictor for all 20 models to the sum of the R 2 values of all 20 models for each catchment and month. These fractions are then multiplied by the mean R 2 values of the best 20 models. This mean predictor importance can be compared to the predictor importance of 30 the best model in order to analyse the stability of the predictor selection within the best 20 LOOCV models.

Results
In order to test the suitability of the LOOCV for the seasonal streamflow forecast the autocorrelation and partial autocorrelation of the streamflow time series was calculated and plotted (Annex 2). Any autocorrelation in the discharge time series could lead to artificial over-estimation of the forecasts skill by the LOOCV. Hardly any autocorrelation at α = 0.05 could be detected.
Only for the Ulba some significant autocorrelation for lag 1 and 2 is shown just above α = 0.05, but by the partial autocorrelation 5 only. No autocorrelation was found at α = 0.01. Therefore it can be stated that autocorrelation does not exist in the discharge time series of all catchments, and thus the proposed LOOCV is an appropriate validation method.
The MLR fitting with LOOCV (cf. 3.2) was applied for different forecast dates ranging from January 1 st to June 1 st for all catchments. The best 20 models according to the PREMS resulting from the LOOCV were retained for the forecasts. The tests for possible violations of the formal MLR assumptions showed, that 89.5% of all models for all catchments and prediction 10 months fulfilled the criteria of normal distributes residuals, 95.8% of the models passed the test for independence of the residuals, and 99.5% of the models have homoscedastic residuals (cf. Annex 3). In summary the formal requirements of MLRs are fulfilled by almost all models, and the use of linear models for seasonal discharge forecast is justified also from a formal point of view.
In general the performance of the linear models increases from January to June, with the best models reaching adjusted R 2 15 (adj. R 2 ) values in the range of 0.76 -0.89 in April and 0.84 -0.99 in June. For most of the catchments high adj. R 2 values in the range of 0.57 -0.78 were already obtained in January. Only for Amudarya, Chu and Chirchik the performance is unsatisfyingly low in January, but increases to adj. R 2 > 0.59 already one month later in February. Table 2 lists the adj. R 2 values of the best LOOCV models for all catchments and forecast months. Note that the adj. R 2 in the table are calculated using the coefficients of the linear models fitted to the whole data set, i.e. they are not cross validated. 20 While for most of the catchments, the performance of the models gradually increases with decreasing lead time, the performance for Chilik, shows significant decreases and increases. This is caused by a comparatively large number of missing discharge and predictor data. The automatic fitting algorithm takes advantage of this by finding models able to explain the fewer data points better compared to the full time series despite the use of PREMS instead of PRESS. However, these models can already represent an overfitting and are thus less reliable or stable in time compared to models fitted to longer time periods . 25 In order to get a more comprehensive picture of the model performance, Figure 4 shows the temporal evolution of the adj. R 2 evaluated for the complete time period of the single best LOOCV model, the minimum adj. R 2 of the best 20 models, the mean adj. R 2 of the best 20 models, the root mean square error (RMSE) of the single best LOOCV model calculated for the full data set normalized to the mean seasonal discharge (cf. Table 1), the normalized mean absolute error MAE, and the PREMS value of the best model. Note that the highest adj. R 2 value is not necessarily the adj. R 2 of the single best model, because the best 30 model is selected according to the lowest PREMS in the LOOCV, and not the best adj. R 2 evaluated using the whole time series. Therefore the mean adj. R 2 in January is occasionally higher than the adj. R 2 of the best LOOCV model, i.e. the most robust model. In general, Figure 4 shows that the different adj. R 2 , RMSE, MAE and PREMS values are similar in their evolution in time, i.e. increase (adj. R 2 ), resp. decrease (RMSE, MAE, PREMS) with later forecast months. This indicates that for all best 20 models the performance is improving with later forecasts. Furthermore, the difference between min adj. R 2 and mean adj. R 2 to the adj. R 2 of the single best LOOCV model is typically larger in the early prediction months. This indicates a wider spread of model performance within the selected 20 models for the predictions with longer lead times. This difference decreases with shorter lead times, meaning that more models with 5 similar high performance can be found, and thus uncertainty of the model ensemble is reduced. To a certain extent this is likely caused by the larger number of possible predictors for later prediction months, but it is also well justified to assume that the later predictors have more predictive power: data from the late winter months can better describe the snow coverage and water content compared to predictors from the previous autumn. This issue will be discussed further in Section 4.3. Figure 4 shows that the RMSE as well as the MAE of the best model of the LOOCV is at maximum about 35% of the long 10 term seasonal mean discharge (Talas in January). However, for most catchments the normalized RMSE and MAE is below 20% in January already. For the important April forecast they fall generally below 10%, except for Talas and Murgap, where it remains at 20%. These values state the high performance of the linear forecast models in terms of actual discharge, and are thus a useful information for practitioners in order to assess the value of the forecasts. PREMS values generally decrease (i.e. improve) with decreasing lead time. However, occasionally increases can be observed for later forecast months. This can be also seen in the adj. R 2 values, but less pronounced because of the scale of the left yaxis. This phenomenon is caused by the changing predictor sets from forecast month to forecast month. Particularly multimonthly predictors change for each prediction date according to the parameter selection outlined in Section 3.1. As this phenomenon of increasing PREMS values usually occurs in April or May, it can be hypothesized that the information of the 20 late winter/early spring months used in the later forecasts does not contain better information about the snow cover as the previous months. With respect to a practical application, the better performing forecasts from the previous months can be used, which is equivalent to an extension of the predictor set by the predictors of the previous month.
This general reduction of PREMS also means that the models become more robust for later prediction months. To illustrate this more clearly, Figure 4 also shows the relation between the mean adj. R 2 of the LOOCV for all 20 models to the mean adj. 25 R 2 of the full model fit. The mean adj. R 2 of the LOOCV is calculated from the LOOCV residuals used to calculate the PREMS.
According to the rationale of the LOOCV, a model is more robust and less prone to overfitting, if the LOOCV-R 2 is very close to the overall R 2 . Figure 4 shows that this is generally the case for the catchments with very high adj. R 2 values, and also for later prediction months. This means that the selection of the predictors is likely stable even if additional data is added to the time series in future. However, there are some catchments for which comparably less robust models could be derived even for 30 later prediction months (5. Ala-Archa, 6. Chu). For these catchments it is likely that the predictor selection will change with additional data. -values of the best performing prediction models from the LOOCV for all catchments and prediction months. "best" indicates the single best model according to the LOOCV, "mean" indicates the mean percentage over the best 20 models according to the LOOCV. The adjusted R 2 values are associated with indicators for significance levels.

January
February March April May June best mean best mean best mean best mean best mean best mean  Figure 4: Performance of the prediction models for the different catchments and prediction months. Adj. R 2 best model is the adjusted R 2 of the single best LOOCV model, mean adj. R 2 is the mean adj. R 2 of the best 20 LOOCV models, min adj. R 2 is minimum adj. R 2 of the best 20 LOOCV models, robustness is mean LOOCV-adj. R 2 of the best 20 models divided by the mean adj. R 2 , RMSE/MAE norm. is the root mean squared error/mean absolute error of the single best model normalized to mean multi-annual seasonal discharge, mean RMSE/MAE norm is the mean root mean square error/mean absolute error of the best 20 LOOCV models normalized to the multi-annual seasonal discharge; PREMS is the predictive residual sum of squares (PRESS) of the single best model, divided by the number of prediction months (i.e. mean of squares).

5
In addition to the performance metrics Figure 5 plots the temporal dynamics of the best LOOCV models for all six prediction months. It can be seen that the models can map the high variability of the observed seasonal discharges very well, often already in January or February. This graphically corroborates the findings derived from the performance metrics and underlines that the good performance of the models is not a statistical artefact.

Figure 5: Forecasts of the seasonal discharge by the single best model selected by the LOOCV for the individual catchments and all prediction months. The blue lines show the observed seasonal discharges. Note that some models do not provide forecasts for every year due to missing predictor data.
In order to set the performance of the presented models in the context of the routines and guidelines of the Central Asian Hydromet Services, the performance of the models was also estimated according to the performance criteria used by the Hydromet Services. This is defined by: 5 With |res| denoting the absolute value of the residual of an individual forecast, and σQs the standard deviation of the seasonal discharge (here calculated for the discharge time series used, i.e. for the period [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015]. According to the protocols of the Hydromet Services an acceptable ("good") forecast is defined by Sσ < 0.675. Table 3 shows how often this criteria was fulfilled during the analysis period 2000-2015 for the best model, and on average by the best 20 models. For the critical forecast month 10 April the criteria was fulfilled for at least 82% of the years (13 out of 16 years) by the model ensemble for most of the catchments. For Ala-Archa and Chu the numbers were lower, but still as high as 75% and 77%. For all catchments the percentages increase further for the later forecast months. These findings are also valid for the ensemble of the best 20 models, as the very similar percentages of the mean of all models compared to the best model indicate. This means that the developed models would provide acceptable forecasts for the Hydromet Services in the range of 80%-90% for the important forecast 15 month April.

Predictive uncertainty
In order to quantify the predictive uncertainty the empirical 10% and 90% percentiles of the residuals of the forecasts ensembles consisting of the up to best 20 models according to PREMS were calculated for every prediction month. Note that for the early prediction months occasionally less than 20 models passed the significance test. The tables in Annex 3 indicate when this was 5 the case. The quantiles of the residuals were then added to the actual model predictions, thus providing an 80% predictive uncertainty band, i.e. an interval in which the true value of the seasonal discharge should lie with a probability of at least 80%. Figure 6 shows the predictive uncertainty bands for every catchment along with the observed seasonal discharge. The predictive uncertainty for the different prediction months are shown in shades of orange. In general it can be seen, that the predictive uncertainty bands narrow with later prediction months, illustrating the better prediction during later prediction 10 months described above. While this is perfectly visible for most catchments (e.g. 3. Chirchik, 7. Karadarya), it is not the case for some others (5. Ala-Archa, 6. Chu, 10. Naryn). The main reason for this is the larger difference between the predictions and performance of the best 20 models compared to the other catchments, as indicated by the difference between the best and mean adj. R 2 shown and listed in Figure 4 and Table 2, respectively. This causes a wider distribution of the residuals of the best 20 models and thus higher predictive uncertainty. However, if only the best or a smaller selection of the best 20 models 15 are used for a forecast, the predictive uncertainty would also be reduced. This means, that the uncertainty bands derived depend on the subjective choice of the number of models to be kept in the model ensemble. Another reason for wider predictive uncertainty bands for later months is the observed decline in performance during later months in some catchments due to the changed predictor set (e.g. for 6. Chu). This causes again higher predictive uncertainty bands, which overlay the narrower band from the previous month. 20 From a formal point of view the uncertainty bands correctly include at least 80% of the observed seasonal discharges, even for very narrow bands (e.g. in June for 3. Chirchik or 9. Karadarya). This indicates that the uncertainty estimation derived from the regression residuals provide a reliable estimation of the uncertainty associated to model selection, and can be used to derive decisions based on the forecasts given by the MLR model ensembles. However, it must be noted that the derived uncertainty bands represent the predictive uncertainty of the MLR models fitted to the available time series. They do not account for any 25 uncertainty stemming from a possible lack of representativeness of the rather short time series used. Longer discharge time series might show a different variability of seasonal discharge, which would then not be covered by the derived models.
However, as the models can be updated every year in future, this potential problem is expected to decrease with further use of the approach in the Central Asian Hydromet Services.
However, it has to be noted that the estimated uncertainty cover only the model selection uncertainty. Other uncertainty sources 30 are: • model structure, which is assumed to be rather low given the high explained variances; • data sources, which is not quantifiable, but might be high, particularly the discharge data; • and the performance criteria for selecting the best models.
The last aspect has been tested. Using other performance criteria as PREMS can result in a slightly different selection of best models, but more often just in a different order of the best models. The best PREMS model is not necessarily the best cross validated R 2 model, or the best MAE or RMSE model. However, as this mainly affects the ordering of the best models, the results in terms of ensemble predictions and predictive uncertainty, if unweighted as presented, would be very similar. 5 Figure 6: 80% predictive uncertainty bands for all catchments and forecasts months. The blue lines indicate the observed seasonal discharges.
In addition to the predictive uncertainty also the reliability of the forecasts was quantified by PIT diagrams and PIT scores. Figure 7 shows the PIT diagrams for every catchment and all forecast months using the forecasts of the selected ensemble models. The PIT diagrams show that the model ensemble predictions are in most cases close to the 1:1 line, i.e. provide reliable forecasts. However, in some cases the predictive uncertainty is under-estimated (PIT diagram lines with pronounced vertical component around the 50% quantile). This means that some of the predictive uncertainty bands presented in Figure 6 are too 5 narrow to reliably quantify the predictive uncertainty. This is mostly the case for the late forecasts with high skill, where the models in the ensemble often produce very similar forecasts. In addition to the diagrams a PIT score was calculated as the area between the PIT curve and the 1:1 line as a summarizing indicator for the reliability (Renard et al., 2010). The theoretically least reliable model has a score of 0.5, a perfect model a score of 0. The highest score, i.e. the lowest reliability, of all models is 0.2, with the majority of the models being in the range of 0.07-0.15. Interpreting the scores with the curves in the PIT 10 diagram it can be deduced that the reliability of model ensembles with PIT scores ≤ 0.1-0.14 is acceptable. For higher scores the predictive uncertainty derived from the model ensemble is likely to be underestimated.  Figure 8 illustrates the importance of the predictors of the selected MLR models as absolute fractions of the R 2 values, whereas it is not differentiated between individual predictors, but rather between predictor classes described in 3.1. The left panel of Figure 8 shows the importance for the single best LOOCV model, while the right panel shows the average importance of the 5 predictors for the best 20 LOOCV models. A comparison of the left and right panels shows that the predictor selection and importance for the different catchments and prediction months of the best model is quite similar to the mean of the best 20 models. This indicates that the predictor selection for the models in the ensemble is quite stable, and hence that the predictor selection is not random, but rather has some hydrological meaning. However, an interpretation of the contributions of the different factors is complicated by the use of the composites, which are almost always selected as one or more predictors in 10 the MLR models. Nevertheless, some general features can be identified from Figure 7:

Predictor importance (Is there some hydrological process information in linear models?)
• Typically there is no single factor dominating the explained variance, with the exception of Karadarya, where the composites have an exceptionally large share on the explained variance. But as the composites are comprised of the other predictors (except antecedent discharge), this statement is actually valid for all catchments. This indicates that the winter snow accumulation providing the bulk of the seasonal discharge is best described by a combination of the 15 factors determining the extent and water equivalent of the snow pack in the catchments (precipitation, temperature, snow coverage). Omitting one of these predictors leads in fact to a reduction in model performance.
• There is a general and plausible trend for higher importance of antecedent discharge in the later prediction months.
In this period it can be expected that antecedent discharge has higher predictive power of the seasonal discharge compared to the winter months, i.e. during the accumulation phase, because it directly indicates the magnitude of the 20 discharge generation from snow melt. This finding is valid for most catchments except Chirchik, Ala-Archa and Chilik. For Chirchik the importance of antecedent discharge is almost constant throughout the prediction months, both for the best model and on average. Contrary to this, antecedent discharge has very little importance for Ala-Archa and Chilik. For Ala-Archa this observation can be explained by the very small catchment size and thus the quick response of discharge to precipitation events and lower transit times, but also with the high proportion of glacier 25 melt during the summer months. The high importance of precipitation, which is higher than in any other catchment particularly in the later prediction months, also supports this reasoning. For Chirchik and Chilik, however, no plausible explanation can be derived from the basic catchment characteristics presented here.
• The importance of the snow coverage predictors indicate a regional differentiation of the predictor importance. For the two catchments in the Altai region (1. Uba, 2. Ulba, cluster 1 in Figure 3) snow coverage as an individual factor 30 is of less importance compared to the other regions. This is due to different snow cover characteristics in these catchments, which have comparatively lower altitudes compared to other catchments in this study. Therefore, snow accumulation in these catchments is comparably low and quickly responds to increasing temperature already in the spring months. Seasonal snow cover variations obtained from MODSNOW-Tool  for these catchments also illustrate sudden snow cover depletion in the month of April for both catchments, and for Uba with multiple depletions also in winter months until April (analysis not shown in this manuscript). Thus, snowmelt is not important in these catchments for seasonal summer discharge, although it may be of high importance for spring discharge which is beyond the focus of this study. The reverse line of argument can be applied for the relatively high 5 importance of snow coverage for the high altitude central Tien Shan catchments (10. Naryn & 11. Upper Naryn) with mean annual temperatures below zero (cf . Table 1), where snow coverage alone explains up to almost 40% of the explained variance by the MLR models, in addition to the share of snow coverage contained in the composites. For these catchments snow coverage alone is thus already a good indicator of the seasonal discharge.
• In terms of predictor importance no obvious differences can be detected on average (right panel in Figure 8) between 10 the Tien Shan and Pamir discharge regimes identified in the cluster analysis (Figure 3), with the exception of the Naryn catchments as stated above. The mean predictor importance figures for those catchments are all very similar.
This can indicate that although the variability of seasonal discharges varies with geographical location, the runoff generating processes seem to be similar.

15
This general interpretation of the predictor importance shows that the selection of the predictors, particularly the change of predictors with prediction months and geographic region, has some hydrological meaning. Due to the simplicity of the approach and the simple linear relationship between the predictors, it is unlikely that more hydrological process information and understanding can be extracted from the MLR results. If this can be achieved at all, then on individual catchment basis only and by the interpretation of the exact predictors, i.e. not aggregated by classes as above. This is, however, beyond the 20 scope of this study. But nevertheless, the observation described above indicate that the general runoff generation processes can be described by linear models, and that the presented forecast results are unlikely obtained by pure chance only.

Potential of operational application
A lot of management and strategic decisions are based on seasonal forecasts of water availability in CA. The main consumer of water resources in the Aral Sea basin is the agricultural sector, which is based on one of the world's largest irrigation systems (Dukhovny and de Schutter, 2011). Very important decisions based on water availability forecasts are the planning of 10 agricultural production crop types and water allocation through the irrigation network. Also the estimation of agricultural yield is related to water availability and is needed for country income planning that heavily depends on agricultural export in some countries. Therefore reliable forecasts of seasonal water availability is essential for the economies of Central Asian states.
In order to design a generic and readiliy applicable forecats tool the presented method was designed according to the needs and data availability of the Central Asian Hydromet Services, which are responsible for the seasonal forecasts. The method is based on station data readiliy available to the state agencies, thus fulfilling a core prerequisite for an operational 5 implementation of the method. Moreover, the procedure for deriving forecast models is fairly simple and implemented in the open source software R. Therefore no limitations due to licence issues exist. The model development is automated requiring only some basic definitions as e.g. the formatting and provision of the predictor data as ASCII text files, and the specification of the prediction month. Therefore the code can be applied by the staff of the Hydromet Services after a short training. However, it has to be noted that the provided predictor data should be as complete as possible in order to avoid 10 spurios model fitting results (overfitting). Due to the automatic model fitting the algorithm may find best performing models fitted to a few years only, if too many predictor data are missing. The chances of overfitting are then greatly increased as the degree of freedom of the linear models, i.e. the ratio of the years used for fitting to the variables in the prediction models, decreases. The presented model system can also be run with alternative predictor data. For example, it has been tested using gridded 15 ERA-Interim re-analysis data for precipitation and temperature, averaged monthly over the individual catchment areas.
Similar, if not better results as presented were obtained. However, due to the latency of at least two months until the data is released, an operational use of the model system with ERA-Interim data is not feasible at the moment.

Discussion and Conclusions 20
The presented study aimed at the development of a flexible and generic forecast model system for the prediction of the seasonal (April-September) discharge in Central Asian river basins, with the final goal of operational use at the Central Asian Hydromet Services. In order to achieve this the data requirements were kept as low as possible, using only monthly precipitation and temperature data from a single station in the individual catchments, accompanied by operationally processed monthly MODIS snow coverage data and monthly antecedent discharge. Based on this core predictor data set, a variety of monthly, multi-25 monthly and composite predictors were automatically derived for different prediction dates. The predictors were then used for predicting the seasonal discharge with Multiple Linear Regression models (MLR). In order to avoid overfitting, restrictions were set on the selection and number of predictors in each MLR, and the models were tested for robustness by a Leave-One-Out Cross Validation (LOOCV). An ensemble of prediction models was then selected based on the best Predictive Residual Error Mean of Squares (PREMS) of the LOOCV. 30 The prediction model system was tested for the period 2000 -2015 on a selection of 13 different river basins in different geographic and climatic regions, and with different catchment characteristics. It could be shown that the models provided good to excellent predictions for all catchments and for all defined prediction dates, resp. lead times. For the first prediction on January 1 st , i.e. for a lead time of three months, the explained variance (expressed as adjusted R 2 ) is already high in the range of 0.46 -0.86 for 9 catchments. For the following prediction on February 1 st the explained variance is above 0.59 for all catchments, and increases further with the following months. For the important prediction date for the planning of water resources in the region on April 1 st just before the high flow season, adj. R 2 values of the best models for each catchment are 5 in the range 0.68 -0.97, indicating exceptional high performance for a seasonal forecast.
The automatic selection of the predictors and their importance revealed some geographic or temporal patterns. Geographically the northern Altai catchments differ in the predictor selection of the best LOOCV-MLR models from the other regions as snowmelt in this region has less contribution to seasonal discharge (April -September), with snow cover often reduced to zero already in early spring. For all catchments the importance of antecedent discharge is increasing with progressing prediction 10 dates. This is plausible from a hydrological perspective: While during the winter months the discharge is dominated by groundwater contribution, the discharge in April and later contains information about the snow melt process, and thus predictive power. Moreover, for predictions following April 1 st the antecedent discharge represents already part of the predictand, and has thus an even higher predictive power. This means in summary that the selected predictors and their importance have some hydrological meaning, thus supporting the validity of the forecast models derived by the model system. 15 However, it has to be noted that specific features of runoff generation in the catchments cannot be detected and discovered by the rather abstract level of predictors, predictor importance and the very basic catchment characteristics.
Overall, the presented simple forecast system proved to be able to provide robust, very skilful, and reliable forecast models for Central Asia. Moreover, it also provides a generic and flexible tool for the development of seasonal discharge forecast models for Central Asian rivers. This tool can be used by the responsible Hydromet Services without the need for larger investments 20 in hardware, software, and education and training of staff. In fact, the model system is already tested in four Central Asian national Hydromet Services. The forecasts provided by the MLR models for the summer discharge of 2017 is benchmarked against existing forecast routines and finally the measured discharge in late fall 2017.
The reason for the high performance is surely the separation of the largest share of annual precipitation (snow in winter), and the runoff generation (snow melt in spring and summer). Due to this temporal separation there is no need to perform a seasonal 25 forecast of the precipitation for the summer period, which is very difficult and uncertain in Central Asia (Gerlitz et al., 2016).
The forecast is rather based on an estimation of the snow pack accumulated in winter and its snow water equivalent, for which the predictors and their combinations provide proxy information. Moreover, the proxy information is not forecasted, but measured, thus providing more reliable information compared to forecasted predictors. As the timely separation of precipitation and runoff is a unifying feature of all Central Asian headwater catchments encompassing high-mountain ranges, 30 the model system is able to perform exceptionally well for all tested catchments, though with different predictor combinations.
It is thus also very likely, that the model system will also work well in the Central Asian catchments not included in this study, with some limitations for very small catchments. Moreover, it can be reasoned that it is likely that the model system will also work well in other regions with similar climatic settings, e.g. the South American dry Andes or the Western U.S. (e.g. the Sierra Nevada). The provided information of seasonal water availability could also be used in dam operation and dam safety procedures, and strategic flood hazard management plans.
Further research using the same method could focus on the use of near-real time satellite based data for the forecasts. As indicated by the successful experiment using ERA-Interim data as predictors, spatially aggregated temperature and precipitation data for whole catchment areas could further improve the forecasts. Additionally the method could be extended 5 to ungauged basins. A potential alternative source for precipitation data could be the near-real time TRMM rainfall estimates.
This needs to be accompanied by a satellite based temperature product. MODIS based surface temperature estimates could serve as data source for this. Another completely different data source with prediction potential are the water storage variations derived from gravity records of the GRACE mission. Ongoing research aims at the provision of a near-real time product of GRACE (EGSIEM project, egsiem.eu), which could then be used for operational forecasts. Considering that the gravity based 10 water storage variation should actually map the overall winter accumulation, it can be expected that this data could well serve as predictor for the seasonal discharge. However, due to the low spatial resolution, GRACE data are currently applicable for larger river basins only. 10.1029/2006wr005653, 2007.

Annex
Annex 1: Correlation of seasonal discharge to sub-seasonal discharge Figure A1: Comparison of seasonal discharge for the whole vegetation period April to September to sub-seasonal discharge time series taking the Naryn basin as example. The sub-seasonal series are highly correlated to the seasonal time series. Numbers in the legend provide the linear correlation coefficient of the sub-seasonal discharges to the seasonal discharge of the whole vegetation period.

Annex 2: Predictors used for the different prediction dates
The following paragraphs list the predictors created and used for the different forecasts dates, ranging from January 1 st to June 1 st . The predictors are abbreviated, with snowcov and sc denoting the snow coverage in the catchment derived by the MODSNOW-tool, precip the station records of precipitation, temp the station records of temperature, Q the discharge recorded at the river gauges. Catchment characteristics and the locations of the gauges are listed in Multi-monthly values are mean values of the monthly values spanning over several months, whereas the range of the months included is indicated by the concatenation of the indicators of the months, e.g. janapr means multi-monthly means for the period January to April, or febmar indicates the mean of the months February and March. The predictor abbreviations are 20 combined with the indicators for the months. snowcov_apr thus stands for the mean snow coverage of the catchment in April, or precip_janmar for the mean of the monthly precipitation sums for the months January to March.
For the composites the predictors included are listed by their abbreviations, followed by the indicators for the months. For calculating the composites, the monthly values of the predictors denoted by the month indicators are multiplied. E.g.
sc_temp_mar thus means the product of the mean snow cover in March and the mean temperature in March, or 5 sc_temp_precip_janmay denotes the product of the multi-monthly means January to May of snow coverage, temperature and precipitation.  The residuals of the models are tested for normality by the Shapiro-Wilk test for normality. Doing so, one has to bear in mind that this test is based on a sample size of maximal 16 values for each model only, so the test may not provide meaningful results. The table below shows the test result for every model, catchment, and forecast month. A "1" indicates normal 5 distributed residuals, "0" not normal distributed residuals. "NA" indicates that no more models with significant predictors could be found. For every forecast month up to 20 indices are given according to the set of best 20 models to be retained.

Predictors used for prediction on January 1 st
The table shows that for most of the models (89.5%) the test was positive, i.e. the residuals are normally distributed, even for this rather low and possibly not representative sample size.
Furthermore it was tested if the residuals are independent applying a test for autocorrelation with lag 1 at significance level p = 0.05. In the table below a "0" indicates independence, a "1" dependence. It shows that 95.8% of the models have independent residuals.
Last the Breusch-Pagan test for heteroscedasticity was applied to the residuals. This test shows that 99.5% of the models have homoscedastic residuals. In the table below a "1" indicates homoscedastic residuals, a "0" heteroscedastic residuals according to the test.