Spatio-temporal assessment of annual water balance models for upper Ganga Basin

The upper Ganga Basin in Uttarakhand, India, has high hydropower potential and plays an important role in the development of the state economy. Thus, an accurate knowledge of annual water yield is of paramount importance to this region. This paper deals with use of contemporary water yield estimation models such as the distributed Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model and the Lumped Zhang model and their validation to identify the most suited one for water yield estimation in the upper Ganga Basin. In previous studies utilizing these models, water yield was estimated by considering a single value of some important model parameters for the entire basin, which in fact show distributed variation at a finer (pixel) scale. Therefore, in this study, pixel-level computations are performed to assess and ascertain the need for incorporating the spatial variation of such parameters in model applications. To validate the findings, the observed sub-basin discharge data are analyzed with the computed water yield for 4 decades, i.e., 1980, 1990, 2001 and 2015. The results obtained are in good agreement with the water yield obtained at the pixel scale.

1 Introduction

An accurate assessment of key ecosystem services (ES) such as water yield has gained focus in recent years in ES modeling, as fresh-water availability in a region is essential for agriculture, industry, human consumption, hydropower, etc. (Redhead et al., 2016). Hydrological ecosystem services generally include drinking water supply, power production, industrial use, irrigation and many more. The accurate estimation of water yield further facilitates in the identification of hotspots for storm-water harvesting in order to fulfill fresh-water demand in the region (Pathak et al., 2017). The hydrological ES are dependent on different factors, such as watershed characteristics (e.g., topography, land use and land cover – LULC, soil type and climatic condition. To incorporate these parameters into assessment and decision-making, there has been a proliferation of ecosystem-modeling tools and methods. Models for ES evaluation often focus on using globally available data, accepting large number of spatially explicit inputs producing spatially explicit output, and limiting the model structure to key biophysical processes involved in land use change (Guswa et al., 2014). The precise estimation of ES using these models is a complicated task owing to spatial variability and dependence of ES on various topographical and climatic factors. Further, the validation and uncertainty assessments in model outputs have proven to be key obstacles to the application of ES models. In the literature, studies focusing on comparison of different ES models have projected some light over the model output validation issues; however, a lack of studies highlighting the validation of these models for Indian river basins still exists. The benefits that can be derived from ES should be analyzed and quantified in a spatially explicit manner (Sánchez-Canales et al., 2012). The uncertainties involved in the determination of spatial and temporal distribution of the climatic variables, especially precipitation, constitute a major obstacle to the understanding of hydrological behavior at the catchment scales (Milly and Dunne, 2002).

The Integrated Valuation of ES and Tradeoffs (InVEST) model, developed by Natural Capital Project (Tallis et al., 2010), is a tool that provides a framework for planners and decision makers to assess trade-offs among ES and enables their comparison in various climate and land use change scenarios. The model includes a biophysical component, which facilitates the provision of fresh water or water yield from different parts of the landscape, and a valuation component, representing the benefits of water provisioning to people. The model works on simplified Budyko theory, which has a long history and still continues to receive attention in the hydrological literature (Budyko and Ronov, 1979; Zhang et al., 2001; Zhang et al., 2004; Ojha et al., 2008; Zhou et al., 2012; Donohue et al., 2012; Xu et al., 2013; Wang and Tang, 2014). The InVEST model applies a one-parameter formulation of the Budyko theory in a semi-distributed manner (Zhang et al., 2004). The model is capable of quantifying the water yield of a catchment under the influence of change in different drivers, viz. climate variables and catchment characteristics (e.g., land use change). Various studies have been carried out in the past demonstrating the application of the InVEST model to different river basins around the world. Sánchez-Canales et al. (2012) carried out a sensitivity analysis of three parameters, i.e., z (seasonal precipitation coefficient), precipitation (annual) and ET₀ (annual reference evapotranspiration), and using the InVEST model for a Mediterranean basin, they found precipitation to be the most sensitive parameter for the study region. Later, Terrado et al. (2014) applied the InVEST model for the heavily inhabited defined as Llobregat river basin. The model is applied for both extreme wet and dry conditions, and the role of climatic parameters is emphasized. Hoyer and Chang (2014) applied this model in the Tualatin and Yamhill basins of northwestern Oregon under a series of urbanization and climate-change scenarios. The results show that the climatic parameters have more sensitivity than other inputs for a water yield model. Hamel and Guswa (2015) applied the same water yield model for the Cape Fear catchment, North Carolina, and concluded that the precipitation is the most influencing parameter. Goyal and Khan (2017) employed the InVEST water yield model for the hilly catchment by considering two catchments, i.e., the Sutlej River catchment and Tungabhadra River catchment. The climate parameters, i.e., precipitation and ET₀, are observed to be the most influencing parameters for water yield in both the river basins. With the aforementioned studies, certain factors exist that limit the application of InVEST model such as the absence or inadequate comparison with observed data, the calibration of the model without prior identification of sensitive parameters and a lack of validation of the predictive capabilities in the context of land use and land cover change (Bai et al., 2012; Nelson et al., 2010; Su and Fu, 2013; Terrado et al., 2014).

The InVEST model operates on the principle of the Budyko theory (Budyko, 1958, 1974). Based on works of Schreiber (1904) and Ol'Dekop (1911), Budyko proposed formulations explaining the relationship between precipitation and potential evapotranspiration (PET) in order to couple water and energy balances, defined as the Budyko hypothesis. Several attempts have later been made to obtain an analytical solution of the Budyko hypothesis (Schreiber, 1904; Ol'Dekop, 1911; Turc, 1954; Mezentsev, 1955; Pike, 1964; Fu, 1981; Choudhury, 1999; Zhang et al., 2001, 2004; Porporato et al., 2004; Yang et al., 2008; Donohue et al., 2012; Wang and Tang, 2014; G. Zhou et al., 2015; S. Zhou et al., 2015). Among these studies, solutions provided by Fu (1981) called Fu's equation, gained significant attention as the work represented the effect of catchment properties on water balance components by incorporating an addition parameter “w”. Fu's equation can provide a full picture of the evaporation mechanism at the annual timescale. Therefore, Fu's equation can be used through a top-down analysis for providing insight into the dynamic interactions among climate, soils, vegetation, and their controls on the annual water balance at the regional scale (Yang et al., 2007).

Considering the lack of studies on analysis and validation of ES on the Indian subcontinent, especially for Himalayan catchments, and to assess the applicability of various water-balance models to Himalayan catchments, the present work attempts to compute and analyze water yield in the upper Ganga Basin using a semi-distributed InVEST model and a Lumped Zhang model. The work primarily considers, in detail, the spatial variation of InVEST model parameters and uses different strategies to compute water yield. Accordingly, water yield is estimated for 4 years, i.e., 1980, 1990, 2001 and 2015 and the most appropriate strategy is identified. The parameters that are adopted as lumped at the basin scale in previous studies are estimated at the pixel scale in order to avoid the dependence of the model parameters on size of the catchment. In addition, pixel-level estimations of water yield are expected to be more accurate than output obtained using the conventional approach with basin-lumped output. The term “finer scale” in the paper represents the incorporation of spatial variations through the pixel-level estimation of parameters involved in InVEST model, which are otherwise taken as lumped. The work also compares the outcomes of spatially distributed water yield models and the conventionally used Lumped Zhang model.

2 Background theory

2.1 Water yield models

In this section, two water yield models, i.e., the InVEST water yield model, which is a distributed model, and the Lumped Zhang model, are described.

2.1.1 InVEST model

The InVEST water yield model (Tallis et al., 2010) is designed to provide information regarding the changes in the ecosystem that are likely to alter the flow. It is based on the Budyko theory, which is an empirical function that yields the ratio of actual to potential evapotranspiration (PET) (Budyko, 1979). To describe the degree to which long-term catchment water balance deviates from the theoretical limits, a number of scholars have proposed one-parameter functions that can replicate the Budyko curve (Fu, 1981; Choudhury, 1999; Zhang et al., 2004; Wang and Tang, 2014). To observe and represent pixel-level changes to the landscape, InVEST model incorporates, explicitly, the spatial variability in precipitation, PET, soil depth and vegetation. The model operates at the grid scale and acquires the inputs in the raster format into a GIS environment such as ArcGIS.

The InVEST water yield model is based on an empirical function known as the Budyko curve (Budyko, 1974). Annual water yield, Y(x), is determined at each pixel of a landscape as follows;

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }Y\left(x\right)=\left(\mathrm{1}-{\displaystyle \frac{\mathrm{AET}\left(x\right)}{P\left(x\right)}}\right)\times P\left(x\right),\end{array}$$

where AET(x) is the actual annual evapotranspiration per pixel x and P(x) is the annual precipitation per pixel x. Actual evapotranspiration (AET) is essentially determined by climatic factors (precipitation, temperature, etc.) and is mediated by catchment characteristics (vegetation cover, soil characteristics, topography, etc.). On the other hand, potential evapotranspiration (PET) represents the evaporating potential of the climate system at a specific location and time of year without the consideration of catchment characteristics and soil properties (Allen et al., 1998). Several attempts have been made in the past to establish a relationship between AET and PET, among which the solution provided by Fu (1981) has been adopted worldwide. Fu (1981) provided an analytical solution to the Budyko hypothesis and related AET with PET by incorporating a dimensionless parameter “w”, which denotes the effect of catchment characteristics.

The InVEST model uses the expression of the Budyko curve proposed by Fu (1981) and Zhang et al. (2004). The ratio of mean annual PET to annual precipitation, known as index of dryness, is expressed as

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\mathrm{AET}\left(x\right)}{P\left(x\right)}}=\mathrm{1}+{\displaystyle \frac{\mathrm{PET}\left(x\right)}{P\left(x\right)}}-{\left[\mathrm{1}+{\displaystyle \frac{\mathrm{PET}\left(x\right)}{P\left(x\right)}}\right]}^{\left(\frac{\mathrm{1}}{w}\right)},\end{array}$$

where PET(x) is the annual potential evapotranspiration per pixel x (mm), and w(x) is a non-physical parameter that influences the natural soil properties. The PET(x) is calculated using the following expression;

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\mathrm{PET}\left(x\right)=K_{\mathrm{c}}\left(x\right)\times {\mathrm{ET}}_{\mathrm{0}}\left(x\right),\end{array}$$

where ET₀(x) is the annual reference evapotranspiration per pixel x, which is computed based on evapotranspiration from alfalfa grass grown at that location using Eq. (6). K_c(x) is the vegetation evapotranspiration coefficient that is influenced by the change in characteristics of land use and land cover at every pixel (Allen et al., 1998). The values of ET₀(x) are adjusted by K_c(x) for each pixel over the map of land use and land cover. w(x) is an empirical parameter, and the expression given by Donohue et al. (2012) for the InVEST model has been applied to define w(x), which is expressed as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }w\left(x\right)=z\times {\displaystyle \frac{\mathrm{AWC}\left(x\right)}{P\left(x\right)}}+\mathrm{1.25}.\end{array}$$

Thus, the minimum value of the parameter w(x) is 1.25, corresponding to bare soil where the root depth is zero (Donohue et al., 2012). The Donohue model was originally developed for Australia, however, the online documentation on InVEST model states its application globally. The parameter z is known as the seasonality factor whose value varies from 1 to 30. It represents the nature of local precipitation and other hydrogeological parameters. The parameter AWC(x) depicts volumetric plant available water content expressed in depth (mm), which can be expressed by following formula for each pixel x:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathrm{AWC}\left(x\right)=\mathrm{Min}.(\mathrm{Restricting}\mathrm{layer}\mathrm{depth},\mathrm{root}\mathrm{depth})\times \mathrm{PAWC}.\end{array}$$

The root-restricting layer depth is defined as the depth of the soil up to which the soil can allow the penetration of the roots, and root depth is defined as the depth where 95 % of the root biomass occurs. Plant available water content (PAWC) is generally taken as the difference between the field capacity and the wilting point. It depends upon the soil properties and can be computed by the Soil-Plant-Air-Water (SPAW) software. In the study, PAWC is calculated using the method described by McKenzie et al. (2003). The modified Hargreaves method and Hargreaves method were employed for computing reference evapotranspiration for the study area at pixel scale.

The modified Hargreaves method is expressed as

$$\begin{array}{ll}{\displaystyle }{\mathrm{ET}}_{\mathrm{0}}& {\displaystyle }=\mathrm{0.0013}\times \mathrm{0.408}\times \mathrm{RA}\times \left(T_{\mathrm{avg}}+\mathrm{17.0}\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }\times (\mathrm{TD}-\mathrm{0.0123}\times P{)}^{\mathrm{0.76}},\end{array}$$

where ET₀ is reference evapotranspiration, T_avg is the average daily temperature (^∘C) defined as the average of mean daily maximum and mean daily minimum temperature, TD (^∘C) is the temperature range computed as the difference between mean daily maximum and mean daily minimum temperature, and RA is extraterrestrial radiation (MJ m⁻² day⁻¹).

According to the Hargreaves method,

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{\mathrm{ET}}_{\mathrm{0}}=\mathrm{0.0023}\times \mathrm{0.408}\times \mathrm{RA}\times \left(T_{\mathrm{avg}}+\mathrm{17.8}\right)\times {\mathrm{TD}}^{\mathrm{0.5}},\end{array}$$

where terms involved in the equation means same as those in the modified Hargreaves method.

For computing the extraterrestrial radiation (RA), the following equation is used;

$$\begin{array}{ll}{\displaystyle }\mathrm{RA}& {\displaystyle }={\displaystyle \frac{\mathrm{24}\left(\mathrm{60}\right)}{\mathit{\pi }}}\times G_{\mathrm{sc}}\times d_{\mathrm{r}}\times \left[w_{\mathrm{s}}\sin \left(\mathit{\phi }\right)\sin \left(\mathit{\delta }\right)\right.\\ \text{(8)}& {\displaystyle }& {\displaystyle }\left.+\cos \left(\mathit{\phi }\right)\cos \left(\mathit{\delta }\right)\sin \left(w_{\mathrm{s}}\right)\right],\end{array}$$

where RA is extraterrestrial radiation (MJ m⁻² day⁻¹), d_r is the inverse Earth–Sun relative distance, G_sc is the solar constant equal to 0.0820 MJ m⁻² min⁻¹, w_s is sunset hour angle (rad), δ is the solar declination (rad) and φ is latitude (rad).

Determination of the parameter “w”

The dimensionless parameter w depends upon the local climatic variables such as the hydrological characteristics of the area, its rainfall intensity and topography. In the InVEST water yield model (Tallis et al., 2010), parameter w can be computed in three different ways. The first method is suggested by Donohue et al. (2012), in which parameter w is computed using Eq. (4) and where sensitivity parameter z is adopted as one fifth of the number of rain events per year. The second method is suggested by Xu et al. (2013), which compares w with latitude, the NDVI (normalized difference vegetation index), area, etc. The third method experiments with various selections of w (one value of w for the entire study region) until there is a good match between observed and computed water yield. Unfortunately, this method is not suited for a pixel-based analysis, as the number of pixels will be extremely large, making the method computationally intensive.

2.1.2 Lumped Zhang model

In this model, the mean value of different parameters is used as an input to compute the average value of the water yield for the whole watershed. The average actual evapotranspiration, potential evapotranspiration, w, precipitation, etc., are described by Zhang et al. (2004).

3 Study area and data

3.1 Study area

The Ganga river in India is ranked amongst the world's top 20 rivers in regards to the water discharge. The Ganga river is segregated into three zones, viz. the upper Ganga Basin, middle Ganga Basin and lower Ganga Basin. The area chosen for the present study, i.e., the upper Ganga Basin, is situated in the northern part of India within the geographical coordinates 29^∘48^′–31^∘24^′ N and 77^∘49^′–80^∘22^′ E, covering an area of 22 292.1 km² and reaching up to Haridwar. The altitude of the study area varies from 275 m in the plains to 7512 m in the Himalayan terrains. A region of approximately 433 km² of the basin is located under glacier landscape, and 288 km² of the region is located under a fluvial landscape. About 60 % of the basin is utilized for agricultural practices, and 20 % of the basin is in the forest area, especially in the upper mountainous region. Nearly 2 % of the basin is permanently covered with snow in the mountain peaks. The most predominant soil groups found in the region are sand, clay, loam and their compositions. In the upper Ganga Basin, the average annual rainfall varies from 550 to 2500 mm (Bharati et al., 2011), where a major fraction of total annual rainfall is received during monsoon months (June–September). The geographical location and other information of the upper Ganga Basin are represented in Fig. 1.

Figure 1 Graphical representation of the study area, the upper Ganga Basin.

3.2 Data

3.2.1 Precipitation and temperature

The daily time series of precipitation and temperature for the study area are acquired from India Meteorological Department (IMD) at a grid size of 0.25^∘ and 1^∘, respectively. The upper Ganga Basin comes within the dataset latitude, which ranges from 29.5^∘ N to 31.5^∘ N, and its longitude, ranging from 77.75 to 80.25^∘ E . The daily time series of precipitation was aggregated to obtain the annual time series at each grid point. Various analyses in the study are carried out for 4 years, i.e., 1980, 1990, 2001 and 2015.

3.2.2 Soil map

Spatial maps of soil were collected from the National Bureau of Soil Survey and Land Use Planning (NBSSLUP) at 1:250 000. Digital maps of soil available at a resolution of 1200 m × 1200 m were resampled to the resolution of land use data, i.e., 30 m × 30 m, using “resample” tool in ArcGIS in order to maintain the scale homogeneity. The attribute table of the raster layer contains fields like soil depth, soil texture, carbon content percentage, drainage, slope, erosion, soil temperature and mineralogy. The relevant features, i.e., soil depth and soil texture are converted into the raster image for the upper Ganga Basin.

Table 1 Value of K_c corresponding to the classes of land use and land cover.

3.2.3 Map of land use and land cover

Satellite images were acquired from different sensors of Landsat, viz. Landsat 3/4 Multispectral Scanner and Thematic Mapper (MSS/TM), Landsat 4 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper (ETM) and Landsat 8 Operational Land Imager (OLI) sensors for the years 1980, 1990, 2001 and 2015, respectively. The images are available at different resolutions and in several wavelength bands, from which green (G), red (R) and near-infrared (NIR) band images are combined to create a false color composite (FCC) for the study area in ERDAS Imagine. FCCs are then classified using supervised classification in ERDAS in six different classes, i.e., forest, water, agricultural, wasteland, snow and glacier, and built-up land. The classification of the area is based on their similar response under different bands. Each class is then recognized with the help of ground-truth and high-resolution satellite images.

4 Methodology

In the present work, five different strategies are employed to compute water yield. For the ease of presentation, these strategies are referred to as A–E. In strategy A, an average value of precipitation, temperature, extraterrestrial radiation and parameter w is used for the entire basin. This strategy is essentially based on Lumped Zhang model. Strategies B–E are designated, corresponding to a particular variation of the InVEST model where water yield is computed using different approach for estimating parameter w. For computing parameter w, relationships for large basins and for the global model from Xu et al. (2013) are given by Eqs. (9) and (10), respectively.

For large basins,

$$\begin{array}{ll}{\displaystyle }w& {\displaystyle }=\mathrm{0.69387}-\mathrm{0.01042}\times \mathrm{lat}+\mathrm{2.81063}\times \mathrm{NDVI}\\ \text{(9)}& {\displaystyle }& {\displaystyle }+\mathrm{0.146186}\times \mathrm{CTI}.\end{array}$$

For the global model,

$$\begin{array}{ll}{\displaystyle }w& {\displaystyle }=\mathrm{3.50412}-\mathrm{0.09311}\times \mathrm{slp}-\mathrm{0.03288}\times \mathrm{lat}+\mathrm{1.12312}\\ \text{(10)}& {\displaystyle }& {\displaystyle }\times \mathrm{NDVI}-\mathrm{0.00205}\times \mathrm{long}-\mathrm{0.00026}\times \mathrm{elev},\end{array}$$

where, “slp” is the slope gradient, “lat” is the absolute latitude of basin center, “CTI” is the compound topographic index, “NDVI” is the normalized difference vegetation index, “lat” is latitude, “long” is longitude and “elev” is elevation.

In strategy B, the entire basin is considered for computing the parameter w for large basins, using Eq. (9), which is given by Xu et al. (2013). In strategy C, the parameter w is computed for entire basin using Eq. (10), which is given by Xu et al. (2013). In strategy D, parameter w is computed at each pixel in order to incorporate the spatial distribution of the hydrologic variables involved in the computations. In Strategy E, parameter z is computed according to the number of rain events in a year; subsequently, Eq. (4) is used to compute the parameter w.

For all the strategies, the extraterrestrial radiation (RA) parameter is computed for each month using Eq. (8), and a raster layer is generated. Precipitation data are obtained from Indian Meteorological Department (IMD) at a grid size of 0.25^∘ for the study area. It has been interpreted and converted to the raster format by using the inverse distance weighted (IDW) interpolation technique in the ArcGIS environment for obtaining the values for all pixels at a resolution equal to the resolution of the Landsat satellite images. The temperature dataset is obtained from the IMD at a grid size of 1^∘ × 1^∘ for the study area and has also been converted to a raster format by using the IDW interpolation technique for obtaining the values for all pixels. Subsequently, the mean monthly value of average temperature (T_avg) and the difference between the mean daily maximum and mean daily minimum (TD) are obtained. The climate datasets used in the present study are of the finest resolution available so far for the study region. Gridded datasets of temperature and precipitation used in the present study have been developed using quality-controlled stations and well-proven interpolation techniques. Further details about the datasets of precipitation and temperature are given in Srivastava et al. (2009) and Pai et al. (2014), respectively.

The modified Hargreaves method is applied for obtaining the value of reference evapotranspiration at each pixel for each month (Droogers et al., 2002). To compute potential evapotranspiration, the yearly values obtained for the reference evapotranspiration are multiplied by the vegetation evapotranspiration coefficient (K_c), which depends on the LULC characteristics, as expressed in Eq. (3). The value of K_c is taken from Allen et al. (1998), as shown in Table 1. In this study, K_c is taken in the same was for all 4 years, as shown in Table 1, and is used to obtain potential evapotranspiration, which is subsequently used to obtain annual water yield at each pixel of the study area.

Figure 2 Reference evapotranspiration (mm) of the upper Ganga Basin for the years 1980, 1990, 2001 and 2015.

5 Results

5.1 Reference evapotranspiration, ET₀(x)

Reference Evapotranspiration (ET₀) is computed for the upper Ganga Basin using a high-resolution monthly climate dataset. The modified Hargreaves method is applied for obtaining the values of reference evapotranspiration at each pixel for each month (Droogers and Allen, 2002). ET₀ is a function of RA, precipitation, T_avg and TD, which are computed pixel-wise for each month of the years 1980, 1990, 2001 and 2015. Some of the months, i.e., July, July and August 1990; June, July and August 2001; and June, July and August 2015, showed negative values of reference evapotranspiration from applying the modified Hargreaves method. For these months, Hargreaves method is applied for obtaining the positive values. Subsequently, all mean monthly values are added up to get the mean annual values of evapotranspiration for the years 1980, 1990, 2001 and 2015, as represented in Fig. 2.

5.2 Potential evapotranspiration, PET(x)

The annual values obtained for the ET₀ are multiplied by the vegetation evapotranspiration coefficient (K_c), which varies with the characteristics of land use and land cover, as expressed in Eq. (3). The value of the K_c is taken from Allen et al. (1998). The values of the vegetation evapotranspiration coefficient are taken from Table 1. Thus, the potential evapotranspiration is computed for upper Ganga Basin for the years 1980, 1990, 2001 and 2015, as represented in Fig. 3.

Figure 3 Potential evapotranspiration (mm) of the upper Ganga Basin for the years 1980, 1990, 2001 and 2015.

5.3 Water yield, Y(x)

As described in the methodology, five different strategies, viz. A–E, are used to estimate water yield for the upper Ganga Basin.

Strategy A: water yield computed using the Lumped Zhang model

Here, the basin average values of all the input parameters are considered, and water yield is computed for the upper Ganga Basin for the years 1980, 1990, 2001 and 2015, which are obtained as 658.52, 925.68, 603.71 and 1194.25 mm, respectively.

Strategy B: water yield obtained by taking the single weighted mean value of parameter “w” from Xu et al. (2013) for large basins

In this strategy, water yield is computed by considering a single value of the parameter w for the whole basin using Eq. (9). The weighted mean value for parameter w for the years 1980, 1990, 2001 and 2015 are obtained as 1.507, 1.541, 1.403 and 1.507, respectively. The spatial distribution of the water yield for the upper Ganga Basin computed using strategy B is represented in Fig. 4. The mean values of water yield as obtained using this method for the years 1980, 1990, 2001 and 2015 are 755.65, 959.48, 742.39 and 1131.42 mm, respectively.

Table 2 Comparison of model-estimated PET and AET with a global dataset from different sources.

Figure 4 Water yield obtained by taking the single weighted mean value of parameter w from Xu et al. (2013) for large basins.

Strategy C: water yield obtained by taking a single weighted mean value of parameter “w” from Xu et al. (2013) for the global model

In this strategy, water yield is computed by considering a single value of parameter w for the entire upper Ganga Basin using Eq. (10). The weighted mean value of parameter w for the years 1980, 1990, 2001 and 2015 are obtained as −0.967, −0.955, −1.010 and −0.968, respectively. The spatial distribution of the water yield for the upper Ganga Basin as computed using strategy C is shown in Fig. 5. The mean values of water yield for the years 1980, 1990, 2001 and 2015 are 1239.92, 1549.46, 1149.93 and 1754.59 mm, respectively.

Strategy D: water yield obtained using the pixel-level estimation of parameter “w” from Xu et al. (2013)

In this strategy, the values of parameter w are estimated at the pixel level. The water yield computed for the years 1980, 1990, 2001 and 2015 for upper Ganga Basin is shown in Fig. 6. The mean values of water yield as computed using strategy D for the years 1980, 1990, 2001 and 2015 are 1240.02, 1549.44, 1149.89 and 1754.62 mm, respectively.

Strategy E: water yield obtained using the pixel-level estimation of parameter “w” from Donohue et al. (2012)

Equation (4) represents the parameter w as a function of parameter “z”, AWC and precipitation. The parameter w in the equation used in strategy E has been proposed by Donohue et al. (2012), which is also cited in online documentation of InVEST model; however, the final equation used for estimating water yield is obtained from the InVEST model. Considering this fact, Donohue et al. (2012) has been cited in strategy E. The water yield as computed using strategy E for the upper Ganga Basin for different years is shown in Fig. 7. The mean values of water yield for the years 1980, 1990, 2001 and 2015 are 1241.09, 1552.38, 1153.95 and 1753.53 mm, respectively.

Table 3 Observed vs. computed water yield for various proposed strategies for Rishikesh sub-basin.

Figure 5 Water yield obtained by taking the single weighted mean value of parameter “w” from Xu et al. (2013) for the global model.

5.4 Validation of ET and water yield estimates

For validation of model outputs, the basin's average annual values of PET and AET estimated using various strategies are compared with the corresponding basin average values obtained from available global datasets (Table 2). Model-simulated AET values are obtained from the Global Land Data Assimilation System (GLDAS) ET dataset from Noah model outputs. Basin average values of PET are obtained from the Climate Research Unit's (CRU's) PET datasets (CRU TS v. 4.01) available at resolution of 0.5^∘. From the comparison, both AET (GLDAS) and PET (CRU TS) values are found to be in fair agreement with the globally estimated values (Table 2). Spatial maps of global datasets of AET and PET are shown in Figs. 8 and 9, respectively.

The validation of water yield obtained from various strategies is performed at the Rishikesh gauging site of the upper Ganga Basin (Fig. 10). The discharge data of the basin are obtained from Irrigation Department of the state of Uttarakhand. The discharge observed in the basin is generated from precipitation as well as snowfall in the region, where 32 % of the discharge has been removed, because it is contributed to by glacier ice melt, as explained by Maurya et al. (2011) for our study area. The aforementioned fraction of discharge had been quantified using an isotope study that separates the contribution of glacier melt in quantifying discharge (Maurya et al., 2011). A comparison of the water yield computed and observed for the study region for different years by various proposed strategies is shown in Table 3.

Figure 6 Water yield obtained by computing pixel-wise value of parameter w from Xu et al. (2013).

As can be seen in Table 3, values of water yield estimated using strategies A to E are systematically increasing but are not steady in nature, as water yield estimated using strategy A and B lies in the range 650–750 mm, whereas water yield from strategies C–E lie in range of 1229–1231 mm for the years 1980 (see Table 3). Similar results are also evident for other years, too. Also, water yield estimated using strategies C–E are more or less the same for a given year, because these strategies involve pixel-based estimations of water yield considering spatial variation in the Budyko parameters. The parameters involved in the Budyko model, such as w, are dependent on various factors, such as catchment characteristics, vegetation cover, etc., as well as climate seasonality (Li et al., 2013). Ahn and Merwade (2017) have analyzed the relationship between basin characteristics and parameter w for 175 stations spread across the USA. Considering their study, no precise conclusion can be drawn regarding relationship between basin characteristics and the value of parameter w, especially in the case of basin-area characteristics. Moreover, no definite relationship has been yet identified between basin characteristics and model parameters, and this is a subject matter for further study.

6 Discussion

The study aimed to apply the InVEST water yield model to compute the water yield for upper Ganga Basin having highly variable topography consisting of hilly, plain and snow-covered areas. The InVEST model is based on the Budyko theory, which requires low amounts of data and low levels of expertise, thus making it acceptable worldwide. The mean monthly precipitation, temperature, monthly value of difference of the mean daily maximum and mean daily minimum, and extraterrestrial radiation parameters for the upper Ganga Basin of all 4 years, i.e., 1980, 1990, 2001 and 2015, are converted into the raster format for various analyses. The monthly reference evapotranspiration is thus computed using input parameters in GIS environment by applying the modified Hargreaves equation for all the months, except for a few months in which the modified Hargreaves equation gives negative results for the reference evapotranspiration. For those months, the Hargreaves method is applied to obtain the positive value of reference evapotranspiration, as also suggested by Goyal and Khan (2017). Reference evapotranspiration when multiplied with K_c gives the potential evapotranspiration. All monthly values are added up to obtain the annual value of reference evapotranspiration. K_c is a function of land use and land cover; thus, supervised classification is done to prepare the raster map of land use and land cover for the upper Ganga Basin. Subsequently, the annual value of potential evapotranspiration is obtained for the study area for the years 1980, 1990, 2001 and 2015.

Figure 7 Water yield obtained by computing pixel-wise value of parameter “w” from Donohue et al. (2012).

Figure 8 Spatial distribution of AET obtained from GLDAS Noah output datasets.

Figure 9 Spatial distribution of PET obtained from CRU datasets.

Figure 10 Graphical representation of sub-basin Rishikesh.

The paper employs various methodologies for water yield estimation, as discussed in the methodology section for the upper Ganga Basin. Thus, water yield is computed both from the InVEST model as well as the Lumped Zhang model. The value of the parameter w is computed using four different approaches, i.e., the mean single value obtained from Xu et al. (2013) for large basins, mean single value obtained from Xu et al. (2013) for the global model, pixel-level estimated value of parameter w from Xu et al. (2013) and pixel-wise value of parameter w from Donohue et al. (2012). Although the upper Ganga Basin lies in large basin category as per the definition from Xu et al. (2013), the yield computed using global model is in good agreement with the observed data for the region. In the study, the pixel-level estimation of parameter w is made in order to incorporate the spatial variability of the parameter involved in water yield estimation. Thus, two pixel-wise values of parameter w are computed for the upper Ganga Basin for years 1980, 1990, 2001 and 2015 by considering two approaches given by Xu et al. (2013) and the approach given by Donohue et al. (2012). Also, the basin-lumped water yield is computed using Lumped Zhang model, which considers the single mean values for entire basin of all the parameters involved in the computation of water yield. The water yield is computed in five different ways for the upper Ganga Basin for the years 1980, 1990, 2001 and 2015.

At the Rishikesh gauging site, surface runoff data are obtained by extracting the snowmelt from the discharge data, as the melting snow contributes about 32 % of total runoff in the Himalayan basins (Maurya et al., 2011). For validating the water yield obtained from different strategies, the observed yield is compared with the computed water yield based on different proposed strategies for the years 1980, 1990, 2001 and 2015, as represented in Table 3. The results obtained from Donohue et al. (2012) and Xu et al. (2013) are computed at pixel level (Strategy C–E); thus, they exhibit better performance than other approaches and are in good agreement with the observed data. These results exhibit the superiority of pixel-level computation to hydrological analyses for a watershed. The parameters involved in the Budyko model are dependent on various factors, such as basin characteristics (size, topography, stream length, slope, etc.), climate seasonality, etc. (Li et al., 2013). Again, the factors affecting model parameters vary both spatially and temporally. Moreover, the relationship between these factors and model parameters are not yet well defined (Ahn and Merwade, 2017). In such scenarios, adopting a hypothesis by assuming either of these controlling factors (such as “w”) to be spatially or temporally constant is inappropriate. Considering these facts, the present study attempts to incorporate the spatial variability of model parameter for estimation of water yield at the pixel level. As the computations are made at pixel level (on a grid of size 30 m × 30 m), the assumption of dependence of model parameters on the size of the catchment may also be disregarded. The computations made in the present work are based on empirical equations; however, the application of these equations has been well documented worldwide for estimations of various water balance components at various basin scales (Zhang et al., 2008; Ma et al., 2008; Ning et al., 2017; Rouholahnejad Freund and Kirchner, 2017; Wang and Zhou, 2016). Hence, it is recommended that for such a large basin, it is required to compute all the parameters involved in the computations of water yield at the pixel scale rather than adopting mean values for entire watershed.

7 Summary and conclusions

The present study aimed to apply the InVEST annual water yield model, a tool that is gaining interest in the ecosystem services community, in the upper Ganga Basin. While such simple models have low requirements for data and level of expertise, practical applications of such a model with single representative values of the model parameter for the entire basin do not provide accurate estimates of water yield. Performing pixel scale computation of water yield in the study indicates a better performance, and the results obtained show better agreement with the observed water yield. As far as parameter w is concerned, the global model works better than other representations of parameter w available in literature, especially in the upper Ganga Basin. In the study, the water yield is computed using five different strategies, and results are validated with the observed data at the outlet of the upper Ganga Basin. The present study attempts to quantify annual water yield at the pixel level, making the computations independent of the size of catchment. Therefore, the proposed methodology is expected to perform well for a catchment of any given size. Changes in catchment water storage over time are required to be quantified in order to validate the applicability of Budyko's model to long-term data for the studied catchment. Earlier, some of the important parameters defining water yield used to be computed at a basin-level scale, which caused errors in the results.

The study attempts to incorporate the spatial variability of parameters involved in the model through the pixel-level estimation of parameters that are otherwise taken as lumped in the previous studies. Study results show that the estimated water yield, considering spatial variability in model parameters, is in better agreement with the observed water yield compared to the water yield estimated when considering the parameters to be lumped over the study region. Further, the computations of various parameters are made at the pixel level; therefore, the estimates of water balance components using this approach are expected to be independent of the assumption of dependence of parameters on catchment size. As the relationship between Budyko's model parameters and their controlling factors has not been well defined (Ahn and Merwade, 2017), the study emphasizes water yield estimation using pixel-based computations. The study outcomes can be summarized as follows: (i) between two approaches used in the study, i.e., considering the entire basin and pixel-level approach, the pixel-level approach is found to provide better results; and (ii) in pixel-based computations, results are further improved with the use of a parameter w based on a global model rather than regional models of parameter w, especially for large basins in the Himalayan region.