Comment on hess-2021-306 Anonymous Referee # 1 Referee comment on " Quantifying the impact of land cover changes on hydrological responses in India

This is my second review of the study by Naha et al in which the authors attempt to quantify the changes in the hydrological cycle due to changes in the land cover under future climate change scenarios. I am pleased to see that the authors have taken onboard all my suggestions/criticisms of the previous version. Overall, the manuscript has significantly improved and the results and analysis are appropriate for a regional study of land cover impacts on hydrology. I especially appreciate the inclusion of a comprehensive sensitivity analysis of the model parameters. I only have a few minor comments which are easily addressable.

with different land use scenarios. However, these models through their complex interactions among the model parameters to generate hydrological processes, can introduce significant uncertainties to the hydrological projections. Therefore, we seek to further understand the uncertainties associated with model parameterization in those simulated hydrological responses due to different land cover scenarios. We performed a sensitivity-guided model 15 calibration of a physically semi-distributed model, the Variable Infiltration Capacity (VIC) within a Monte Carlo Framework to generate behavioural models for subcatchments of the Mahanadi river basin. These behavioural models are then used in conjunction with historical and future land cover scenarios from the recently released, Land use Harmonisation (LUH2) to generate hydrological predictions and related uncertainties from behavioural model 20 parameterisation. The LUH2 dataset indicates a noticeable increase in the cropland (23.3% cover) at the expense of forest (22.65% cover) by the end of year 2100 compared to the baseline year, 2005. As a response, simulation results indicate a median percent increase in the extreme flows (defined as the 95th percentile or higher river flow magnitude) and mean annual flows in the range of 1.8 to 11.3% across the subcatchments. The direct conversion of 25 forested areas to agriculture (on the order of 30,000 km 2 ) reduces the Leaf Area Index and which subsequently reduces the Evapotranspiration (ET) and increases surface runoff.

Context and Background
Land use and land cover change (LULC) induced by the rapid anthropogenic activities, is one of the major causes of change in hydrological and watershed processes (Rogger et al., 2016).
Alterations of existing land cover types and land management practices in a catchment can 40 thereby, significantly modify the rainfall path into runoff by changing the hydrological dynamics such as surface runoff, baseflow, Evapotranspiration (ET), water holding capacity of the soil, interception and groundwater recharge, thus reflecting a change in the water demand (Berihun et al., 2019;Bosch and Hewlett, 1982;Costa et al., 2003;Foley et al., 2005;Garg et al., 2017;Hamman et al., 2018;Mao and Cherkauer, 2009;Rogger et al., 2016;Zhang 45 et al., 2014). For instance, developing countries like India are facing rapid growth in population which has prominent effects on LULC dynamics through deforestation, rapid urbanization, and agricultural intensification, subsequently modifying the hydrological cycle in many river basins of India. A recent analysis on global land cover change have shown a significant increase of 82% in the croplands in India (Chen et al., 2019a;IPCC, 2019). 50 Therefore, a comprehensive understanding and evaluation of land cover change impacts on hydrological processes are essential for decision makers to plan environmental policies which focuses on water resource allocations, riparian ecosystem protection and river restoration (Chen et al., 2019b;Chu et al., 2013).
Many studies have attempted to evaluate the hydrological responses to different LULC 55 patterns on specific geographic locations (Abe et al., 2018;Chu et al., 2013;Eum et al., 2016;Li et al., 2015;Ma et al., 2010;Rodriguez and Tomasella, 2016;Viola et al., 2014;Woldesenbet et al., 2017) including Indian river basins (Babar and Ramesh, 2015;Dadhwal et al., 2010;Das et al., 2018;Gebremicael et al., 2019;Wilk and Hughes, 2002). Most of these studies used physically distributed hydrological models (e.g., SWAT, VIC, MIKE-SHE) to simulate the 60 complex hydrological processes and to examine the impact of LULC changes on those processes. Conventionally, this is done by calibrating and validating the hydrological model against the observed data and then setting up that single calibrated model for a baseline land https://doi.org/10.5194/hess-2021-306 Preprint. Discussion started: 23 June 2021 c Author(s) 2021. CC BY 4.0 License. cover scenario. The calibrated model is then run for different land use scenarios and subsequently the differences in simulations are compared. However, it is widely recognised 65 that hydrological predictions obtained from single calibrated model can be biased and therefore the measure of their reliability is always questionable (Beven and Binley, 1992;Huang and Liang, 2006). There may exist 'equally probable parameter set' or 'behavioural set' that can yield equally good or acceptable model predictions, due to the complex interactions among the model parameters to represent the complex hydrological processes. This is known 70 as equifinality and is considered as one of the main sources of uncertainty in hydrological modelling (Her et al., 2019). Recent climate change studies have acknowledged the uncertainties stemming from model parameters, and therefore they take into account these uncertainties while predicting the hydrological responses due to climate change (Chaney et al., 2015;Feng and Beighley, 2020;Her et al., 2019;Huang and Liang, 2006;Joseph et al., 75 2018; Mockler et al., 2016;Singh et al., 2014). However, little is known about the contributions of model parameter uncertainties to the land use change impacts and thus, very few studies exist (Breuer et al., 2006;Chen et al., 2019b) which reported that uncertainties associated with the model parameters could significantly influence land cover change impacts and hence should not be overlooked while modelling hydrologic responses to LULC change. 80 This paper specifically focusses on the Mahanadi River basin, an easterly flowing river basin in India. Eastern part of India is amongst the most rapidly changing landscape over the country, specifically, Mahanadi river basin has undergone drastic land cover changes in the last decades (Behera et al., 2018;Dadhwal et al., 2010). In this study, we address the science  (Liang et al., 1994) and historical and future land cover scenarios from the Land Use Harmonisation 2 (LUH2) database (Hurtt et al., 2011) are used to simulate the discharge and other hydrological components at daily time scales in the Mahanadi river basin. The ability of VIC to simulate the impacts of LULC changes on hydrology are well documented in various research articles (Garg et al., 2017(Garg et al., , 2019Hurkmans et al., 2009;Mao 100 and Cherkauer, 2009;Patidar and Behera, 2019;Zhang et al., 2014).
We first perform sensitivity analysis of the model parameters and calibrate the hydrological model using Monte Carlo simulations to identify behavioural model simulations that implicitly account for the uncertainties from model parameterisation. Behavioural models are then used to predict the hydrological impacts due to different LULC scenarios. The land cover 105 scenarios used in this study are most up-to-date scenarios which represents future changes in the LULC based on Shared Socioeconomic Pathways (SSPs). Previous studies (Breuer et al., 2006;Chen et al., 2019b) have focussed only on the historical land use scenarios to evaluate the hydrological impacts, however and to our knowledge , this is the first study that uses applications of the VIC model in conjunction with future land cover datasets produced under 110 combined SSP and RCP scenarios. While most past studies in other catchments used aggregated (monthly) time steps to model the change, we use daily time steps to capture the dynamics of daily flow variability. Moreover, analysis in most land use impact studies is limited just with the streamflow, missing an overall picture of the hydrological processes. The Mahanadi River basin is located in the eastern part of India ( Figure 1) and drains an area of 141,589 km 2 , which nearly accounts for 4.3% of the total geographical area of India. The 120 basin has a varying topography with its lowest elevated area (-17 m) lying in the coastal reaches and the highest elevated area (1323 m) in the northern hills. The basin is characterized by tropical climate zone and receives rainfall from southwest monsoons which commences in June and lasts till October. The average annual rainfall is 1200 mm, with 90% of the total annual rainfall occurring during the monsoon months (Jin et al., 2018). The mean 125 annual discharge is 1895 m 3 /s. The basin is also subjected to spatial variability in terms of https://doi.org/10.5194/hess-2021-306 Preprint. Discussion started: 23 June 2021 c Author(s) 2021. CC BY 4.0 License.
receiving rainfall which has resulted in floods in some parts of the basin and drought in others.
Notice that about 65% of the basin is placed upstream of the Hirakud dam. The Hirakud dam with a gross storage capacity of 8.136 km 3 is the major hydro project in the river basin insignificant. In addition, loamy and clayey are the major soil types covering roughly 53% and 42% respectively of the total basin area (NBBSS-LUP, India). Approximately 90% of the basin has moderately shallow to deep soil having depths greater than 50 cm. The VIC model is a semi-distributed, land surface hydrologic model which solves both water and energy balance within the grid cells (Cherkauer and Lettenmaier, 1999). VIC maintains sub-grid heterogeneity in land cover classes i.e. divides each grid into tiles based on number of land cover classes, and also considers sub-grid variability in the soil moisture storage capacity (Liang et al., 1994). Surface runoff in VIC is generated through an infiltration excess 155 by using the Xiangjiang formulation (Zhao et al., 1980) in the upper two soil layers. Baseflow is generated from the third soil layer by applying the Arno formulation (Franchini and Pacciani, 1991). Actual Evapotranspiration of each grid cell in VIC is obtained by summing up three types of evaporation: Evaporation from bare soil, evaporation from canopy layer for each To obtain the discharge at the basin outlet, the VIC model is coupled to a stand-alone routing model (Lohmann et al., 1996). This routing model follows a simple river routing scheme where 165 runoff and baseflow are first routed to the edge of the grid cells using an instantaneous unit hydrograph and finally transported to the river/channel network using a linearized St.
Venant's equation. More details about the structure and formulations of the model can be found in the literature (Gao et al., 2010;Liang et al., 1994).
In this study, we implement VIC model, version 4.2.d in the water balance mode at a daily

Datasets
The key input data required by the VIC model are meteorological forcings (precipitation, maximum temperature, minimum temperature and wind speed), soil type, land cover information and topographic features. Topographical features are determined using the 30- The observed discharge at daily scales at multiple gauges (Fig 1) for the simulated time  2010) are obtained from the Central Water Commission (CWC), India, for validating the simulated discharge.

Model parameters
We have selected 16 VIC model parameters (Table 1) Table 1). Typical calibration in VIC involves only streamflow related parameters as also recommended by VIC model developers (Gao et al., 2010;Gou et al., 2020;Xie et al., 2007). However a few studies have reported that some vegetation parameters are sensitive to the runoff in the VIC model (Demaria et al., 2007;Joseph et al., 2018). Parameters subjected to SA in this study include, among others, rarely 205 implemented soil properties: Bulk Density (BD), Fractional water content at wilting point (Wpf) and at critical point (Wcrf); vegetation properties: architectural resistance (rarc) Stomatal Resistance (rmin); and routing parameters: velocity (v) and diffusion (diff). A multiplier of Wcrf is used to compute wpf, to meet the criteria that soil moisture at wilting point should always be less than soil moisture at critical point and the multiplier is tested for sensitivity rather 210 than the actual parameter. Similar approach is followed by Rosolem et al., (2012) while testing sensitivity of parameters in a land surface model. Feasible ranges (minimum and maximum values) of soil parameters (BD, Wcrf, Ksat, Exp) are obtained based on average hydraulic properties of USDA soil textural classes (Cosby et al., 1984;Rawls et al., 1998;Reynolds et al., 2000) considering only the dominant soil textures within the basin. Ranges for the rest of the In addition, the LAI is an important vegetation factor, having substantial control over the water balance by directly influencing the ET rates (Gao et al., 2010;Matheussen et al., 2002).
LAI is specified at a mean monthly basis in VIC. We compared the monthly mean LAI averaged over the time period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) from MODIS AQUA/TERRA with the LAI values from GLDAS 225 database for the river basin. We observed that the monthly mean LAI of all the LULC types from MODIS captures the phenological characteristics more realistically than the GLDAS LAI ( Fig 2b) which shall have further implications on water balance. We find the range of MODIS LAI obtained for each LULC type are well in agreement with the LAI values obtained in a nearby Indian river basin (Patidar and Behera, 2019). 230 Another important factor linking vegetation characteristics to hydrological processes in VIC is the root zone distribution. Typically, root zone allocation in VIC requires user-defined rootzone depths and fractions for each land cover types that are kept fixed during the calibration process. We derived root zone depths and estimated the fractions of roots in each zone following (Zeng, 2002) for each vegetation type, and used a simplified approach to vary the 235 root zone distributions with respect to the soil depths during calibration. This ensures root zone properties vary for different model calibration with a reduced number of parameters, hence providing a more manageable calibration strategy. For details on our root allocation approach, please refer to the Supplement (Sect. S1) 3.4 Experimental design Where r is the linear correlation between observed and simulated discharge, is an estimate of flow variability error and β is a bias term. and are standard deviations in simulated and observed discharge, respectively. µ and µ are mean of simulated and observed discharge, respectively. 260 We first visually inspect SA results and assume a screening threshold value for the sensitivity index, below which the parameters can be regarded as either completely insensitive or less influential. This is a common practice followed in previous SA studies (Gou et al., 2020;Sarrazin et al., 2016;Tang et al., 2007;Vanrolleghem et al., 2015). Next, to achieve a more objective screening convergence result, we compute the width of the 95% confidence interval 265 of the sensitivity indices (Herman et al., 2013;Wang and Solomatine, 2019) and then use maximum width of the 95% confidence interval, as a statistic (Sarrazin et al., 2016) of the total simulations (Chaney et al., 2015;Mockler et al., 2016). This is relevant in our study as choosing model simulations based on a particular KGE score is subjective given that the 285 behavioural performance, as well as the behavioural parameters, vary across the subcatchments. Therefore, we first assess the performance of top 10%, 5% and 2% of model  We used the behavioural models to simulate discharge for the baseline scenario using land cover map from LUH2 of year 2005 so as to attain more confidence in the future scenarios.
We compare LULC maps, NRSC2005 and LUH2005 ( Figure 3) and observe spatial patterns of the most dominant land-use classes, Cropland (CL) and Forest (F), shows a similar spatial distribution and having comparable aerial coverage. The only notable difference in both maps is that the Barren ground (BG) class is missing in LUH2005. Table 2 shows the percentage of area covered by each land use classes in the basin. Note that we will refer to DBF as Forest (F) henceforth. Therefore, any changes observed in the hydrological components in these scenarios will be only due to the change in land use. The percent areas covered by each land use classes at all 340 subcatchments across the scenarios are shown in Table 4.

Sensitivity Analysis, Model Calibration and Validation
It is to be noted that SA is conducted for all subbasins individually, hence the Morris screening results obtained for each subbasin are independent of each other. However, we observe that 345 the non-influential parameters match closely with each other across subbasins (Fig S2).   This underestimation can be attributed to the absence of 12% Barren Ground in the baseline 400 land cover, which is replaced by croplands (4%), forests (5.02%), grasslands (4.57). The increase in flows due to the increase in cropland is compensated by the decrease in flows due to the increase in forest. Therefore, the underestimation in the simulated flows using LUH2005 may result from the increasing grasslands which increased LAI, thus resulting in an increase in ET and decrease in surface runoff, respectively. Contrarily, a slight positive bias of 405 3% is observed at the smallest subcatchment (Sa) in the baseline simulation, compared to +55% in the validation simulation. KGE values obtained across calibration, validation and baseline period indicates an overall good performance of the basin as per the existing studies using KGE as a performance metric (Knoben et al., 2019). In overall, baseline land cover map,   (Table 4). Among all the scenarios, maximum uncertainty is 450 observed in the hypothetical 'All Forest' scenario followed by 'All Cropland' scenario. In overall the uncertainty of hydrological model parameterization is observed at the largest subcatchment Basantpur and decrease with respect to the decrease in the catchment size.

LULC impacts and uncertainties
We analysed the water balance components to understand the factors causing changes in the streamflow. We notice that the model is able to estimate all the water budget components 455 and maintain proper closure of the water balance in all the scenarios across the subcatchments. In overall, we found that the increase in the mean annual flows is caused due to the increment in runoff and reduction in ET across all subcatchments. Positive median changes are observed in runoff in scenarios (NF, FF and CL) ranging between (2.8 to 14) % and negative changes of (-4 to -37) % in F scenario. Negative median changes are observed in ET 460 in scenarios (NF, FF and CL) ranging between (-1.4 to -3.4) % and positive changes of (1.9 to 7.8) % in F scenario. Removal of forests at the expense of cropland decreases the LAI of the natural vegetation and hence decreases ET. Moreover, the removal of forest cover reduces the root water uptake by plants which increases the water content of the second and third layer of the soil. The top thin soil layer in VIC model helps in partitioning the rainfall amount 465 into direct runoff and the amount entering the soil. Therefore, the increase in the cropland results in more direct runoff thus reducing the soil moisture content in the first soil layer. The increase in runoff is not significant, despite the occurrence of major deforestation in the future scenarios. This is because the decrease in ET due to forests removal is compensated as increment in croplands also leads to a major increase in ET rates, which is why we do not see  (Dadhwal et al., 2010), neighbouring basins (Das et al., 2018;Kundu et al., 2017) and elsewhere (Abe et al., 2018;Berihun et al., 2019;Cornelissen et al., 2013;Costa et al., 2003). Kundu et al., (2017) found an increase in runoff and decrease in ET 510 due to the expansion in projected agricultural land in Narmada river basin in India. Das et al., (2018) predicted that deforestation, urbanization and cropland expansion in eastern river basins of India, in the future would increase runoff and baseflow and decrease ET%. We found a small change in mean annual discharge as well as in water balance components despite a major change in land cover. This correlates with research (Ashagrie et al., 2006;Fohrer et al., 515 2001;Hurkmans et al., 2009;Kumar et al., 2018;Patidar and Behera, 2019;Rogger et al., 2016;Viglione et al., 2016;Wagner et al., 2013;Wilk and Hughes, 2002) wherein they have reported that the impacts of land cover change on water balance components in a large-scale river basin are too small to be detected due to the compensation effects. Wilk and Hughes, (2002) showed that removal of large forests led to little or no changes in annual runoff in large 520 heterogeneous catchments in South India. Patidar and Behera, (2019) in a recent study in a large river basin in India, reported that the conversion of forest to agriculture may not alter the water balance significantly as the impacts on ET and runoff cancels out at the basin scale.
The range of these hydrological estimates (Figures 7, 8 and Table 4)

865
Parameter names in bold are sampled on log domain. "a" indicates parameters that are suggested by VIC model developers as the most sensitive parameters (Gao et al., 2010). "b" indicates parameters suggested in the literatures to be tested for sensitivity (Demaria et al., 2007;Gou et al., 2020;Joseph et al., 2018;Yanto et al., 2017 Table 4 Ranges of percent change, change in flows, and uncertainty (i.e., difference between max. and min. predicted flow) in extreme and mean annual flows in all the scenarios with respect to the baseline scenario.