Hydrological modeling as a tool for water resources management of the data-scarce Brahmaputra basin

Expanding over four different countries, the Brahmaputra River basin occupies an area of 5,42,450 km² with the outlet (23.576N: 89.442E) in Bangladesh. The Brahmaputra is one of the most highly sediment-laden rivers of the world (Goswami 1985). Also, this river basin is located in an area of high structural instability as a large number of earthquakes have been evidenced in the Himalayan catchment through which it flows (Lahiri & Sinha 2012). The Brahmaputra carries the largest amount of water and silt of all Indian rivers. However, this is probably the least exploited basin, even though it has enormous potential regarding water resources management and hydropower generation. Hydrologic models are among the available tools used to estimate the parameters required for water resource planning and management.

There have been many case-specific studies on definite parts of the Brahmaputra basin (Goswami 1985; Sarma 2005; Akhtar et al. 2011; Ghosh & Dutta 2012; Gogoi et al. 2012; Apurv et al. 2015; Sahoo & Sreeja 2015). These studies addressed the local issues only, rather than the whole basin. On the other hand, only a few studies have so far been carried out for this important river basin regarding bank line migration (Singh et al. 2004), ion chemistry (Hren et al. 2007), and applicability of Lacey's regime equations (Swamee et al. 2008). None of these studies addressed the hydrologic assessment of the basin. However, Rao et al. (2009) attempted a hydrologic assessment, but the entire Brahmaputra basin was not covered in their study. In another study, although the whole basin was considered by Aktar et al. (2015), they failed to establish whether the hydrologic model developed would be applicable for deriving outputs across the basin. Here, the model was calibrated at only one location and no consideration was paid to spatial outputs. As such, we have attempted in this study to understand the hydrology of the entire Brahmaputra River basin so as to facilitate the water resource planners with science-based information at the desired locations.

Understanding the Brahmaputra basin has always been challenging due to its complex characteristics, both hydraulically and hydrologically (Rao et al. 2009). Although hydrologic models like SWAT, VIC (Variable Infiltration Capacity), WMS (Watershed Modelling System), HEC-HMS (Hydrologic Engineering Center–Hydrologic Modeling System) serve the purpose of hydrological analysis, they require correct input information. Even ANN (artificial neural network) can effectively be used to model rainfall–runoff processes (Hasanpour et al. 2014). The Brahmaputra basin is characterized by high spatial variations of topography, land use, soil properties, and weather components. As such, the challenge to model this basin lies in hydro-climatic information across and beyond the national boundary. This is because the data scarcity situations widely impact the hydrologic model results, as evidenced by many studies (Yu et al. 2011; Ploeg et al. 2012; Tabari et al. 2012; Valdivieso & Sendra 2014). Here, the spatial scarcity (Brath et al. 2004; Chen et al. 2018) as well as the temporal scarcity (Chen et al. 2018) of data are the major concerns for obtaining poorer simulation results. In the absence of available data, however, the missing variables may be generated (Nyeko 2015) and/or explicitly trained as constraints in mathematical programming (Pande et al. 2012). Unfortunately, these approaches often provide less satisfactory simulation results. In this study, we, therefore, have tried to incorporate the presently available data at many sources locally or globally.

The weather variables are the main components to widely impact the model results, but sufficient observed data are unavailable to model this river basin. Indeed, the observed weather records in large parts of the world are often scarce, discontinuous, and frequently contain discrepancies (Hughes 2006; Koutsouris et al. 2016). Even if longer records are available in the gauge observations, they do not afford a reliable spatial depiction (Yatagai et al. 2009) and, as such, evaluations, especially for rainfall characteristics, are important prior to applications in water resource schemes (Nyatuame et al. 2014). In this study, evaluations of various datasets are carried out based on the simulation results of the SWAT model, set up independently for different input datasets. Subsequently, a proper dataset suitable for the data-scarce Brahmaputra basin is identified.

The process of calibration and validation of a hydrologic model is not possible at the ungauged sections where observed data are not available. Also, the calibration of the large Brahmaputra basin at a single site (Aktar et al. 2015) is not judicious due to its heterogeneous characteristics. However, the present study adopted a multi-site calibration technique for addressing better intra-watershed simulations.

The flow production in the present basin is estimated by implementing the SWAT model, which has previously been successfully applied by many researchers across the globe. Some such studies include flow and sediment simulation (Rosenberg et al. 1999; Srinivasan et al. 2005; Zhang et al. 2008; Valdivieso & Sendra 2014; Daggupati et al. 2015; Dahal et al. 2016), as well as quantitative and qualitative assessment of water quality (Fohrer et al. 2001; Grizzetti et al. 2003; Bekiaris et al. 2005; Abbaspour et al. 2007; Guse et al. 2015; Chen et al. 2018). Another reason for using SWAT is that it is computationally efficient and capable of continuous simulation over long periods of time (Borah & Bera 2004; Arnold et al. 2012). SWAT is a continuous, semi-distributed process-based river model, and describes interactions between the elementary units of various hydrological processes (Lempert & Ostrowski 2002). As the model developers need to make many critical decisions during development and calibration (Zhang et al. 2009; Arnold et al. 2015), the present study is done by utilizing various combinations of datasets in order to ensure accurate, real-world situations.

The Brahmaputra River basin has huge potential with regard to water resources projects, although severe hazards like flood and bank erosion that result in huge loss of life and property every year are also much pronounced. Developing a comprehensive management policy is the current priority, but requires adequate data regarding flow, sediment, soil moisture, water quality, etc. We, therefore, adopted hydrologic models as a tool in order to estimate these data for the Brahmaputra basin, although their establishment was a major challenge due to the lack of input data. Furthermore, we could not find any observational records for the majority of the tributaries and sub-tributaries that join into the main stem of the Brahmaputra River, within Indian territory in particular. Therefore, modeling these tributary basins individually is really a challenging task, which can probably only be solved through modeling practice at a comprehensive level covering the entire transboundary Brahmaputra River basin, for which certain data are at least available. The present study presents a modeling-based approach to understand the dominant processes controlling the water balance to bridge the gap of base-line knowledge of water resources in the data-scarce Brahmaputra basin. The main aim of this modeling exercise is to identify a suitable dataset for carrying out hydrologic assessment across the basin, including locations having no historical records. This study would provide information required for management of water resources across the basin including its sub-basins.

The proposed study is intended for the transboundary Brahmaputra River basin (Figure 1). The Brahmaputra River originates in southern Tibet at an elevation of 5,300 m. Out of its total length of 2,880 km, the Brahmaputra covers a major part of its journey in Tibet as ‘Tsangpo’ River (also popularly known as ‘Yarlung Zangbo’ in China) flowing 1,625 km in Tibet parallel to the main range of the Himalayas before entering India through Arunachal Pradesh where it is popularly known as the ‘Siang’ River. It flows for about 35 km and is joined by two other major tributaries, Dibang and Lohit, at a place in the west of Sadiya, Assam, India. From this confluence point, the river is known as the Brahmaputra until it crosses Assam, traveling 720 km and entering Bangladesh as Jamuna. After traveling a distance of around 95 km from this location, the Brahmaputra is joined by Surma, another major river of Bangladesh. From this point, the Brahmaputra is popularly pronounced as Meghna which then finally merges into the Bay of Bengal to end its journey. Although the main river flows through China, India, and Bangladesh, its basin boundary encompasses Bhutan as the fourth nation (Figure 2(a)).

In this approach, the results of the SWAT model obtained by using various weather datasets are presented. Gassman et al. (2007) reported on SWAT as a physically based model that simulates the physical processes through input parameters like topography, land use, climate variables, and soil properties. In this study, the basin boundary is delineated through the use of Shuttle Radar Transmission Mission (SRTM) digital elevation model (DEM) of 90 m spatial resolution. In order to form a single DEM, mosaic in ArcGIS 10.4 was carried out for 20 DEMs so as to cover the entire basin. The stream network (Figure 2(b)) is generated for the threshold drainage area of 20,000 km².

In SWAT, the Brahmaputra basin watershed is divided into 41 sub-watersheds, which are then further subdivided into 1,578 hydrologic response units (HRUs). The HRUs represent an area consisting of dominant land use, soil characteristics, topography, and management practices. The present study utilizes 0.5 km resolution MODIS-based global data for the land use and to capture the land cover of the present basin. A soil map having 0.9 km spatial resolution, and provided by the Food and Agriculture Organization (FAO) has been used in the present study. Based on this soil map, we found as many as 39 soil types in the entire basin. The Brahmaputra basin is characterized by a very high variation in topography (Figure 3(a)).

Water balance is the driving force behind all the processes in SWAT because it impacts plant growth and the movement of sediments, nutrients, pesticides, and pathogens (Arnold et al. 2012). The hydrologic cycle depends upon weather variables such as rainfall, air temperature, solar radiation, wind speed, and relative humidity, which ultimately control the water balance of a watershed. The present study uses four sets of climate data obtained from various sources, as shown in Figure 3(b). The first set of weather variables is the Climate Forecast System Reanalysis (CFSR) data, provided by Texas A&M University (TAMU) as gridded global, high resolution coupled atmosphere–ocean–land surface–sea ice system weather data. These weather data consist of six variables: precipitation (pcp), maximum/minimum temperature (tmp), relative humidity (hmd), solar (slr), and wind (wnd), and are available for 36 years (1979–2014). Out of numerous grid points available within (21°–33°) North latitude and (82°–98°) East longitude, only 35 stations are selected so that at least one station falls on each of the major sub-basins, based on its spatial extent. The second weather datasets include observed data at six gauge stations in Indian territory, and collected from the Indian Meteorological Department (IMD). These datasets consist of four variables (pcp, max/min tmp, wnd), and were collected for a period of 25 years (1991–2015). The third dataset includes the measured weather data at a few stations in the Tibet region of China (see Results and discussion section), and contains all the six variables spanning records of 22 years (1991–2012). Finally, the IMD gridded weather data generated by Dr. Balaji Narasimhan, Associate Professor, IIT Madras have been used as the fourth weather dataset. These data contain only three variables (pcp, tmp min/max) and are available up to the year 2005. Any variables if not available in the later three datasets are copied from the nearest corresponding CFSR dataset. The missing values are filled with negative 99, which means SWAT is directed to generate the missing values based on the available records.

The ArcSWAT ArcGIS extension is a graphical user interface for the SWAT, which is a semi-distributed, time-continuous watershed simulator operating on a daily time step (Arnold et al. 1998). Here, a total of six SWAT models have been developed by using the same DEM, land use map, soil map, and slope classifications in ArcSWAT2012. However, different weather datasets were utilized for each model. The first SWAT model is developed through the use of only the 35 global CFSR stations' data; whereas, the others are developed by using a combination of weather datasets, as indicated in Table 1.

The first four SWAT models were simulated for 22 years (1991–2012) on a monthly time step basis whereas the last two models (MODEL_5 and MODEL_6) were restricted to simulation for only 15 years (1991–2005), due to unavailability of the IMD gridded weather data beyond 2005. Therefore, the calibration and validation periods of these two SWAT models were chosen differently than the previous ones.

SWAT-CUP is an interface to SWAT, and any calibration/uncertainty or sensitivity program can easily be linked to the later models. It provides a decision-making framework that incorporates a semi-automated approach (SUFI2) using calibration and incorporating sensitivity and uncertainty analysis (Arnold et al. 2012).

The key parameters are necessary to be identified prior to calibration (Ma et al. 2000), and the present study identifies nine sensitive parameters (Table 2) during the iterations of each SWAT model. As the watershed models suffer from various uncertainties in regard to model, input, and parameters, the present study has adopted an uncertainty analysis for the parameters identified through the sensitivity analysis. Table 2 shows the range of sensitive parameters to mark the uncertainty level along with the best fit values at which calibration and validation results are obtained. Normally, three to five iterations are sufficient (Abbaspour et al. 2015) to conclude the best fit of the calibrated flows with respect to the measured ones. Out of four locations (Figure 2(b)) of available discharges, the outlets on the main stem, namely, (i) Bhomoraguri (26.606417N; 92.843865E), (ii) Pandughat (26.176872N; 91.687088E), and (iii) Pancharatna (26.201282N; 90.575692E) are used for model calibrations and validations. Golaghat (26.50278N; 93.95194E), on the Dhansiri tributary basin, is used for cross verification of the model results. The evaluations of model performance are carried out using certain statistical approaches: (a) the coefficient of determination (R²), which indicates the strength of the relationship between the observed and simulated flow values; (b) Nash–Sutcliffe coefficient (NS), which shows how well the observed values fit the simulated values; (c) P-factor, which represents the percentage of measured data bracketed by the 95% prediction uncertainty (95PPU); and (d) R-factor, which is the measure of uncertainty.

Initially, all the models were simulated, calibrated, and validated at a single-outlet basis, for discharges (collected from Central Water Commission, (CWC)) at the three main stem outlets individually. Based on these results, MODEL_6 was identified as the suitable model for the present condition and thus forwarded for the spatial calibration/validation. This is because the performance of the default model should not be too drastically different from the measurement. If so, often calibration can be of little help (Abbaspour et al. 2015). Finally, the model calibrated and validated at the main stem outlets is cross-verified for outputs at the tributary basin to check whether the same would produce satisfactory results at the sub-basin level or not.

The Brahmaputra is one of the largest river systems in the world and, as such, researchers have faced a huge challenge to understand the physical processes (Pahuja & Goswami 2006). It is, therefore, the need of the moment to develop a comprehensive water management cooperation mechanism for water sharing between the co-basin nations – India, China, Bhutan, and Bangladesh (Sharma 2014). In one approach, Scholars from Bangladesh, China & India (2014) articulated some issues about water resources of the Himalayan region and were of the opinion that these could be brought before the rest of the world to forward constructive discussions and information exchange by international experts, media, and civil society. Thus, sharing data between the riparian countries can only be achieved through an inter-governmental understanding (Thu & When 2016); but, there is no such agreement known to us between the co-nations. Therefore, obtaining flow data of the river from the upper riparian country was not possible. However, for modeling purposes, we have used some precipitation data along with a few other weather components of China (Tibet), obtained through academic collaboration. This section of the paper highlights the results of the present modeling study along with the subsequent analyses.

It is very important to use adequately correct data while establishing a hydrologic model, without which reasonably valid outputs can rarely be achieved. Figure 4 shows a scatter plot between the simulated and observed flows for the models at Pancharatna (Figure 2(b)), as an example. Similar trends were also observed at the other two outlets, namely, Bhomoraguri and Pandughat. This plot provides an idea of how the simulated flow is different from the observed flow. A perfect match between the two represents a straight 1:1 line.

It is evident from Figure 4 that the MODEL_1 which uses only the CFSR global stations' datasets exhibits the least correlation (R² = 0.607) value, indicating the highest disagreement between the measured and simulated flows out. The flows out of MODEL_2 and MODEL_3 that use observed stations' data of IMD and China, respectively, also did not show a good correlation with the observed flow data. Here, although there are improvements in model outputs, significantly wide changes have not been observed, due to non-capture of all the observed weather datasets by the second and third SWAT models. This is because only the weather station close to the geometric center of sub-watersheds is picked up by SWAT. On combining all the observed weather datasets with the global CFSR dataset as used in MODEL_4, however, the correlation further improved with R² value being 0.65.

The simulation result of MODEL_5 has surprisingly improved as the IMD gridded data are added to the global CFSR weather data. Furthermore, MODEL_6 provided the best correlation (R² = 0.816) value among all models developed in this study. This analysis indicates MODEL_6, which uses the global CFSR dataset in conjunction with the IMD gridded, as well as the observed data over China, would provide the best performances in simulating the discharges, and it, therefore, may be taken up as the default model for the present basin.

The calibration and validation of the SWAT models are done in the SWAT-CUP platform that tries to capture most of the measured data within 95% prediction uncertainty (95PPU) during iteration.

The two indices P-factor and R-factor are used for the goodness of fit between the simulated and the observed values. The P-factor is the fraction of observed data bracketed by the 95PPU band and varies from 0 to 1, where 1 indicates 100% bracketing by 95PPU. On the other hand, R-factor is the thickness of the model envelope (i.e., 95PPU) and should be around 1, but a value of <1.5 is acceptable (Abbaspour et al. 2007) for this index.

Initially, all the SWAT models are calibrated and validated on a single-outlet basis. In the present analysis, we show the calibrated and validated results at Pancharatna outlet only, as a continuous graph (Figure 5) for simplicity and clarity of presentation. Here, we show the time series plot for ‘non-dimensional flow’ rather than the actual values due to the data secrecy policy of the concerned authority. The first four models, i.e., MODELs_1, 2, 3, and 4, were calibrated for ten years (1999–2008) and validated for four years (2009–2012) by skipping eight years (1991–1998) for initiation of the hydrologic response. Meanwhile, the last two models, i.e., MODELs_5 and 6, were calibrated for nine years (1993–2001) and validated for four years (2002–2005). Although the more the skipping periods the better the initial response, we could not make it longer for the last two models due to limitations of model input data, and a two-year (1991–1992) period only was provided for them.

It is understood from the results that MODEL_1 provides the least satisfactory output, as evidenced by the calibration statistics in Table 3. Here, the value of R² (0.78) is fair, but the other values like NS (0.68), P-factor (0.60), and R-factor (0.63) are relatively inferior for acceptance of the model. On checking the validation statistics, however, slightly better performances were noticed. The statistics of MODEL_2 and MODEL_3 are found to comply satisfactorily. However, MODEL_4 provided the best statistics among the first four models, during both the calibration and validation. Surprisingly, MODEL_5 which showed good agreement during the simulation (Figure 4) failed to produce good statistics during calibration and validation, although they were carried out by using the same sensitive parameters (Table 2) as the previous ones. Among all models, MODEL_6 is found to perform the best, since it provides a very good strength of statistics for all the parameters, even during both the calibration and validation (Table 3). This fair result is due to using a proper combination of input datasets in this model (MODEL_6), as compared to those used in the other models (MODELs_1–5). The R-factor (0.97) value during calibration (Table 3) is nearly equal to the desired value (1.00). It indicates a very good measure of parameter uncertainty. As such, this model has only been extended for subsequent analyses.

The calibration at a single site is judicious for the small watershed only, because of homogeneous watershed characteristics. On the other hand, a large watershed needs a multi-site calibration in order to represent heterogeneous characteristics. The single-site calibration if done on a large watershed may result in a combination of under/overestimation for values (Qi & Grunwald 2005), due to consideration of homogeneous characteristics over the entire basin. To have a better intra-watershed spatial accuracy, the present study includes multi-site calibration/validation at the three main-stem outlets simultaneously. This spatial accuracy is done only for the best SWAT model as identified earlier (sections Simulation results and Single-outlet calibration/validation), i.e., MODEL_6, assuming that the other models would not lead to better performances, even during the multi-site calibration and validation.

The time series plots are shown in Figure 6 as a continuous graph for simplicity in the presentation. Due to the non-symmetry of available discharge data length, different calibration periods were chosen for the outlets (Figure 2(b)). However, the same validation periods for all the three sites were chosen (2002–2005).

Due to its complex characteristics, a large watershed may not be expected to produce good statistics against each outlet during multi-site calibration; however, a balance among the statistical parameters is required for reasonable acceptance of a model. During this spatial analysis (Figure 6(c)), the most downstream outlet (i.e., Pancharatna) is found to conform with good statistics during calibration (R² = 0.84; NS = 0.74) and validation (P-factor = 0.94; R-factor = 0.91). Similarly, the model is found to produce satisfactory statistics at the other two outlets, namely, Bhomoraguri (Figure 6(a)) and Pandughat (Figure 6(b)). However, the calibration results at Pancharatna, especially the value of R² (0.84) are relatively low as compared to Bhomoraguri (0.86) and Pandughat (0.91). This is probably due to different lengths of warm-up period maintained for initiation of watershed responses by the hydrologic models. The warm-up period for Pancharatna is only two years (1991–1992) against eight years (1991–1998) for the later two outlets. Overall, the spatial calibration/validation results of MODEL_6 may be termed as satisfactory.

The present SWAT models were developed for the large Brahmaputra River basin, but whether this model built for a basin area of more than 5,42,000 km² would be suitable for its sub-basins occupying smaller areas, or not, needs to be evaluated. This is done through cross verification of the model results. This analysis is possible for several outputs (discharge, sediment, evapotranspiration, groundwater substances, etc.) of the SWAT model, provided historical records are available. Unfortunately, the majority of the tributaries have hardly any observational records, and we could collect only one set of discharge data for a tributary. Here, we compared the final MODEL_6 results pertaining to Golaghat, an outlet at Dhansiri tributary (Figure 2(b)), with the measured discharges. It is evident from the plot (Figure 7) that the model simulated values follow a fair correlation (R² = 0.77) with the measured discharges, in terms of magnitudes and pattern. As such, the present model is expected to provide satisfactory outputs at the other tributaries too. Unfortunately, the model verification was not possible at any other tributary since no historical records are available to us. It is evident from the analyses that the SWAT model could provide satisfactory outputs at both the main stem outlets as well as tributary outlets of the data-scarce Brahmaputra basin.

The SWAT model provides location-specific hydrologic information across a river basin which may be utilized by the water resource managers for planning and designing hydraulic structures. This information belongs to the quantitative and qualitative assessment of river flows. A list of certain parameters is enumerated in Table 4 corresponding to a location (latitude: 29.15485N; longitude: 95.006979E) termed the ‘Indo-China Border’ (Figure 2(b)). Here, monthly values are obtained from the outputs of the final run of MODEL_6 which are then averaged during 1991–2005. The SWAT hydrologic model output for river flows at this location varies widely from a very low value (28.0 m³/s) to a very high value (8,364.8 m³/s). Similarly, a broad range of outputs has also been noticed for the sediment concentration with a maximum value of 445.815 mg/L during the month of June. High values of sediment concentration may lead to loss of aquatic habitats, wetlands, and recreation attributes. Moreover, it becomes a concern for human health and erosion hazards. Water quality parameters like nitrogen, phosphorus, nitrates, ammonium, dissolved oxygen, etc. may be a concern for human and plant health. Based on the qualitative analysis for those parameters, the processes and cost of treatment of water, if derived for water supply projects, can be ascertained. Although the present model is capable of providing quite a lot of information regarding water quantity and quality at all desired locations across the basin, the present analysis is shown corresponding to one outlet only, as an example. Such information provided by the hydrologic models at a desired location may be utilized by water resource managers. As such, hydrologic models serve the purpose of deriving information at locations with no historical records, especially at those locations beyond the national boundary for such a transboundary river basin.

There has been a huge challenge for researchers to model the transboundary Brahmaputra River basin due to its complex characteristics and data-scarce situations. However, we can succeed in modeling this important basin, as evidenced by the results discussed above. Here, the situations arising from the limited information and reliability of the available data are handled by developing various SWAT models, in order to decide the proper dataset suitable for the Brahmaputra basin. The conclusions from this study can be drawn as follows. (1) The CFSR weather data provides an underestimation of flow values in the Brahmaputra basin. (2) The observed weather stations data (IMD gridded + China weather) when used in conjunction with the CFSR data provided near-real hydrologic assessment for the basin. This dataset that produces acceptable model-generated stream flow series in the case of both single and multi-outlet validation is considered as the best combination. Therefore, this dataset along with the model parameters can now be used with confidence, for estimating outflow at any desired location of this large data-scarce basin. (3) The SWAT hydrologic model developed for the large Brahmaputra River basin can provide satisfactory monthly outputs for discharges, even at its tributaries. Calibration and validation may be possible for many variables like sediment and water quality parameters, and corresponding to more locations distributed throughout the basin, and even for the daily time series. However, that was not possible in this study because of the non-availability of records to us. This may be extended as scope for future study. (4) For regions lacking historical records, as in the present case, hydrologic models can serve the purpose of providing information regarding flow, sediment, and water quality across the basin. This would enable the water resource planners and managers to decide their plans and actions for extracting the utmost benefits from the basin. (5) The suitable hydro-meteorological dataset identified in this study can be used for GCM downscaling and hence for climate change projection of the Brahmaputra basin.

We acknowledge the support of the CWC, Govt. of India, for providing necessary discharge data to establish model. We also acknowledge the Indian Meteorological Department for providing meteorological data.