2.1 Study domain, river gauging stations and catchment boundaries

The study domain includes the river basins that span India, including the catchment areas that fall outside of the political boundaries of India (Fig. 1). The boundaries of the river basins used here are generally consistent with those used by India's Central Water Commission (CWC). Consistent with the CWC, adjacent watersheds in some regions were pooled to create composite river basins, such as west-flowing rivers (WFR) north and south, east-flowing rivers (EFR) north and south, and west-flowing rivers of Kutch (WFR Kutch). The catchment boundaries used here are from the GHI dataset (Goteti, 2023), a quality-controlled dataset on India's river gauging stations, catchment boundaries and hydrometeorological time series. However, the GHI dataset is limited to Peninsular India. The catchment boundaries for the Northern Indian watersheds were derived using the HydroSHEDS suite of products, using the same procedures as the GHI dataset. Station descriptions available from the CWC were validated using online maps (e.g., Google Maps). Stations were then relocated to the closest point on the river network. The watershed draining into this relocated point and all of the upstream watersheds were recursively identified using geographic information system (GIS) software. Catchment areas for the delineated watersheds were validated against those reported by the CWC.

Figure 1(a) Major river basins spanning India. The shaded basins are transboundary basins. Stations with daily streamflow data (grey dots) are from the GHI dataset (Goteti, 2023), and stations with only annual streamflow data (black dots) are from CWC-19 (2019). (b) River basins of Peninsular India. Some basin names have been shortened for ease of display within the map, with the complete names shown next to the map.

The river basins of Peninsular India, the non-shaded region in Fig. 1a, have daily streamflow data that are available through India's Central Water Commission (CWC). Limited streamflow data are available for the river basins of Northern India (shaded regions in Fig. 1a). The stations used here were chosen such that the catchment area discrepancy between the GHI and CWC datasets is less than 5 %, and there was at least 5 years of observed streamflow data with minimal missing records. A total of 242 stations are used in this analysis, with 213 of these stations located in Peninsular India and 29 in Northern India (dots in Fig. 1a). The number of stations within each basin and other pertinent information are summarized in Tables S2.1 and S2.2.

2.2 Precipitation

The selected P datasets used here are outlined in Table 1 and are briefly described here. In addition to these datasets, the PBCOR dataset is used as a reference climatology in certain parts of this analysis. Additional information on the PBCOR dataset is provided in Sect. S1, while additional information on P datasets is presented in Sect. S3.

Table 1The P datasets analyzed here and relevant information. The input used in the creation of each dataset is indicated by one or more of “G” (gauge), “O” (observation-based data product), “R” (reanalysis) and “S” (satellite). In the time span, the end year is left blank when a dataset extends to the present.

^a Asian Precipitation – Highly-Resolved Observational Data Integration Towards Evaluation of Water Resources. ^b Climate Hazards Group InfraRed Precipitation with Station data. ^c Climate Prediction Center Morphing Technique. ^d European Centre for Medium-Range Weather Forecasts (ECMWF) land component of the fifth generation of European ReAnalysis (ERA5). ^e Global Satellite Mapping of Precipitation. ^f Indian Monsoon Data Assimilation and Analysis reanalysis. ^g Integrated Multi-satellitE Retrievals for Global precipitation measurement. ^h Multi-Source Weighted-Ensemble Precipitation. ⁱ Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks (Cloud Classification System–Climate Data Record). ^j Soil Moisture to Rain (Advanced Scatterometer, v1.5).

Download Print Version | Download XLSX

The P datasets used here were often identified in the recent literature as being reasonable representations of the observed P, and they range in spatial resolution from about 4 to 25 km and in temporal frequency from half an hour to a month. Datasets included here are based on rain gauges (e.g., the IMD dataset), reanalysis (e.g., ERA5-Land), satellites, or a combination of sources (e.g., CHIRPS). The IMD gauge-based dataset is of primary interest here, since it is an often-used benchmark that is employed in a number of studies.

The reader should note that while the IMD dataset is limited to India's political boundaries, the rest of the P datasets are not. However, certain river basins of India extend beyond India's boundaries and are part of this analysis. To enable an appropriate comparison between datasets, the IMD dataset is complemented, where needed, with the APHRODITE dataset (Yatagai et al., 2012). The APHRODITE dataset was chosen for several reasons: it is also based on rain gauge data, similar to the IMD dataset; its spatial and temporal resolution are the same as the IMD dataset's resolution (0.25° or ∼25 km and daily); and studies in the literature have found that APHRODITE compares reasonably well with the IMD dataset across many parts of India (e.g., Prakash et al., 2015b). While limitations with APHRODITE are discussed by such studies, it is assumed to be the best gauge-based alternative to the IMD dataset.

For those regions where data from the IMD are unavailable, grids from APHRODITE were identified, and then the data from those grids was interpolated to align with the IMD dataset grid. Finally, a blended product called IMD-APHRO which spanned the entire study domain was created. In the remainder of this paper, unless otherwise stated, IMD-APHRO refers to the blended product created here, and IMD refers to the product confined to India's political boundaries. Also, in the remainder of this paper, each P product is referred to by its alias (see Table 1).

2.3 Evapotranspiration

Two datasets were considered for this analysis, based on their usage in studies across India (Table 2) – the Numerical Terradynamic Simulation Group (NTSG) at the University of Montana (Zhang et al., 2010) and the Global Land Evaporation Amsterdam Model (GLEAM) (Martens et al., 2017; Miralles et al., 2011) datasets. GLEAM provides estimates of the different components of ET, including transpiration, bare-soil evaporation, interception loss, open-water evaporation and sublimation. A comparison of NTSG and GLEAM indicates that they are generally consistent with each other across several basins (Sect. S4). However, estimates from GLEAM tend to be lower than those from NTSG. GLEAM was the primary dataset used here because of its longer time span and its availability up to the present time.

Table 2The evapotranspiration datasets analyzed here and relevant information. In the time span, the end year is left blank when a dataset extends to the present.

^a Numerical Terradynamic Simulation Group (NTSG) Process-based Land Surface evapotranspiration/Heat (PLSH) fluxes algorithm. ^b Global Land Evaporation Amsterdam Model.

Download Print Version | Download XLSX

2.4 Other data

2.4.1 Elevation, land cover and land use

Figure 2a shows the variability in elevation across the study domain. The dominant features include the Himalayas in the northern and northeastern regions, the mountains (or ghats) along the western and eastern coasts of India, the plains of the Ganga and Brahmaputra basins, and the Deccan Plateau in Peninsular India.

Figure 2(a) Elevation based on HydroSHEDS 500 m topographic data. For ease of display, 500 m data are aggregated to 10 km resolution and the maximum value within each 10 km grid is displayed. (b) Major land-cover and land-use types based on PROBA-V 100 m data. For ease of display, 100 m data are aggregated to 10 km resolution and the dominant value within each 10 km grid is displayed.

The high-resolution (100 m) global dataset based on the PROBA-V satellite (Buchhorn et al., 2020) was used to identify the dominant land-cover and land-use types. Figure 2b shows the dominant land-cover and land-use types. For the purpose of this analysis, the land-cover types of grass, shrubs, and trees/forest were pooled into one category.

2.4.2 Water management

The water management considered here includes groundwater extraction, diversions (imports and exports), and reservoir storage, which are summarized in Fig. 3. A detailed description of these data is provided in Sect. S5. Groundwater extraction and recharge estimates are available from India's Central Ground Water Board (CGWB) for select years. The extent of annual groundwater extraction is quantified as a fraction of the annual P. Similarly, basin-scale imports and exports from CWC-19 (2019) are expressed as a fraction of the annual P. Information on large dams and reservoirs in India was obtained from the National Registry of Large Dams (NRLD, 2019). For each of the 242 watersheds used here, the cumulative live storage capacity from all dams present within the watershed is expressed as a fraction of the annual P.

Figure 3(a) Districtwise annual groundwater extraction as a percent of the annual P; (b) basin-scale imports and exports as a fraction of the annual P; (c) density of dams and reservoirs, represented as the number of dams per 0.25° (about 25 km) grid; and (d) cumulative maximum live storage capacity of each watershed expressed as a fraction of the annual P for WY 2019.

Figure 3a shows that groundwater extraction can be a substantial fraction of the annual P in certain parts of India, but it is minimal in the mountainous and wet regions of India – the western coast of India, northernmost India and Northeastern India. Similarly, water diversions are highest in the agricultural regions of the Ganga Basin and the interior parts of Peninsular India. The highest density of dams is in arid Western India, while the lowest density occurs in the plains and mountains of Northern India. There are some watersheds in coastal Peninsular India, where reservoir storage is a significant portion of the annual P, but most of the other watersheds are minimally impacted by such storage.

2.4.3 Changes in terrestrial water storage (TWS)

Changes in terrestrial water storage (TWS), inferred from the Gravity Recovery and Climate Experiment (GRACE) satellite mission (Tapley et al., 2004), are useful for identifying regions where large-scale water management is causing substantial changes to the natural hydrologic cycle (e.g., Famiglietti, 2014; Rodell et al., 2009). TWS includes water stored below the ground, on the ground and above the ground. GRACE-based TWS anomalies from the Center for Space Research (CSR) (Save et al., 2016; Save, 2020) were used to estimate the change in annual TWS (or ΔTWS) as a fraction of the annual P (see Sect. S6). Figure 4 shows the maximum and minimum ΔTWS over the period WY 2002–2014. The magnitudes of such changes for most of the study domain are within +20 % or −20 % of the annual P. However, there are regions, such as Northwestern, Northern and Eastern India, where the magnitudes of such changes are larger than 20 % of the annual P.

Figure 4(a) Gridwise maximum values of annual ΔTWS for WY 2002–2014. (b) Same as (a) but for the minimum values.

3.1 Water imbalance scenarios

In order to take advantage of the datasets on TWS anomalies and water management discussed in Sect. 2.4, the traditional annual water balance equation is formulated in two different ways in the following discussion.

Under natural circumstances, one could express the annual water balance of a watershed by assuming that the net change in terrestrial water storage (Δ_TWS) is the imbalance between the total actual P (P_act), the output fluxes of R and ET, and inter-watershed groundwater flow (IGF):

TWS is the sum of all the potential water reservoirs – groundwater, soil moisture, snow water equivalent, surface water, land ice and water in the biomass (Humphrey et al., 2023). Watershed boundaries do not always coincide with underlying aquifer boundaries, and IGF could play an important role in the watershed's water balance (e.g., Fan, 2019; Liu et al., 2020). However, in the absence of field data on groundwater flow pathways, it is not possible to quantify the effect of IGF on the water balance. IGF is assumed to be negligible. The implications of this assumption are discussed in Sect. 5.1.2.

Figure 6Schematic illustrating potential spurious scenarios derived from Eqs. (3) and (4) when P_act is underestimated. The terms in grey text within the annual water balance equations are relatively small compared to P_act.

Download

If one were to account for the effects of management, Δ_TWS would represent changes due to both natural and human-related causes such as groundwater extraction, reservoir storage and diversions. Under such circumstances, after ignoring IGF, one could reformulate Eq. (1) as Eq. (2). Net surface water diversions are represented by two terms: Exports (net loss of water) and Imports (net gain of water). The terms P_act, R, ET, Exports and Imports are non-negative. Δ_TWS is positive if there is a net increase in TWS and negative if there is a net decrease in TWS.

Rearranging Eq. (2) results in Eq. (3). The equality in Eq. (2) has been replaced with an approximation in Eq. (3) because the data needed, if available, are often not at the spatial or temporal resolution required to accurately balance the water budget.

There is another way, although more approximate than Eq. (3), of formulating the annual water balance. Management of water is present in many parts of India and includes groundwater extraction, reservoir storage and diversions (CWC-19, 2019). To take advantage of these data on management, the annual water balance is approximated as Eq. (4). Groundwater storage changes – both natural (Δ_GW natural) and human-caused changes (Δ_GW human) – are included. Changes to reservoir storage (Δ_Reservoir) and diversions (Exports and Imports) are also explicitly included. In Eq. (4), both Δ_GW terms are positive if there is a net aquifer recharge and negative if there is a net aquifer depletion. Thus, the groundwater extraction presented in Fig. 3 would be a negative quantity. Δ_Reservoir is positive if there is a net increase in reservoir storage and negative if there is a net decrease in storage.

The reader should note that Eqs. (3) and (4) are separate, but useful, ways of analyzing the water budget. While Δ_TWS in Eq. (3) includes changes in all potential water reservoirs, Eq. (4) is an approximation and does not adequately capture the effect of snow processes, does not include water stored as soil moisture, and does not capture all the effects of management. The reader should also note that hydrologic analyses often make the a priori assumption that the net annual change in storage (Δ_TWS or Δ_GW) is negligible. This study does not make such an assumption within Eqs. (3) and (4).

If UoP is absent (i.e., P_obs≈P_act), then, based on Eqs. (3) and (4), it is reasonable to expect R to be only a portion of P_act, regardless of the extent of management. If the effects of management – the two rightmost terms in Eq. (3) and the four rightmost terms in Eq. (4) – are relatively small compared to P_act, then it is also reasonable to expect P_act to approximately equal $R+\mathrm{ET}+\mathrm{\Delta }$, where Δ is either Δ_TWS oder Δ_GW natural. As discussed later in this section, for most watersheds in the study domain, a reasonable upper bound on the magnitude of Δ_TWS (and Δ_GW natural) is 20 % of P_obs. The above expectations are illustrated by the “likely scenarios” in Fig. 6.

If UoP is present (i.e., P_obs≪P_act), one could potentially realize the “spurious scenarios” of P_obs≤R and $P_{\mathrm{obs}}\ll R+\mathrm{ET}$ (see Fig. 6) when the extent of management is minimal. If, on the other hand, management is moderate to extensive, it is difficult to generalize the relationship between the relative magnitudes of P_act, R and ET since R and ET are no longer constrained by the natural water balance. The spurious scenarios in this situation only include the case where P_obs≤R. The word “minimal” is used in a relative sense when the overall effect of annual management at the watershed scale relative to the annual P_obs is minimal. It should not be interpreted as the effect of local management on specific storm events.

The two specific scenarios investigated here are based on the spurious scenarios in Fig. 6. The annual P_obs is less than or equal to R in Scenario I.

Thus, the annual runoff coefficient is at least 1 in Scenario I. This scenario could be realized regardless of the extent of management outlined in Fig. 6. Such a scenario was also used by other studies (e.g., Beck et al., 2020) to identify UoP. However, instances where the annual runoff coefficient is less than 1 but still spuriously high (e.g., 0.95) are excluded by Scenario I. Scenario II attempts to include such instances.

Moreover, Scenario II is also intended to capture instances where the sum of R and ET greatly exceeds P_obs. If UoP is present, relatively high values of R combined with reasonable estimates of ET result in the sum of R and ET greatly exceeding P_obs. The formulation of Scenario I exactly follows the first of the spurious scenarios in Fig. 6. The second of the spurious scenarios in Fig. 6, $P_{\mathrm{obs}}\ll R+\mathrm{ET}$, is not an exact mathematical relationship. It is made exact by the use of heuristics. The rationale behind such heuristics is presented in the following discussion.

The typical wet season runoff coefficient for the whole of India was estimated to be about 0.38 by Gupta et al. (2016). The basin-scale average annual runoff coefficient was estimated by Xiong et al. (2022) to range from 0.10 to 0.40 for several large river basins of India, with higher coefficients for the Indus and Brahmaputra basins. Considering the magnitudes of those estimated runoff coefficients, a coefficient of 0.70 was assumed to be a reasonable lower bound for identifying spuriously high annual runoff coefficients.

As shown in Fig. 4, for regions that have hilly terrain or are covered by forests, the magnitude of Δ_TWS is typically within 20 % of the annual P_obs. Watershed management is represented by the four rightmost terms in Eq. (4): Δ_{GW net recharge}, Δ_Reservoir, Exports and Imports. In the regions that have hilly terrain or are covered by forests (Fig. 2), where management can be assumed to be minimal, the magnitude of the individual effect of each type of management is typically less than 5 % of P_obs (Fig. 3). A reasonable upper bound on the cumulative effect of the four rightmost terms in Eq. (4) is also 20 % of the annual P_obs. Thus, when management can be considered minimal, it is reasonable to expect R+ET to have a maximum value of 1.20×P_obs. This is the justification for the heuristic of 1.20 in Scenario II. As mentioned earlier, “minimal” management is used in the context of the overall effect of annual management at a watershed scale relative to annual P_obs. For instance, a 20 % management effect of the annual P_obs in a watershed with 0.4 as the runoff coefficient translates to a 50 % ($\mathrm{20}\mathit{\%}/\mathrm{0.4}=\mathrm{50}\mathit{\%}$) effect on the annual R. Thus, minimal management could still have a substantial effect on the annual R.

This study identifies UoP by first identifying individual years within watersheds where Scenario I or II is realized. Then, it proceeds to investigate the extent of management and extent of Δ_TWS within those imbalanced watersheds. If management and Δ_TWS are deemed minimal relative to annual P, then it is concluded that the likely cause of the spurious imbalance is UoP.

3.2 Time series compilation

In order to investigate the abovementioned scenarios, annual time series of all the relevant terms need to be compiled. All of the variables needed are expressed in the same units of volume. Observed daily streamflow (R), available in units of m³ s⁻¹ was aggregated to cumulative monthly and annual volumes of million m³ s⁻¹ (MCM per month and MCM per year, respectively). Gridded data on P_obs and ET, available in units of depth per unit area per month (e.g., mm per month), were also aggregated to watershed-scale monthly and annual volumes. The process of aggregating grid-based products to a watershed involves identifying the spatial overlap between the grids and the watershed. Such relationships were identified using a GIS analysis. Grid-specific fractional areas were used in the process of aggregation. A schematic illustrating the process of aggregation is shown in Sect. S7. The time series needed were compiled for each of the 242 watersheds analyzed here. P datasets are often available up to the current year, but the latest year for which observed R data are available is WY 2017, and ET data since WY 1980 are available. The time span of the data compiled here is WY 1980 to 2017 (38 WYs), whenever data are available.

4.1 Imbalanced watersheds identified using IMD-APHRO

An example of Scenario I is shown in Fig. 7 for the Bantwal station on the Nethravathi River in the WFR south basin. The annual time series is shown in Fig. 7a, while the monthly time series for select years is shown in Fig. 7b. There are several WYs where the total annual volume of P is less than the total observed R, such as WYs 2011–2013 in the recent past. The total live capacity of all upstream reservoirs is 0 since there are no dams in this watershed. The monthly time series is also shown in for select years (WYs 2011–2015; Fig. 7b). The strong seasonal pattern imposed by the summer monsoon is evident, with the months of June to September having the highest values of P and R. There are several months within each year where the observed R is greater than P. It is useful to note that the above spurious relationship in which the annual R exceeds the annual P for the Bantwal watershed was also tabulated by CWC-19 (2019) (see their Appendix R, Table R-2), based on the same P and R data sources as those used here.

Figure 7(a) Example showing the imbalance scenario P≤R (Scenario I) for Bantwal station on the Nethravathi River in the WFR south basin. The annual R (grey bars) exceeds annual P (blue line) in certain years. ET (green bars) and cumulative reservoir storage capacity (red line) are also shown for reference. (b) Monthly volumes instead of annual volumes for select WYs. Months are indicated by initials and follow the June–May WY convention.

Download

Figure 8(a) Imbalanced watersheds from Scenario I (grey-shaded regions) and the gauging stations (blue dots) at the outlets of those watersheds. Watersheds shaded in pink have at least 1 % of their catchment area outside of India. (b) Same as (a) but for Scenario II.

Watersheds where either Scenario I or II was realized were identified by analyzing the annual P, R and ET for all 242 watersheds. Figure 8 shows the catchment areas corresponding to these imbalanced watersheds (grey areas) and the gauging stations at the outlets of these watersheds (blue dots). These watersheds are located along the western coast of India, in the forested and hilly regions of central India, and within the Himalayan mountains and their foothills. The locations of these imbalanced watersheds coincide with the regions receiving the highest annual P (see Fig. S1.2). Most major river basins have at least one imbalanced watershed. Some of these watersheds have catchment areas that are outside of India's political boundaries. Such watersheds with at least 1 % of the total catchment area outside of India are shown in pink in Fig. 8. Due to the limited availability of observed R data in Northern India, only a small number of imbalanced watersheds could be identified. In contrast, Peninsular India has many more imbalanced watersheds.

The watersheds identified above are based on a specific set of heuristics (0.70 and 1.20) within Scenario II. In order to understand the impacts of changing the heuristics, three other sets of heuristics were tried. Instead of 0.70, values of 0.60, 0.80 and 0.90 were used, and instead of 1.20, values of 1.10, 1.30 and 1.40 were used. Figure S9.1 shows the imbalanced watersheds resulting from the use of each set of heuristics. By lowering these heuristics, one would expect more watersheds to be categorized as imbalanced, while raising them would result in fewer watersheds. As expected, lower values of the heuristics (e.g., 0.60 and 1.10) result in a larger number of watersheds, and higher values (e.g., 0.90 and 1.40) result in a lower number of watersheds compared to the watersheds shown in Fig. 8. However, the general locations of these watersheds remain the same – the western coast of India, the forested and hilly regions of central India, and the Himalayan mountains and their foothills.

4.2 Characteristics of imbalanced watersheds

The dominant physical characteristics associated with these imbalanced watersheds are summarized in Fig. 9. The sizes of these watersheds can range from more than a 100 000 km² in the northern portion of the study domain to less than a 1000 km² in Peninsular India (Fig. 9a). The maximum elevation within such watersheds is about 2000 m (Fig. 9b), which is much higher than the average elevation of India – about 600 m (estimated in this study). The statistics on fractional land cover and land use indicate that most of these watersheds are predominantly covered by natural land cover types (grass, shrubs, trees/forest, or bare/snow), followed by crops (Fig. 9c).

Figure 9Characteristics of the imbalanced watersheds identified using IMD-APHRO. The watersheds are sorted in ascending order of catchment area (left to right) in all of the panels. The thick broken vertical line in each plot separates the watersheds of Northern India from those of Peninsular India. (a) Catchment area; (b) maximum elevation (blue line shows the average elevation for the whole of India: 600 m); (c) land-cover and land-use fractions (blue line shows the 50 % fraction for reference); (d) average annual P (blue line shows the average annual P for the whole of India: 1100 mm); (e) maximum runoff coefficient (blue line indicates a value of 1.0); (f) cumulative maximum live storage capacity expressed as a fraction of the annual P; and (g) maximum (shaded bars) and minimum (unshaded bars) ΔTWS as fractions of the annual P.

Download

The average annual P for these imbalanced watersheds is typically around 2000 mm yr⁻¹ (Fig. 9d) – about twice the average annual P for the whole of India (about 1100 mm, see Table S3.1). Thus, such watersheds are typically wetter than the rest of India. Moreover, what is presented here is the observed P, which is potentially affected by UoP, and the actual P could be much higher. The maximum annual runoff coefficient for these watersheds typically exceeds 1 (median value of 1.15; maximum value of 3.33; Fig. 9e). The extent of reservoir storage is quantified as the cumulative sum of the maximum live storage capacity of all reservoirs present in the watershed, expressed as a percentage of the average annual P (Fig. 9f). While most watersheds have relatively minimal storage, some of them could have more than 50 % of the annual P captured in the reservoirs. However, the P data used here is the observed P (affected by UoP), not the actual P. Therefore, the actual effect of reservoirs is expected to be smaller than what is represented here. Finally, the minimum and maximum watershed-averaged values of ΔTWS expressed as fractions of the annual P (Fig. 9g) indicate that the magnitude of ΔTWS is less than 20 % for most of these watersheds.

Based on these physical characteristics, the imbalanced watersheds identified using the IMD-APHRO dataset are typically forested (or minimally impacted by agriculture) and located in relatively wet regions and at relatively high elevations, often have annual runoff coefficients exceeding 1.0, and, in general, are minimally impacted by reservoir storage. Moreover, based on a visual comparison of the extent of large-scale management shown in Fig. 3 and the locations of imbalanced watersheds in Fig. 8, the imbalanced watersheds can be considered to be minimally affected by groundwater extraction and diversions. Furthermore, based on a visual comparison of annual ΔTWS in Fig. 4 and the locations of imbalanced watersheds in Fig. 8, the imbalanced watersheds are typically in regions not affected by relatively large annual changes in TWS.

4.3 UoP within IMD-APHRO versus other datasets

Similar to the earlier analysis in which watersheds potentially affected by UoP were identified using the IMD-APHRO dataset, the potential UoP within other P datasets is analyzed in this section. For each P dataset, Table 3 shows the number of station-years across all imbalanced watersheds according to that dataset. The number of station-years by scenario are tabulated separately for the watersheds of Northern India and Peninsular India. Since the different P datasets have differing time spans, the total number of WYs varies by P dataset. ERA5, IMD-APHRO and TERRA have the longest time spans (782 station-years in Northern India and 6153 station-years in Peninsular India), while SM2RAIN has the shortest time span (195 station-years in Northern India and 1784 station-years in Peninsular India).

Table 3Total number of station-years analyzed and the imbalanced station-years by scenario for all the watersheds of Northern India (left) and Peninsular India (right) for each P dataset.

Download Print Version | Download XLSX

The total number of imbalanced years for which either UoP scenario is realized is expressed as the percentage of the total analyzed station-years. This percentage acts as proxy for the extent of UoP, and can vary from about 2 % to 29 % in Northern India and from 5 % to 19 % in Peninsular India, depending on the P dataset. The APHRO dataset is consistent with IMD-APHRO in Peninsular India but not in Northern India. Across the entire study domain, the satellite-based GSMAP, PERSIANN and CMORPH datasets typically have the highest percentages of imbalanced station-years, while the reanalysis-based datasets of ERA5 and IMDAA have the lowest percentages. While ERA5 and IMDAA are consistent across both Northern and Peninsular India, the MSWEP and TERRA datasets have the lowest percentages in Peninsular India but do not have such low percentages in Northern India. The reanalysis-based datasets of ERA5 and IMDAA outperform IMD-APHRO as well as the high-resolution satellite products such as CMORPH and PERSIANN. The GSMAP dataset has the highest percentage of imbalanced watersheds in both Northern and Peninsular India.

The statistics presented in Table 3 are based on a specific set of heuristics (0.70 and 1.20) used within Scenario II. In order to understand the impacts of changing the heuristics, three other sets of heuristics were tried. Instead of 0.70, values of 0.60, 0.80 and 0.90 were used, and instead of 1.20, values of 1.10, 1.30 and 1.40 were used within Scenario II. Tables S9.1–9.3 show the new set of statistics (similar to Table 3) for each set of heuristics. It is evident from these tables that the performance of the datasets remains similar to that shown in Table 3. ERA5 and IMDAA outperform IMD-APHRO consistently across both Northern and Peninsular India, while the MSWEP and TERRA datasets have the lowest percentages in Peninsular India but do not have such low percentages in Northern India.

The metrics presented in Table 3 are associated with watersheds where adequate hydrometeorological data are available. Since these watersheds are limited to only certain portions of India, these metrics do not accurately reflect the spatial distribution of UoP present within each P dataset. In order to assess the spatial distribution of UoP, potential correction factors (CFs) are estimated for select datasets in Sect. 4.4. The ERA5, IMDAA, MSWEP and TERRA datasets are chosen for further analysis because of their potential ability to be less affected by UoP than IMD-APHRO.

4.4 Potential correction factors (CFs) for specific datasets

Correction factors (CFs) represent ratios of actual and observed P. Since it is not possible to estimate them without knowing the actual P, they were estimated assuming that select datasets from the above analysis are reasonable proxies for the actual P. These estimated CFs are referred to as potential CFs to distinguish them from true CFs. As mentioned in Sect. 4.3, the ERA5, IMDAA, MSWEP and TERRA datasets suffer less from UoP than IMD-APHRO. Using these datasets, potential CFs were estimated using Eq. (7):

For each dataset, data were first aggregated to IMD's resolution of 0.25° (∼25 km). Then, for the 30-year common data period of WY 1985–2014, the gridwise average annual P was estimated. The ratio of gridwise 30-year average annual P between each dataset and IMD is presented in Fig. 10. The spatial domain is limited to the political boundaries of India where IMD data are available.

Figure 10Gridwise potential CFs for select datasets based on the period WY 1985–2014.

The spatial maps of potential CFs shown in Fig. 10 can be compared to those presented in Figs. S1.1 and S1.2. High CFs are present in the mountainous western coast of India for all four datasets and in the mountainous parts of Northern India for only fERA5 and IMDAA. This is consistent with the percentage of imbalanced station-years associated with each of these datasets (see Table 3). Another feature that is evident from Fig. 10 is that the highest CFs occur in the wettest parts of India (Fig. S1.2). If these potential CFs are reasonably accurate, then one could conclude that UoP is a substantial problem in the wettest parts of India. A CF of at least 1.5 (yellow-, green- or blue-shaded areas in Fig. 10) indicates that the actual P is at least 50 % higher than the observed P. There are wide swaths of mountainous and hilly regions of India with such CFs. In order to identify the river basins of India that are most affected by UoP, basin-aggregated potential CFs are analyzed.

Table 4Basin-aggregated potential CFs for select datasets based on the period WY 1985–2014. Average and maximum values for each dataset are shown for each river basin. The upper half of the table shows the basins of Northern India, while the lower half shows the basins of Peninsular India.

Download Print Version | Download XLSX

Table 4 shows the basin-aggregated potential CFs for the above four P datasets. An additional table for all of the P datasets analyzed here is shown in Table S3.1. The potential CFs shown here were estimated as the ratio of annual P for each dataset and IMD. The average and maximum values for the 30-year period of WY 1985–2014 are shown in Table 4. Since IMD is the main P dataset of interest and is limited to the political boundaries of India, only that portion of each river basin falling within India's boundaries was included when estimating these potential CFs.

Across the whole of India, ERA5, IMDAA and MSWEP are on average 9 %, 26 % and 3 % higher than IMD, respectively, while TERRA is 2 % lower than IMD. However, the maximum values indicate that ERA5, IMDAA, MSWEP and TERRA can be up to 19 %, 37 %, 10 % and 6 % higher than IMD, respectively. There is substantial variability across basins and datasets. For instance, for the Brahmaputra Basin, ERA5 and IMDAA are 56 % and 90 % higher than IMD; however, MSWEP and TERRA are 5 % and 9 % lower than IMD. Similarly, for the Ganga Basin, on average, ERA5, IMDAA and MSWEP are 9 %, 36 % and 8 % higher than IMD; however, TERRA is 1 % lower than IMD. Similarly, for the Indus Basin, on average, ERA5 and IMDAA are 6 % and 26 % higher than IMD; however, MSWEP and TERRA are 33 % and 43 % lower than IMD. This pattern of ERA5 and IMDAA being higher than IMD while MSWEP and TERRA were lower than IMD in the basins of Northern India is consistent with the potential CFs shown in Fig. 10. ERA5 and IMDAA have CFs exceeding 1 in many regions of Northern India, while MSWEP and TERRA do not have such high CFs to the same extent.

Table 4 also shows that for most basins of Peninsular India, potential CFs from the four selected P datasets are almost always greater than 1. This implies that P is underestimated across most of Peninsular India, regardless of which of the four datasets is used as a proxy for actual P. The Godavari and Krishna basins are the two largest basins of Peninsular India. In the Godavari Basin, on average, the four datasets are 4 % to 13 % higher than IMD. In the Krishna Basin, on average, the four datasets are 13 % to 19 % higher than IMD. The wettest basins of Peninsular India are the WFR north and WFR south basins. In these two basins, MSWEP and TERRA are higher than IMD, while ERA5 and IMDAA tend to be similar to or lower than IMD. This is consistent with the percentage of imbalanced station-years associated with each of these datasets in Peninsular India (see Table 3).

5.1 Limitations

The watersheds affected by UoP were identified by analyzing the extent of the annual water imbalance. As such, the results are dependent on the quality of the data and strength of the assumptions used. The limitations of the datasets used here and also the limitations imposed by the assumptions made within this analysis are discussed here.

5.2 Spurious patterns within IMD

During the course of this analysis, several potential issues with trends in the IMD dataset were encountered. The following is a discussion on basin-scale trends in the IMD dataset and those present in other datasets. For the purposes of this discussion, the spatial domain is limited to the political boundaries of India where IMD data are available. Basin-scale aggregation of gridded P was performed only using the grids within India's boundaries.

Trends in the four datasets identified in Sect. 4.3 are compared against those in IMD. Trends in basin-aggregated annual P for WY 1985–2014 were estimated using the non-parametric Theil–Sen slope (Helsel et al., 2020), making use of the R statistical package “RobustLinearReg” (Hurtado, 2023). Table 5 shows the trends for select basins where mountains are present. Table S3.2 shows the trends for all of the basins and all of the P datasets.

Table 5Trends in annual basin-aggregated P (mm yr⁻¹) for WY 1985–2014 for select datasets and select river basins. Values that are statistically significant at the 95 % confidence level are indicated by ^*.

Download Print Version | Download XLSX

Figure 11(a) Annual P for Barak Basin for select datasets. IMD (blue) is compared against select datasets. The periods of interest, WY 1961–1980 and WY 1981–2000, are highlighted. Thin lines show the annual values, while the thick lines show the 9-year running average. (b) Same as (a) but for the Indus Basin.

Download

The annual P from IMD for the whole of India shows a decreasing trend of −1.7 mm yr⁻¹. In contrast to IMD's decreasing trend, all other datasets have an increasing trend. However, none of these trends are statistically significant at the 95 % confidence level. There is substantial variation in regional trends, as is evident in the trends presented for individual basins. For the IMD dataset, the Barak, Brahmaputra, Ganga and Indus basins in Northern India show decreasing trends. However, other datasets do not always have the same magnitude or sign as IMD. For instance, for the Ganga Basin, IMD shows a negative trend of −4 mm yr⁻¹, while ERA5, IMDAA and MSWEP show a positive trend. None of these trends are statistically significant at the 95 % confidence level. For the wettest basin of India, the WFR south basin, all datasets show a positive trend, with most of them being statistically significant. Based on Table S3.2, there appears to be more consistency in trends between IMD and other datasets for the basins of Peninsular India compared to the basins of Northern India.

Another issue which was encountered was abrupt changes in the time series of P from the IMD dataset, particularly in the earlier part of its record. The time period of interest here is the 20-year period of WY 1981–2000, which is compared to the prior 20-year period of WY 1961–1980. The time series of basin-averaged annual P from IMD is compared with the corresponding time series from three P datasets which have data available for these periods – APHRO (gauge-based), ERA5 and TERRA (reanalysis-based).

Figure 11 shows such comparisons for the Barak and Indus basins, while Figs. S3.16 to S3.38 show a similar time series comparison for all of the major basins. Both annual values (thin lines) and the 9-year running average (thick lines) are shown in Fig. 11 to highlight the short- and long-term changes in P in each of the datasets. Figure 11 shows that for the Barak Basin, IMD shows an increase in average annual P of about 22 % for WY 1981–2000 relative to WY 1961–1980. However, APHRO, ERA5 and TERRA show a change of 8 %, −3 % and 5 %, respectively. Also, IMD has a distinct visual increasing trend from low values in the early 1960s to high values in the early 1990s. Such a pattern is not present in APHRO, ERA5 or TERRA. Similarly, for the Indus Basin, IMD shows an increase in average annual P of about 35 % for WY 1981–2000 relative to WY 1961–1980. However, APHRO, ERA5 and TERRA show a change of 5 %, 4 % and 3 %, respectively. IMD, once again, has a distinct visual increasing trend from low values in the mid 1970s to high values in the late 1990s. Such a pattern is not present in APHRO, ERA5 or TERRA.

Overall, the above discussion highlights two related issues with the IMD dataset. First, trends present within the IMD dataset are not always present in other datasets. Second, the conspicuous temporal shifts present in the IMD dataset are not present in other datasets. Lin and Huybers (2019) noted a potentially spurious shift in the IMD dataset over central India. It is not known if, and to what extent, such issues are caused by UoP within these datasets.

5.3 Interim measures

Solving the problem of UoP by increasing the station density in relevant areas, by monitoring and analyzing extreme P events and rainfall–runoff relationships for such events, or by any other means requires significant planning and resources from the relevant government agencies. Such efforts are strongly encouraged by the authors. In the interim, there are several useful and feasible ideas the community could pursue to help address the issue of UoP. The following is a brief discussion of such ideas.

Raw station data from the IMD would be extremely helpful for resolving discrepancies associated with trends and discrepancies with other datasets. However, such data are not publicly available. The IMD could help resolve the issue of UoP by making such raw station data publicly available. Other data from India's water agencies, such as streamflow data, which are currently classified by the CWC for Northern India, would also be valuable for addressing UoP.

Some studies have demonstrated the ability of high-resolution simulation models to capture P in watersheds dominated by hilly or mountainous terrain. For instance, Li et al. (2017) implemented the Weather Research and Forecasting (WRF) Hydro model in a high-resolution setting (3 km grid) across a mountainous watershed of Northern India. They demonstrated that such a system can reasonably simulate P and can overcome the deficiencies of typical gauge-based products and satellite-based products. Hunt and Menon (2020) also used the WRF-Hydro modeling system in a high-resolution setting (4 km grid) to analyze P during the catastrophic flooding of 2018 in the state of Kerala in Peninsular India. Their study was also able to capture the spatial structure and magnitude of the observed P reasonably well. Such modeling studies should be pursued further to better identify and quantify UoP within traditional products.