Impact of climate change on the stream flow of lower Brahmaputra : trends in high and low flows based on discharge weighted ensemble modelling

The paper addresses the issue of climate change impact in Lower Brahmaputra river basin (LBRB). It uses global hydrologic model PCR-GLOBWB under 12 different GCM and ensemble discharge predicted from it to study the trend in the high and low flows in LBRB. Undoubtedly Brahmaputra basin is one the hot-beds to study climate change impact as is clear from many studies and physical evident. Since climate models often have ambiguous results, there is a strong need of ensemble based forecast. Hence the paper tries to addresses a very relevant and important issue which seems to fall pretty well in the scope of the current journal.


Introduction
Climate change is likely to lead to an intensification of the global hydrological cycle and to have a major impact on regional water resources (Arnell, 1999).The IPCC Fourth Assessment Report mentions with high likelihood that observed and projected increases in temperature, sea level rise and precipitation variability are the main causes for reported and projected impacts of climate change on water resources, resulting in an overall net negative impact on water availability and the health of freshwater ecosystems (Kundzewicz et al., 2007).
Among the river systems, the hydrological impact of climate change on Ganges-Brahmaputra Basin is expected to be particularly strong.There are three major reasons for this.First, stream flow is strongly influenced by the melt of snow and ice in the upstream part of the catchment.As 60 percent of the basin area has an elevation of over 2000 m cryospheric processes are deemed important when considering basin hydrology.Projected rise in temperature will lead to increased glacial and snow melt, which could lead to increased summer flows in some river systems for a few decades, followed by a reduction in flow as the glaciers disappear and snowfall diminishes (Immerzeel, 2008).This is particularly true for the dry season when water availability is crucial for the irrigation systems.Immerzeel et al. (2010) stated that the Brahmaputra is most susceptible to reductions of flow, threatening the food security of an estimated 26 million people.Second, the Ganges-Brahmaputra basin is highly influenced by extreme monsoon rainfall and flooding (Mirza, 2002;Warrick et al., 1996).If climate change results in changes of both the intensity and reliability of the monsoon, it will affect both high and low flows leading to increased flooding but possibly also to increased variability of available water, both in space and time (Postel et al., 1996).The latter refers to the fact that discharging water during floods and wet seasons cannot be used during the low flow seasons unless large storage systems are in place (Oki and Kanae, 2006).Third, climate change induced sea level rise results coastal flooding and riverine flooding by causing back-water effect of the Ganges-Brahmaputra basin along the delta (Agrawala et al., 2005).
The objective of this study is to investigate trends in both high and low flow for the Lower Brahmaputra River that may arise as a result of climate change.Compared to previous assessments (Warrick et al., 1996;Mirza, 2002;Immerzeel, 2008;Immerzeel et al., 2010) we do not build a basin-specific hydrological model for this purpose.Instead, we use existing results of a global hydrological model that was forced by data from 12 global climate models (GCMs) (Sperna Weiland et al., 2010) in a weighted ensemble analysis.The noveltly in this approach lies in that GCM-weights are determined based on the proximity of the associated streamflow simulations to observed streamflow (see Sperna Weiland et al. (2011) for a first application of this method).This approach is an improvement of other methods that have previously been applied.Immerzeel (2008) for examples uses a multiple regression model to predict streamflow at Bahadurabad, but in this case the ensemble results of a physical based distributed hydrological model are matched to observed discharges and hydrological processes are likely to be captured more accurately in the results.Also, the method by which we construct transient stream flow time-series can be considered as novel.Based on the constructed time series of transient stream flow (for the years 1961-2100) we then project trends in low and high flow statistics for the A1B and A2 emission scenarios.
In the remaining part of the paper we first describe the methodology of constructing the transient future time series of river flow in detail.We then show and discuss the results related to the analysis of both low and high flow analysis and conclude the paper by reporting and discussing the major findings.

The lower Brahmaputra river basin
The Brahmaputra is a major transboundary river which originates in the glaciated areas of the Kailash range in Tibet (China) at an elevation of 5300 m above the sea level (m a.s.l.).The river has a length of 2900 km, drains an area of around 530 000 km 2 and traverses four different countries (% of total catchment area in brackets): China (50.5 %), India (33.6 %), Bangladesh (8.1 %) and Bhutan (7.8 %).Average discharge of the Brahmaputra is approximately 20 000 m 3 s −1 (Immerzeel, 2008).The climate of the basin is monsoon driven with a distinct wet season from June to September, which accounts for 60-70 % of the annual rainfall.Immerzeel (2008) categorized the Brahmaputra basin into three different physiographic zones: Tibetan Plateau (TP), Himalayan belt (HB), and the floodplain (FP).These zones respond differently to the anticipated climate change.The TP covers 44.4 % of the basin, with elevations of 3500 m and above, whereas, HB covers 28.6 % of the basin with elevations ranging from 100 m a.s.l. to 3500 m a.s.l.The area with an elevation of less than 100 m a.s.l. is considered as FP and comprises about 27 % of the entire basin.This study is focusing on river flow in the lower Brahmaputra River Basin which belongs to the FP (Fig. 1).In the lower Brahmaputra, average temperature in winter is about 17 • C and summer temperatures are on average as high as 27 • C. Total annual precipitation is about 2354 mm concentrated in the monsoon months June, July, August and September (JJAS).The major discharge measuring station of the lower Brahmaputra is in Bahadurabad (Bangladesh).This is the only station in the lower Brahmaputra for which long-term observed records are available through the Bagnladesh Water Development Board.The data are of high quality and used in major hydrological studies for flood forecasting and other planning purposes.Therefore, long-term observed records from this station will be used to weigh the global hydrological model outputs resulting from the different GCMs.

Creating an ensemble of discharge time series for the reference period
To investigate the impact of climate change on hydrology we have to rely on combinations of runs of climate models and hydrological models.When it comes to climate projections, there is no single best model but rather a pool of models or model components that must be interrogated (Knutti, 2008).Projected values of models are inherently uncertain, because a model can never fully describe the physical system and complete confirmation of model output through verification and validation is impossible (Oreskes et al., 1994;Parker, 2006).Therefore, a collection or ensemble of models is preferably used to characterize the uncertainty in projections, while the credibility of projected trends increases when multiple models point in the same direction.Moreover, the average of a multi-model ensemble often outperforms single models when compared with observations (Gleckler et al., 2008;Reichler and Kim, 2008;Knutti, 2008).This study considers multiple outputs of 12 Global circulation models (GCMs).The output of these GCMs was used to force the global hydrological model PCR-GLOBWB.PCR-GLOBWB (van Beek and Bierkens, 2009;Bierkens and van Beek, 2009) calculates for each grid cell (0.5 • × 0.5 • globally) and for each time step (daily) the water storage in two vertically stacked soil layers and an underlying groundwater layer, as well as the water exchange between the layers and between the top layer and the atmosphere (rainfall, evaporation and snow melt).The model also calculates canopy interception and snow storage.Sub-grid variability is taken into account by considering separately tall and short vegetation, open water, different soil types and the area fraction of saturated soil and the frequency distribution of groundwater depth based on the surface elevations of the 1 × 1 km Hy-dro1k data set.Fluxes between the lower soil reservoir and the groundwater reservoir are mostly downward, except for areas with shallow groundwater tables, where fluxes from the groundwater reservoir to the soil reservoirs are possible (i.e., capillary rise) during periods of low soil moisture content.The total specific runoff of a cell consists of saturation excess surface runoff, melt water that does not infiltrate, runoff from the second soil reservoir (interflow) and groundwater runoff (baseflow) from the lowest reservoir.To calculate river discharge, specific runoff is accumulated along the drainage network by means of kinematic wave routing including storage effects and evaporative losses from lakes, reservoirs and wetlands.
In a previous study (Sperna Weiland et al., 2010, 2011) the output of 12 GCMs (Fig. 2 for short names) was used as input to PCR-GLOBWB.Daily precipitation and data to calculate daily reference potential evaporation were collected from the data portal of the Program for Climate Model Diagnosis and Intercomparison (PCMDI), (https://esg.llnl.gov:8443/index.jsp).For each GCM model runs for two scenarios, A2 and A1B, were selected that represent the upper range of possible CO 2 emissions.GCM runs comprised the 20C3M control experiment  and the future scenarios A1B and A2 (2071-2100).When multiple ensemble runs were available for one model, the first run was selected.Although the data portal does not provide all required parameters for the Hadley centre climate models, HADGEM1 has been included for it is frequently used in climate change studies.HADGEM1 data has been retrieved from the CERAgateway, http://cera-www.dkrz.de.

Ensemble weighting based on observed discharge
Rather than statistically downscaling each of the GCMs based on local meteorological data we attached a weight to each of the GCM-PCR-GLOBWB simulated outputs based on a novel method, following Sperna Weiland et al. (2011).Instead of weighting based on similarity of observed GCMbased input (e.g.rainfall), weighting is based on similarity of observed discharge.Using the mean monthly value of observed and simulated discharge during the overlapping period, a weighting factor for each model is computed according to Eq. (1).
Where, w is the weighting factor, j is month number, i is model number, σ i the standard error of discharge observations (m 3 s −1 ), which was assumed to be 25 % of the observed value, y j is the observed average discharge for each month j , and z ij is the mean monthly discharge for model i and month j.The estimate in standard error in discharge obervations is conservative, so as not to unjustifiably discard GCMs.Recent work of Di Baldassarre and Montanari (2009) showed that the overall error in river discharge observations ranges between 4.2 % and 42.8 %, with average of 25.6 % at the 95 % confidence level.The resulting weighting factors for those models with a significant non-zero value are shown in Table 1.It shows that MICRO received the highest value, followed by GFDL, GISS, CCCMA, CGCM, BCCR, HADGEM, NCAR, and ECHAM.We apply constant weights for the entire time series and do not vary weights for different months or flow conditions for robustness of the method.We validate this assumption by comparing the flow duration curve of the observations with the modelled discharges for the period between 1973 and 1995 in Fig. 3.The results show a good match between modelled and observed discharges and therefore we use a single set of weight for all conditions.Using these weighting factors, the daily weighted ensemble discharge ) and variance z ) lated of 1961 to 1990 and 2071 to 2100 according to Eqs. ( 2) and (3). (3)

Construction of a daily transient time series from 1961 to 2100
The 12 GCMs as obtained from the PCMDI used in Sperna Weiland et al. ( 2010) only provide runs for time slices (e.g.1961-1990 and 2071-2100).There are transient runs for some of the GCMs (e.g. at CERA-gateway), but certainly not for all of them.Therefore, to simulate transient time series of discharge for the period 1961-2100, for each of the GCMs the following steps were taken: For each year between 1991 and 2070 a random year is selected either from the reference period or from the projected period.The probability of selecting a random year from the reference period or from the projected period for year i depends on how many years year i is separated from either the reference period or the projected period.For example the probability (P r ) that for the year 2000 a random year is selected from the reference period is 0.88 according to Eq. (4).
Using this approach a complete time series is constructed from 1991 to 2070, resulting in a full time series from 1961-2100.The full time series from 1961 to 2100 is used in the subsequent analysis of trends in high and low flows.Using this approach year to year variability is preserved in the constructed time-series.It should be noted however that, as we sample discharges directly, we may encounter welding problems between subsequent sampling years: jumps between 31 December and 1 January.We have validated our approach by artificially reconstructing a transient time series during the observational period.We constructed a time series from 1980 to 1989 by sampling from the time slices 1970-1979 and 1990-1999 similar to what is described above.We then compare the daily data of the simulated transient ime series with the actual observations during 1980 to 1989 and derive a number of statistics.Our analysis shows that the Pearson correlation coefficient is 0.85, the bias is −2.3 %, the root mean square error is 9323 m 3 s −1 and the Nash-Sutcliffe criterion for model efficiency equals 0.71.These numbers show that our approach is valid and that simulated discharge measure observed discharge well.It should be noted that this may be different for river basins where seasonality in discharge is less pronounced.

Extreme value analysis
The low-flow regime of a river can be analyzed in a variety of ways depending on the type of data availability and the type of output information required (Smakthin, 2001;Pyrce, 2004).Here we use the N-day minima approach.Traditionally, the annual minimum (AM) values have been used for low flow frequency analysis, as droughts particularly become an issue when they persist.We use a 7-day low flow frequency using a moving average for the A1B and A2 scenario from 1961 to 2100.To estimate trends in high flow frequencies we performed a traditional extreme value analysis based on yearly maxima for different time slices.

Trends in discharge
We use linear trend analysis similar to Gain et al. (2007Gain et al. ( , 2008)).Before analysing the trend of the complete data series, we compare the trend between modelled discharge and observed records during the overlapping period of 1973-1995.For the observed records, the trend was 195 m 3 s −1 yr −1 whereas this value of modelled data was 173 m 3 s −1 yr −1 .This result shows that modelled outcomes are consistent with observed trend.Table 2 presents the annual and monthly trends in discharge.From 1961-2100 trends are calculated by first calculating a trend parameter per GCM and then calculating the weighted mean trend and its variance using Eqs.( 2) and ( 3).From this it can be tested whether a trend is significant or not, using a two-sided t-test.Similarly, the goodness of fit coefficient, R 2 is first calculated for each GCM subsequently the weighted average over all models is calculated.This analysis was done on both yearly average discharge as well as on discharge per month.Table 2 shows that on annual basis there is a strong positive trend in stream flow that is mainly caused by a strong increase in monsoon discharge.During the dry seasons a modest increase is observed.The only negative trend is found in May, but the correlation is small and the trend non-significant.This is because the monsoon may have weakend at the onset of the monsoon season and strengthened during the later months (CCC, 2009).
Seasonal average flow for both A1B and A2 scenario of four time slices are compared in the box-whisker plots of Fig. 4. Box plots were obtained by first calculating cumulative frequency distributions per GCM and then constructing a weighted cumulative frequency distribution by weighting values belonging to the same quantile.The statistics in the box plots are thus based on the weighted cumulative  frequency distribution.Figure 4 shows that the strongest increase in both average and extreme discharge is predicted for the summer and autumn periods.It also shows that changes in discharge distributions are quite similar between scenarios, except for summer and autumn (i.e.monsoon) maximum flows, where the increase is more pronounced for the more extreme A2 scenario.It should however be noted that future spring and early summer discharge may be underestimated as the model does not take into account the increase of melt from glaciers in the upstream parts of the basin, which does play an important role in the Brahmaputra (Immerzeel et al., 2010).

Flow duration curves
The Lower Brahmaputra River Basin (LBRB) is characterized by water shortages in the dry season and water excess and flooding during the monsoon months.To further understand the projected change in range of river discharge, we constructed flow duration curves (Smakhtin, 2001).First for each GCM a flow duration curve was estimated for four 20year time slices.Next, for each time slice the weighted flow duration curve was calculated by weighting discharge for a given duration.Figures 5 and 6 provide the results for the A1B and A2 scenarios respectively.As can be seen, the Q90 and Q95 flows, commonly used as low flow indices (Pyrce, 2004), remain relatively constant for both scenarios, while the larger changes occur for the larger discharges, i.e.Q25 and higher.
Figure 7 shows a projected decrease in the likelihood of severe low flow events.This is because due to an increase in precipitation that outbalances the increase in evapotranspiration.The differences between the scenarios and time slices increase over time and the A1B scenario yields a stronger increase in low flows than the A2 scenario, which may be related to a less strong decrease in evapotranspiration due to a smaller projected temperature rise.To show the difference between the 12 models we provide Fig. 8 which shows for the A1B scenario the weighted distribution as a boxplot of yearly average 7-day low flow.Figure 8 shows that the there is a large variation in low flows between model runs but that all model runs show an increase in 7-day low flow.

High flows
The results of the high flow analysis are shown in Fig. 9.The graphs are constructed the same way as Fig. 7, but now based on yearly maxima.Figure 9 shows a very strong increase in annual peak flow, which may have severe impact for flooding in the LBRB.In this case the A2 scenario is the most extreme in line with the steep increase in monsoon precipitation.The 1:10 year discharge is projected to increase from 82 000 m 3 s −1 currently to 140 000 m 3 s −1 by 2100 and a peak flow that currently occurs every 10 yr will occur at least once every two years during the time slice 2080-2099.It is striking that for peak flows with larger return periods the strongest increase already occurs during the first 20 yr.This could most likely be attributed to sampling variability resulting from performing the extreme analysis on relatively short 20 yr time slices resulting in more than the expected number of randomly selected years from the 2071-2100 time slice.This could be corrected for by performing the analysis repeatedly for each model on multiple transients constructed by Eq. (4).

Conclusions and discussions
In this study we applied a new method to construct a daily discharge time series from using a discharge-weighted ensemble based on inputs from 12 GCMs to a global hydrological model.Weighted discharge time series were subsequently used to analyze future trends in average flow and extreme flow results show that climate change is likely to im- prove dry season conditions in the LBRB.For both scenarios (A1B and A2), for all models and for all time slices both average flow and extreme low flow is projected to increase in size.Low flow conditions may even be slightly underestimated as the accelerated glacial melt in the upstream parts of the catchment may, albeit temporarily, further enhance low flow.The A1B scenario projects the strongest increase in low flow.On the other hand, our analysis also shows a large increase in peak flow size and frequency.The impact for the already highly flood prone plains of Bangladesh may be devastating, in particular in combination with the projected sea level rise.The A2 scenario projects the strongest increase in high flow.
For the assessment of streamflow of Ganges-Brahmaputra basin, previous studies (Warrick et al., 1996;Mirza, 2002;Immerzeel, 2008;Immerzeel et al., 2010) applied basinspecific hydrological model.Through statistical downscaling of six GCMs and using multiple regression analysis, Immerzeel (2008) found a sharp increase in the occurrence of average and extreme discharge at downstream for A2 and B2 storylines.Mirza (2002) used climate change scenarios from four GCMs as input into hydrological models and result of the study demonstrates substantial increases in mean peak discharges in the rivers of Ganges-Brahmaputra basin.But in our study, we use existing results of a global hydrological model that was forced by data from 12 global climate models (GCMs) in a weighted ensemble analysis.Through a weighting method, we prioritize these GCMs based on their relative performances and our analysis shows that the observed discharges can be simulated well using results of 4 GCMs.
The results in this paper show that all GCMs point toward an increase in discharge of the lower Brahmaputra river.However, it should be noted that there is quite some uncertainty about the change in South-Asian Monsoon strength, and most climate models have difficulty simulating mean monsoon characteristics and associated inter-annual precipitation variation (Annamalai et al., 2007;Yang et al., 2008).Experiments with regional climate models even show contradictory results (e.g.Kumar et al., 2006vs. Ashfaq et al., 2009).However, given all the evidence, an increase in peak flow and flood frequency is likely and adaptive measures should be seriously considered.
In this paper we performed no model simulations of our own.Instead we made use of a repository of existing runs of a global hydrological model forced by a multi-model ensemble of climate data for both a reference period and 2071-2100 projections.Comparable weighting methods have been applied for GCM ensemble averaging of precipitation and temperature (see Giorgi and Mearns, 2002;Räisänen et al., 2010), but applying the approach to discharge is new.By weighting the simulated discharge with discharge observations a multi-model ensemble analysis of climate change effects could be made for a particular location, in this case the lower Brahmaputra at Bahadurabad station.Through this, a form of implicit downscaling is achieved that also takes account of inter-GCM uncertainty, because an ensemble of GCMs is used in reconstructing observed discharge.Moreover, the method, which allows for a very quick and cheap analysis of the effects of climate change plus uncertainty, is quite generic and can be used at other locations in the world with discharge observations.The method is applicable in any case where a hydrological model is forced with an ensemble of climate models and a sufficient long time series of observed discharges is available.The method can be easily improved to allow for the case that none of the models is doing a good job in reproducing discharge by adding bias-correction methods.
Ideally, the hydrological community could make a repository where the results of combinations of different GCMs and different global hydrological models are stored; both reference runs and projections for future time slices.Analyses by the method presented in this paper could then be done very quickly for any large river in the world, but now also taking the uncertainty about hydrological response into account.To have transient runs would be even better, but given that they are only available for a few GCMs at this time, transients could be constructed similar to our method for rivers with a strong seasonal signals as in our case.Alternatively, instead of interpolating discharge itself, one could also construct a transient of statistics by first estimating discharge statistics for each time slice and then interpolating changes of these statistics between time slices.

Fig. 1 .
Fig. 1.Overview of the Brahmaputra river basin (red polygon), the Brahmaputra river (blue line), the outlines of the lower Brahmaputra river basin (shaded white) and the Bahadurabad gauging station (red dot).

Fig. 3 .
Fig. 3. Comparison of flow duration curve between observed and discharge weighted modelled data for the period 1973-1995.

Fig. 4 .
Fig. 4. Box plots of stream flow for different seasons and for different time slices.Box plot represents the multi-model weighted variation over the season.

Fig. 5 .
Fig. 5. Flow duration curve for observed and multi-model weighted discharge of A1B scenario of four different time slices.

Fig. 6 .
Fig. 6.Flow duration curve for observed and multi-model weighted dischargeof A2 scenario of four different time slices.

Fig. 7 .
Fig. 7. 7-day low flow for different return periods for different scenario's and time slices as obtained from a weighted average of 12 model outputs.

Table 1 .
Computed weighing factors for the different model forcings.

Table 2 .
Trends in monthly annual stream flow from 1961 to 2100.Trends and R 2 are first calculated per GCM and subsequently the weighted average calculated.All trends are significant at the 95 % confidence level, except for the trend in May discharge for the A2 scenario.