Comparative Study of Two State-of-the-Art Semi-Distributed Hydrological Models

Abstract

Performance of a newly developed semi-distributed (grid-based) hydrological model (satellite-based hydrological model (SHM)) has been compared with another semi-distributed soil and water assessment tool (SWAT)—a widely used hydrological response unit (HRU)-based hydrological model at a large scale (12,900 km²) river basin for monthly streamflow simulation. The grid-based model has a grid cell size of 25 km^2, and the HRU-based model was set with an average HRU area of 25.2 km² to keep a balance between the discretization of the two models. Both the model setups are calibrated against the observed streamflow over the period 1977 to 1990 (with 1976 as the warm-up period) and validated over the period 1991 to 2004 by comparing simulated and observed hydrographs as well as using coefficient of determination (R²), Nash–Sutcliffe efficiency (NSE), and percent bias (PBIAS) as statistical indices. Result of SHM simulation (NSE: 0.92 for calibration period; NSE: 0.92 for validation period) appears to be superior in comparison to SWAT simulation (NSE: 0.72 for calibration period; NSE: 0.50 for validation period) for both calibration and validation periods. The models’ performances are also analyzed for annual peak flow, monthly flow variability, and for different flow percentiles. SHM has performed better in simulating annual peak flows and has reproduced the annual variability of observed streamflow for every month of the year. In addition, SHM estimates normal, moderately high, and high flows better than SWAT. Furthermore, total uncertainties of models’ simulation have been analyzed using quantile regression technique and eventually quantified with scatter plots between P (measured data bracketed by the 95 percent predictive uncertainty (PPU) band) and R (the relative length of the 95PPU band with respect to the model simulated values)-values, for calibration and validation periods, for both the model simulations. The analysis confirms the superiority of SHM over its counterpart. Differences in data interpolation techniques and physical processes of the models are identified as the probable reasons behind the differences among the models’ outputs.

1. Introduction

Distributed hydrological models, with varying degree of complexity, are essential tools for modeling the spatial variability effects of basin characteristics and forcing variables (e.g., precipitation) on streamflow [1,2,3,4]. These models divide spatially heterogeneous space (basin or watersheds) into a number of near homogeneous units following various discretization schemes including: representative elementary area (REA) [5], grouped response unit (GRU) [6], representative elementary watersheds (REW) [7], hydro-landscape unit [8], triangular irregular network (TIN) [9], hydrological response unit (HRU) [10] and grid-based approaches [11]. Among these discretization schemes, HRU and square grid approaches are the most commonly used in hydrological modeling.

HRUs are formed by lumping individual areas of similar soil, topography, and land-use altogether within a sub-basin. However, there is no interaction between the HRUs, and these are routed individually to the sub-basin outlet [12,13]. Arnold et al. [13] studied the effect of HRU discretization on streamflow and concluded that many HRUs are too big to resolve into individual topographic positions since they occupy the landscape continuum from the divide up to valley bottom. They also identified that the impact of an upslope HRU management on a downslope HRU cannot be assessed. Furthermore, though the HRU-based approach is simple and computationally efficient, spatial information from high-resolution land-use or soil maps can be lost depending on the scale of the HRUs. On the other hand, grid-based discretization scheme uses aggregated spatial variations over each grid. The use of smaller HRUs, instead of grid cells, may yield similar results but incorporating raster data into the HRU based approach would require data transformation from simple grid geometry to a patchy geometry of irregular polygons. Therefore, a grid-based approach appears better to use to avoid the inconvenience.

To describe the basin topography accurately, the grid size is considered up to an acceptable range while keeping the trade-off between model simulation time and simulation accuracy to a minimum. Though, in theory, modeling with a finer grid cell resolution is expected to yield better results because of better-resolved model input data (e.g., rainfall, topography, land cover, etc.), it may not always happen [14]. Therefore, several studies have focused on examining the impact of grid cell size on model simulation results and model simulation time to find out the optimum resolution of grid cells for a particular modeling study. Finnerty et al. [15] illustrated the changes in water budget with continuous simulations at various spatial scales, ranging from 4 km × 4 km to 256 km × 256 km. Wood et al. [16] used a 1° × 1° gridded structure for modeling continental-scale basins. Kuo et al. [17] applied a variable-source-area hydrological model to grid sizes ranging from 10 to 600 m and observed increasing misrepresentation of the curvature of the landscape with increasing grid size. In modeling the 375,000 km² Senegal River basin, Andersen et al. [18] used grid cell resolution of 4 km × 4 km. Booij [19] compared three versions of Hydrologiska Byråns Vattenbalansavdelning (HBV) model [20] with different spatial resolutions in the Meuse river basin in Europe and found that the version with finer resolution reproduced a slightly improved average and extreme discharge behavior at the basin outlet in both calibration and validation periods. Recently, Haghnegahdar et al. [21] carried out a modeling study in a 2700 km² area with model grid cells of 15 km × 15 km resolution.

The effect of different spatial discretization schemes on streamflow simulation has been studied by researchers. For example, Abu El-Nasr et al. [22] assessed performances of fully distributed grid-based MIKE Systeme Hydrologique Europeen (SHE) and the semi-distributed HRU-based SWAT and showed that MIKE SHE can predict the overall variation of stream flow slightly better. There are more examples of studies investigated utility of different grid-based models and compared results with an HRU based SWAT model [23,24,25,26,27,28,29,30]. Arnold et al. [13] used a modified SWAT model, with landscape routing method, to compare modeling results, under four discretization methods: lumped, HRU, catena, and grid. The comparison showed that a high-resolution grid approach would include the impact of an upslope grid cell on a downslope grid cell and provide accurate spatial detailed output. Comparing SWAT model performances with HRU and grid-based structures, Pignotti et al. [31] concluded that the grid-based model under predicts streamflow from 5% to 50% with respect to the usual HRU-based model. Surfleet et al. [32] compared two HRU-based models namely the precipitation–runoff modeling system (PRMS) [33] and groundwater and surface-water flow (GSFLOW) with the grid-based variable infiltration capacity (VIC) model for future climate change analysis and concluded that the future changes can quantitatively be attributed not only to the scale of the models but also to the ability of models to represent hydrological processes. Findings of these various studies also pointed out that model simulation results also vary depending on several factors other than the spatial discretization scheme. These factors include the physiographic characteristics of the basin, seasonality of precipitation, season of the year, and dominating runoff producing mechanisms and, thus, emphasize the uncertainty of analysis of model simulation results for successful comparisons of different hydrological models in a particular study (e.g., [34,35]).

Keeping this in mind, this study aims at in-depth inter-comparison of simulation results of two state-of-the-art semi-distributed hydrological models, namely the satellite-based hydrological model (SHM) and soil and water assessment tool (SWAT), under similar discretization scale, and uncertainty related to the simulations [36,37] in a large scale (>1000 km²) [13,36] sub-tropical river basin, namely Baitarani. The idea behind the similar discretization scale is to reduce the effect of different discretization schemes of the two models and analyze the effect of other factors on the streamflow simulation.

The remainder of this paper is organized as follows. The following section presents the description of the study basin and data used in the study. A description of the models along with sensitive parameters employed in the study is provided in Section 3. The methodologies of model setup, calibration, and validation procedure, as well as the consequent data analysis (including uncertainty analysis), are outlined in Section 4. The results are presented and discussed in Section 5. The final section, Section 6 provides conclusions.

2. Study Area and Data

The study has been performed in Baitarani river basin (12,900 km²) in India which is bounded between 20°35′ N to 22°15′ N latitude and 85°10′ E to 87°03′ E longitude (Figure 1). It comes within the sub-tropical monsoon climate zone [38] and receives an annual rainfall of about 1450 mm (Annual Report, 2011-12, 2011). Almost 80% of the annual rainfall occurs during the four months of south-west monsoon season (June to September) that generates heavy flow and creates floods in lower reaches [39]. Daily temperature varies from 5 °C to 47.5 °C. The elevation of the basin ranges from 10 m to 750 m above mean sea level. Soils of this area vary from rich red loamy to gravely detritus.

For a consistent comparison of performances, the same datasets were used in SHM and SWAT models. Daily Rainfall and daily maximum and minimum temperature have been obtained from the India Meteorological Department (IMD), Pune at 1° × 1° resolution. Data have been interpolated to 5 km × 5 km resolution by using bi-linear interpolation technique to use as input into the SHM. Soil and land use land cover (LULC) maps were collected from the Food and Agriculture Organization (FAO) website (http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/) at 1 km × 1 km scale. The digital elevation model (DEM) of 30 m × 30 m resolution was taken from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) website (https://asterweb.jpl.nasa.gov/gdem.asp). All the static information (soil map, LULC map, and DEM) have been resampled into 5 km × 5 km resolution to use in the SHM. The weather database of SWAT is developed using the weather generator (WXGEN) model using the closest station scheme [40]. Observed streamflow data, at Anandpur gauging station (21.21° N, 86.12° E), were collected for the period of 1977 to 2004 from the Central Water Commission (CWC), Bhubaneswar, India.

3. Comparative Discussion on SHM and SWAT

In this section, short descriptions of the SHM and SWAT are provided (Section 3.1 and Section 3.2, respectively). Then the identified sensitive parameters of both the models, which have been used to calibrate the models, are discussed in Section 3.3.

3.1. Description of the SHM

The SHM works on 5 km × 5 km spatial grid resolution and properties at the center of a cell are assumed to be the properties of the cell. SHM has five modules: surface water (SW), forest (F), snowmelt (S), groundwater (GW), and routing (ROU). SHM grid cells corresponding to forest and snow land cover are modeled using the F and S modules, respectively; whereas other grid cells are modeled using the SW module.

In the SW module [41,42], the Soil Conservation Service (SCS) curve number (CN) method [43] is used to estimate the surface runoff along with the Hargreaves method [44] to estimate the potential evapotranspiration (PET). Soil moisture is estimated by using the water balance technique. The soil profile is considered as a single-layered zone of 300 mm, and moisture-holding and moisture-transmitting characteristics of the soil layer and underlying layer are considered to account for the soil moisture. Infiltrated water wets the soil layer, and excess water from the maximum capacity (saturation) contributes after percolation to GW module. The soil moisture is depleted by evapotranspiration, at a potential rate or actual rate, depending on soil moisture condition.

The F module serves, based on water balancing and the dynamics of the subsurface, to provide output in the form of runoff, soil moisture, evapotranspiration, and contribution to groundwater using the technique and parameters stated in [45]. Subsurface is reckoned on having soil matrix and macropores of main bypass and internal catchment types. The main bypass directly contributes to groundwater. Soil matrix is considered of having three layers, which are important with respect to water balance and change in soil moisture. After infiltration, the saturation of three layers gets started from the top in batch, and after complete saturation of the three layers, the excess water goes to groundwater. After a precipitation event, runoff generation occurs according to the antecedent moisture conditions in the subsurface.

The S module determines the snow density from snow albedo [46] for estimating snowmelt depth by using two different algorithms, viz., the temperature index algorithm and radiation-temperature index algorithm. Since the study area does not have any snow land cover; the S module is not considered in this study.

The GW module uses the contribution from SW, F, and S modules and generates baseflow following the water level variation process described in [47]. The resultant baseflow along with the surface runoff generated from other modules is routed up to the outlet as streamflow.

In SHM, a distributed routing technique [41], termed as time-variant spatially distributed direct hydrograph (SDDH) travel time method [48], was adopted. It requires the flow path, which is derived from DEM. The downstream cell, in the direction of the steepest descent, is defined from the DEM by the use of the flow direction geographic information system (GIS) function with a unique connection from each cell to the watershed outlet. This process produces a cell network to present the flow paths. The threshold number of upstream cells is set equal to two (based on trial and error) to delineate the channel network for the watershed. Any cell with a number of upstream draining cells equal to or greater than the threshold value is considered to be a channel cell, whereas others are considered as overland flow cells. The key point of this approach is the travel time estimation. SHM uses MySQL (open source software) as a relational database management system (RDBMS).

3.2. Description of SWAT

SWAT is used for simulation of the water cycle and its corresponding fluxes of energy and matter (e.g., sediment, nutrients, pesticides, and bacteria) as well as the impact of management practices on these fluxes at basin scale [49]. SWAT uses Microsoft Access as RDBMS. SWAT, however, first discretizes the watershed into a network of irregular sub-basins and then divides each sub-basin into HRUs. The model includes components for hydrology, sedimentation, crop growth, nutrients, and agricultural management [11]. A detailed description of all components of the model can be found in Arnold et al. [49] and Neitsch et al. [10].

In the present study, SWAT has been used with the Soil Conservation Service Curve Number (SCS-CN) method as a runoff generation technique along with the Hargreaves method to determine PET. SWAT calculates baseflow contribution to streamflow from groundwater depending on the water balance approach in a shallow aquifer [49]. In SWAT, runoff is first computed separately for each of the HRUs within the sub-basin and then routed through the stream network to obtain the total streamflow for the watershed. Since the study area does not have snow-covered land, the snow-melt runoff simulation procedure of SWAT is not discussed here.

3.3. Sensitive Parameters of Both the Models Used for Calibration

The total number of parameters of the two models varies in number for streamflow analysis. Three parameters of SHM and seven parameters of SWAT have been found sensitive for streamflow simulation (

Table 1. Summary of the sensitive parameters, respective hydrological processes with methodology, theoretical range of their values and information on the spatial variation over the study area, used by SWAT and SHM to simulate hydrological components.

SWAT						SHM
List of Sensitive Parameters	Process	Method	Theoretical Range of Parameter Values	Calibrated Range of Parameter Values	Spatially Varied or Not	List of Sensitive Parameters	Process	Method	Theoretical Range of Parameter Values	Calibrated Range of Parameter Values	Spatially Varied or Not
Cn2	Runoff	SCS-CN [52]	0–100	50–100	Yes	CN	Runoff (SW)	SCS-CN [52]	0–100	30–100	Yes
Gwqmn	Baseflow	Threshold value based contribution from shallow aquifer storage	0–5000 (mm)	0–300	Yes
Gw_delay			0–500 (day)	8	No
Alpha_bf			0–1	0.5	Yes
Ch_N2	River flow Routing	Variable storage/Muskingum	−0.01–0.03	0.03	No	no	Routing (ROU)	SDDH [48]	0.01–0.05	0.01–0.03	Yes
Ch_K2			0–500 (mm/h)	0.45	No	nc			0.01–0.05	0.015	No
Esco	To compensate the soil evaporative demand with the depth of soil layers	Water balance	0–1	0.01–0.95	Yes
	Total	7						3

Enlargements: Cn2 (Curve number), Gwqmn (Threshold water level in shallow aquifer for base flow), Gw_delay (Delay time for aquifer recharge), Alpha_bf (Baseflow recession constant), Ch_N2 (Manning’s n for main channel), Ch_K2 (Effective hydraulic conductivity of soil), Esco (Soil evaporation compensation coefficient), CN (Curve Number), n_o (Manning’s roughness co-efficient for overland flow), n_c (Manning’s roughness co-efficient for channel flow).

), in this study.

During calibration of SHM, parameters of SW and ROU modules have been changed manually (since an auto-calibration option is not available). For this purpose CN, Manning’s roughness coefficient for overland cell (n_o), and Manning’s roughness coefficient for channel cell (n_c) have been used as sensitive parameters [50]. The parameters of the F and GW modules have been set at their default values as recommended by the developers. The theoretical ranges of sensitive parameters are given in Table 1. CN is responsible for runoff generation in the SW module, and n_o and n_c affect the routing procedure of generated runoff and baseflow from a grid cell up to the outlet of a basin. Using calibrated values of the sensitive parameters, SHM simulates monthly streamflow at Anandpur gauging station of Baitarani basin.

For the SWAT model, seven sensitive parameters are identified (Table 1) for model calibration based on the analysis of parameter sensitivity using the Latin hypercube-one factor at a time (LH-OAT) method [51]. Curve number (Cn2) and baseflow recession constant (Alpha_bf) are responsible for runoff generation; delay time for aquifer recharge (Gw_delay) and threshold water level in a shallow aquifer for base flow (Gwqmn) are responsible for baseflow generation, and the soil evaporation compensation coefficient (Esco) is responsible for soil evaporation losses. Manning’s n for the main channel (Ch_N2) and Effective hydraulic conductivity of soil (Ch_K2) are responsible for controlling river flow routing. Table 1 summarizes the sensitive parameters of both the models with corresponding hydrological processes, estimation methodology and their theoretical ranges. The table also focuses on the spatial variability of the sensitive parameters.

4. Methodology

4.1. Model Setup, Calibration, Validation

At first, both the models were setup with the same input data. For SHM setup, the Baitarani basin is represented by 498 grid cells of 25 km². The threshold values of LULC, soil, and slope were taken, respectively, 1%, 1%, and 2% for the development of HRUs in the SWAT model so that the average area of HRUs is around 25 km² and two discretization schemes come in a balanced scale. This assumption led to having 312 sub-basins and 511 HRUs in the Baitarani basin. However, the smallest HRU has an area of 1.9 km², and the largest HRU has an area of 52 km² with an average area of 25.2 km². Both the models were then calibrated (1977–1990) and validated (1991–2004) on a monthly basis. The performance evaluation of both the models has been done by comparing observed and simulated streamflows by using graphical interpretation and statistical indices, namely coefficient of determination (R²), Nash–Sutcliffe efficiency (NSE), and percent bias (PBIAS) for the calibration and validation periods, separately. The used statistical analyses are discussed below.

4.1.1. Nash Sutcliffe Efficiency (NSE)

It is defined as one minus the sum of the absolute squared differences between observed and simulated values normalized by the variance of observed values [53]. It varies from −∞ to 1, 1 being the perfect fit. It is chosen because of its extensive use in the field of hydrology, which facilitates comparison between different studies. However, it is highly sensitive to peak flows resulting in negligence of low flows.

$$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{Q}}_{\mathrm{o}}-{\mathrm{Q}}_{\mathrm{sim}}\right)}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{Q}}_{\mathrm{o}}-\overline{\mathrm{Qo}}\right)}^2}$$

(1)

where $\mathrm{Q}$_o is the observed streamflow; $\mathrm{Q}$_sim is the simulated streamflow; $\overline{\mathrm{Qo} }$ is the average observed streamflow and N is the number of events in the time-series of streamflow.

4.1.2. Coefficient of Determination (R²)

The coefficient of determination (R²) describes the proportion of the total variance in the observed data that can be explained by a model. It ranges from 0 to 1, with higher values indicating better agreement, and is given by:

$${\mathrm{R}}^2={\left[\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{Q}}_{\mathrm{o}}-\overline{\mathrm{Qo}}\right)\left({\mathrm{Q}}_{\mathrm{sim}}-\overline{\mathrm{Qsim}}\right)}{{\left\{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}\left({\mathrm{Q}}_{\mathrm{o}}-\overline{\mathrm{Qo}}\right)\right\}}^{0.5}{\left\{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{Q}}_{\mathrm{sim}}-\overline{\mathrm{Qsim}}\right)\right\}}^{0.5}}\right]}^2$$

(2)

where, $\overline{\mathrm{Qsim} }$ is the average simulated value of streamflow.

4.1.3. Percent Bias (PBIAS)

It measures the average tendency of the simulated data to be larger or smaller than their observed counterparts. Its ideal value is 0. A positive value indicates model underestimation bias and a negative value indicates model overestimation bias.

$$\mathrm{PBIAS}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{Q}}_{\mathrm{o}}-{\mathrm{Q}}_{\mathrm{sim}}\right)\times 100}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{Q}}_{\mathrm{o}}}$$

(3)

4.2. Analysis of Results

After calibration and validation, the model-simulated streamflows were analyzed to compare model performances for both the periods with respect to annual peaks. Then, inter-annual variability of simulations of both the models for each month of the year, for the total period of analysis, were analyzed. Eventually, the capability of both models was compared using the five percentile series derived from observed data. Therefore, to understand the difference in the capability of the models to simulate different streamflow ranges in an improved manner, four percentile points of observed monthly streamflow, S5 (5th percentile), S25 (25th percentile), S75 (75th percentile), and S95 (95th percentile), were used to divide the overall flow range into five percentile series: low flows (<S5: <139.95 m³/s), moderately low flows (S5–S25: 139.95 m³/s to <349.4 m³/s), normal flows (S25–S75: 349.4 m³/s to <6590 m³/s), moderately high flows (S75–S95: 6590 m³/s to <20,760 m³/s), and high flows (≥S95: ≥20,760 m³/s). Finally, uncertainty analysis has been performed of the models.

Uncertainty Analysis

Using quantile regression, a stochastic approach [54], uncertainty from all sources was analyzed, as a whole and for monthly simulation of both the models at Anandpur gauging station for both the calibration and validation periods. The observed, simulated, and residual values of streamflow are linked with the following equation:

$$\mathrm{Q} \left(\mathrm{t}\right)=\hat{\mathrm{Q}}\left(\mathrm{t}\right)+\mathrm{e}\left(\mathrm{t}\right)$$

(4)

where Q(t) is the observed daily streamflow, $\hat{\mathrm{Q}}\left(\mathrm{t}\right)$ is the simulated streamflow, and e(t) is the residual.

The method assumes a functional relationship between residuals and estimates in the Gaussian domain, i.e., normalized quantile streamflow (NQS) and normalized quantile residual (NQR). A linear relation between NQS and NQR was also used in previous studies [55,56]. Hence, NQR may be expressed as:

NQR = a × NQS + b

(5)

Different quantile regression lines may be obtained by minimizing the absolute bias by assigning different weights to positive and negative residuals in the Gaussian domain. Absolute bias can be considered for this purpose as an objective function (OF) which is expressed mathematically as:

$$\mathrm{OF} = \mathrm{Min}{\displaystyle \sum }{\mathsf{\rho }}_{\mathsf{\tau }}\left(\mathrm{mod}\left[\mathrm{NQR}-\left(\mathrm{a} \times \mathrm{NQS}+\mathrm{b}\right)\right]\right)$$

(6)

where $\mathrm{a}$ is the slope, b is the intercept, and ${\mathsf{\rho }}_{\mathsf{\tau }}$ is the quantile regression function which pushes the regression line to the desired location.

To estimate the streamflow corresponding to a given confidence limit, the simulated streamflow is transformed to the Gaussian domain as NQS first, and then, the error in the Gaussian domain, NQR is estimated using the regression line (Equation (5)). The estimated error, NQR is transformed back to the original domain using the pre-estimated mean and standard deviation of the residual. Finally, the estimated residual is added to the daily simulated streamflow to obtain the streamflow which includes uncertainty. Regression lines were used to analyze uncertainty in the simulated streamflow for different confidence intervals. The slope and intercept of these lines are estimated by Equation (6) using the calibration period data. Furthermore, to verify the correctness of error models, the models were applied for both the calibration and validation periods.

Moreover, to have quantitative realization of uncertainty, P and R values have been calculated, and P vs. R plot has been generated for both the calibration and validation periods.

P-value represents the measured data bracketed by the 95 percent predictive uncertainty (PPU) band [57]. P-value has been determined by the following equation:

$$\mathrm{P}-\mathrm{value}=\frac{{\mathrm{q}}_{\mathrm{in}}}{\mathrm{N}}$$

(7)

where, ${\mathrm{q}}_{\mathrm{in}}$ are the total number of observed data points bracketed by the 95PPU band, N is the total number of observed data points.

R-value expresses the relative length of the 95PPU band with respect to the model simulated values [57]. R-value has been determined by the following equation:

$$\mathrm{R}-\mathrm{value}=\frac{{\overline{\mathrm{d}}}_{\mathrm{x}}}{{\mathsf{\sigma }}_{\mathrm{x}}}$$

(8)

where ${\mathsf{\sigma }}_{\mathrm{x}}$ is the standard deviation of the model simulation x. ${\overline{\mathrm{d}}}_{\mathrm{x}}$ is the average distance between the upper and lower limit of the 95PPU band. ${\overline{\mathrm{d}}}_{\mathrm{x}}$ has been calculated using the following equation:

$${\overline{\mathrm{d}}}_{\mathrm{x}}= \frac{1}{\mathrm{k}}{\displaystyle \sum }_{\mathrm{l}=1}^{\mathrm{k}}{\left({\mathrm{q}}_{\mathrm{U}}-{\mathrm{q}}_{\mathrm{L}}\right)}_{\mathrm{l}}$$

(9)

where l is counter, k is the total number of simulated data points for streamflow q, ${\mathrm{q}}_{\mathrm{U}}$ and ${\mathrm{q}}_{\mathrm{L}}$ are the upper and lower limit of the 95PPU band.

Both the values vary between 0 and 1. P-value equal to 1 and R-value 0 represent the best model simulation with no uncertainty. In the P-Q plot, this point can be identified as the point of no uncertainty. Since to reach the point of no uncertainty is nearly impossible to achieve for any model simulation as a result of model uncertainties and measurement errors, the simulation nearest to the point may be considered as the simulation with the lowest uncertainty.

5. Results and Discussion

5.1. Calibration and Validation of the Models

Comparison between observed and models’ simulated monthly streamflow are shown in Figure 2a for the calibration period and in Figure 2b for the validation period. Figures show good agreement among observed and simulated streamflow by both the models. However, SHM simulates the temporal pattern of observed streamflow relatively better in comparison to the SWAT model in both calibration and validation periods including the reproduction of peak flows. To strengthen this observation 1:1 scatter plots, between observed and models’ simulations for the calibration and validation periods, have also been used (shown in Figure 2a,b). From the scatter plots it is evident that SWAT simulated streamflow deviates considerably from the observed streamflow with respect to the SHM simulated counterpart during both the calibration and validation periods. Moreover, scatter plots also depict that SWAT underestimates high flow more in comparison to SHM.

The goodness-of-fit statistics of both the models on monthly calibration and validation are shown in

Table 2. Calibration and validation performances of the models at a monthly scale.

Period	Statistics	SHM	SWAT
Calibration	R2	0.93	0.75
	NSE	0.92	0.72
	PBIAS	11.62	2.01
Validation	R2	0.93	0.58
	NSE	0.92	0.50
	PBIAS	8.67	−1.4

. Generally, if R² > 0.6, NSE > 0.5, and −25% ≤ PBIAS ≤ 25%, the model simulation results are judged as satisfactory [58,59]. Thus, both SHM and SWAT models have produced satisfactory model simulations for both the calibration and validation periods in the study area. However, the monthly streamflow simulated by SHM shows better fit with the observed monthly flow in comparison to the SWAT simulated streamflow during the calibration as well as validation periods. SHM shows similarity in results during both the calibration and validation periods with a slightly reduced PBIAS during validation than calibration period, thus, improvement in water balance dynamics. On the other hand, SWAT shows considerable deterioration in results during the validation period in comparison to the calibration period which is evident from the values of R² and NSE (Table 2). The results, thus, show improved performances of both SHM and SWAT simulations in comparison to previous studies performed at the Anandpur sub-basin [60,61,62,63,64,65].

5.2. Analysis to Compare Annual Peaks

To perform comparison of the annual peak simulation capabilities of both the models, observed and simulated annual peaks (from both the models) for the calibration (Figure 3a) and validation (Figure 3b) periods have been plotted against the 1:1 line. Figure 3 depicts that SHM reproduces annual peaks better than SWAT. Therefore, SHM can be a good option for streamflow simulation for extreme rainfall events as well as analyzing flooding possibility in the region. Findings are well comparable with the study performed by Baratti et al. [66] in which they estimated annual flood frequency for the same region.

5.3. Inter-Annual Variability of Model Simulations

To understand the difference in models’ capabilities of producing inter-annual variability of monthly streamflow, comparison between observed and simulated monthly streamflow from both the models have been analyzed and are shown in Figure 4. From the figure, it is evident that SHM performs satisfactorily in simulating streamflow during the months of June to October (monsoon) season with the best simulation identified for the month of August throughout the analysis period. In addition, it is also evident that SHM reproduces observed streamflow better for all the months over the analysis period in comparison to SWAT streamflow.

The differences in the results of two models for inter-annual variability on a monthly scale for the total period of analysis is mainly attributable to two reasons: different input data interpolation schemes and variation in modeling processes. First, meteorological data have been bi-linearly interpolated into 5 km × 5 km to run SHM and SWAT model and have used the meteorological data from the closest IMD grid to simulate monthly streamflow in a sub-basin instead of interpolation [67] (as stated earlier in Section 2). Different input data interpolation schemes and variation in spatial discretization create the difference in the spatial distribution of meteorological input for the models [30]. Second, apart from SCS-CN of the SW module, no modeling process of SHM matches with the SWAT model. However, the modeling process combination of SHM proved to be better in monsoonal months in comparison to SWAT; though modeling processes of both models may require improvement for the low rainfall months. In particular, water level variation approach of baseflow generation, in SHM, and water balance approach of baseflow calculation, in SWAT, may be compared with separate analysis for the non-monsoonal months, for the purpose. Furthermore, better calibration may also improve results for months with low rainfall.

5.4. Comparison of Model Simulations for Percentile Flows

The models’ performance has also been analyzed in simulating streamflow of various magnitudes by considering five percentile classes (Section 4.2). The respective simulated flows by both models have been compared with that of observed streamflow by using scatter plots for the calibration (Figure 5a) and validation periods (Figure 5b). The performance of SWAT in simulating moderately low flows during the calibration period is better than simulating other streamflow percentiles during both the calibration and validation periods. SHM performs better for simulating normal, moderately high, and high flows during both the calibration and validation periods. Overall, both the models show an extremely poor performance in simulating low flows during both the periods and moderately low flows during the validation period.

The variation in percentile flow estimation of the models can also be attributable to different input data interpolation schemes. However, SCS-CN plays a major role in both models. Therefore, streamflow simulation may not be appropriate when the rainfall amount is small [30,49,68]. Similar results have been identified for the non-monsoonal months during analysis of inter-annual variability of the models, in the previous section. In addition, the different runoff generation technique of the F module and baseflow generation technique of the GW module of SHM (stated earlier in Section 3.1) with respect to techniques used in the SWAT model are also responsible for the different results of the months. In particular, the soil matrix and antecedent condition of the F module may play a role in the poor model simulation of SHM for low and moderate low flows. Moreover, the routing technique of SHM seems to be the reason behind the upper hand in simulating high flows in comparison to SWAT, by capturing the travel time of the streamflow in a better manner.

5.5. Uncertainty Analysis of Monthly Simulations

Figure 6a,c present the 95PPU uncertainty band for monthly simulation during the calibration period. Figure 6e,g present the 95PPU uncertainty band for monthly simulation during the validation period. Among them Figure 6a,e are for SHM simulations and Figure 6c,g are for SWAT simulations. In addition, Figure 6b,d present the scatter plot of NQR and NQS along with two regression lines: corresponding to upper and lower limits of 95% confidence interval (CI) and one corresponding to the median for the calibration period. Figure 6b,d are for the SHM and SWAT simulation, respectively. Figure 6f,h present the scatter plot of NQR and NQS along with two regression lines: corresponding to upper and lower limits of 95% confidence interval (CI) and one corresponding to the median for the validation period. Figure 6f,h are for SHM and SWAT simulation, respectively.

Figure 6a,c,e,g depict that most of the observed streamflow falls inside the defined bands, though the amount is higher for SHM simulations (Figure 6a,c). Moreover, from Figure 6a,c,e,g it is also evident that the width of 95PPU band is thinner for SHM simulations in comparison to SWAT simulations. Thus, it can be inferred that SHM has less uncertainty in model simulations in comparison to SWAT simulations. Figure 6b,d,f,h depict the relationship between residual and simulated streamflow in the Gaussian domain and confirm that the simulated streamflow is able to capture 95% of the observed streamflow during the calibration and validation periods for both the models. For estimating the collective uncertainty, Dogulu et al. [69] supported the use of the quantile regression (QR) technique due to its simplicity and linearity which has been used elaborately by Kumar et al. [70].

P and R values based on the uncertainty analysis results of SHM and SWAT simulations for the calibration and validation periods, respectively shown in Figure 7a,b, elaborate that SHM poses less uncertainty in monthly simulation than SWAT model.

Though uncertainties from all sources have been counted in the QR uncertainty analysis technique, spatial distributions of input data are different for the models due to different data interpolation techniques and model structures of the two models. Although the models’ parameters take care of the modeling processes during calibration, the spatial variation of input data may affect the uncertainty of the models’ simulation significantly. The results of the uncertainty analysis also represent this aspect and show that SHM represents the spatial variations of landscape characteristics and input data more accurately.

6. Conclusions

SHM and SWAT models were used to simulate the monthly streamflow at Anandpur gauging station of Baitarani basin for in-depth inter-comparison of the models’ performances. The SWAT model was set to have an average size of the HRUs equal to 25.2 km², (nearly equal to the grid cell resolution of SHM, i.e., 25 km²) so that the two discretization schemes were in similar scale. Results showed that although both SHM and SWAT have produced reasonable results, SHM performed better. To be more specific, SHM performed better in simulating annual peak flows, and reproduced the annual variability of observed streamflow for every month of the year. In addition, SHM estimates normal, moderately high, and high flows better than SWAT. Uncertainty analysis of simulated streamflow of both the models also supports the superiority of SHM model in comparison with SWAT model. Possible impacts of the model structure were also identified for the results.

In summary, SHM produced better results in comparison to SWAT at the monthly scale with proof of better model structure for the large research catchment. However, we cannot draw a conclusion that grid-based hydrological modeling is better than the HRU based. More researches should be carried out for comparing different discretization schemes for other Indian basins and other parts of the world.