Hydrological model parameterization using NDVI values to account for the effects of land cover change on the rainfall–runoff response

The impact of land use and land cover change on runoff generated from a river basin is of major interest as the tendency towards increased change continues. Studies have shown that land use change can have a substantial effect on the hydrological response of a watershed to rainfall, resulting in quicker response times (Huang et al. 2008), greater river flow volumes (Hawley & Bledsoe 2011), and higher recurrence of floods (Braud et al. 2013).

Conceptual unit hydrograph (UH)-based models, derived from linear systems theory, have been used to investigate the effects of land use change on rainfall–runoff process (e.g., see Kang et al. 1998; Cheng & Wang 2002; Huang et al. 2008). Cheng et al. (2010) explored the effects of population density and impervious areas on the Nash model (Nash 1957) parameters for the urbanized Wu-Tu watershed. Subsequently, Huang et al. (2012) derived a power law between the Nash model parameters and population and impervious areas in the Wu-Tu watershed.

Unlike lumped models, semi-distributed geomorphologic models require large amounts of spatial data in which a combination of remote sensing, geographic information systems (GIS), and hydrological concepts can be used to represent spatial hydro-geomorphological data and information at landscape to global scales (Pickup 1995). Since geomorphologic models use the physical properties of the watershed to inform specification of values for the model's parameters, they are less dependent on the availability of adequate historical event data and can therefore be applied in ungauged watersheds with only limited rainfall–runoff records.

In support of geomorphological modeling, remotely sensed data about land cover type are widely used to derive input variables for a variety of hydrologic-response environmental models (Miller et al. 2007). Several studies assessed the effects of land use changes urbanization on rainfall–runoff response using a modeling strategy supported by remote sensing (e.g., see Verbeiren et al. 2013; Miller et al. 2014; Nourani et al. 2015).

As stated by Yang et al. (2011), during the past two decades remote sensing indices have been proposed to represent land cover types of vegetation (Soil-Adjusted Vegetation Index), bare soil (Normalized Multi-band Drought Index), impervious (Normalized Difference Built-up Index), and water area (Normalized Difference Water Index). However, well-known normalized difference vegetation index (NDVI) and leaf area index (LAI) are the two most commonly used indices to monitor and study vegetation responses to climatic changes at different scales (e.g., Wang et al. 2003; Beck et al. 2006; Eckert et al. 2015; Jiang et al. 2015). NDVI is related to the LAI and vegetation coverage (Baret & Guyot 1991; Turner et al. 1999; Goswami et al. 2015); the higher value of NDVI indicates the larger value of the LAI, and thus the higher amount of vegetation coverage. Li et al. (2009) incorporated the remotely sensed (MODIS-LAI) data into Xinanjiang rainfall–runoff model and assessed the model performance on 210 catchments in south-east Australia. The results showed that the inclusion of LAI data improved both the model calibration results as well as the daily runoff prediction in ungauged catchments. Ali et al. (2011) combined an empirical land use change model and an event scale rainfall–runoff model to quantify the impacts of potential land use change on storm–runoff generation in the Lai Nullah Basin. They calibrated the hydrologic engineering center-hydrologic modeling system (HEC-HMS) rainfall–runoff model and validated for five storm events in the study area. Their results showed a good consistency between the simulated and measured hydrographs at the outlet of the basin. Lakshmi et al. (2011) studied the influence of land surface on hydrometeorology and ecology, addressing various relational investigations among vegetation properties such as NDVI, LAI, surface temperature, and vegetation water content derived from satellite sensors. Also, they examined the effects of vegetation and its relationship with soil moisture on the simulated land–atmospheric interactions through the weather research and forecasting model. Li et al. (2012) used the Xinanjiang-ET and SIMplified HYDrology model (SIMHYD-ET) rainfall–runoff models (Zhang et al. 2009) to investigate impacts of vegetation change and climate variability on streamflow in a Southern Australian catchment. The models incorporated remotely sensed LAI data. The results of their study suggested that increase in plantations can reduce streamflow substantially, even more than climate variability. Liu et al. (2012) analyzed a correlation, based on different vegetation types, between time series of monthly NDVI and Palmer drought severity index (PDSI) during the growing season from April to October in Laohahe catchment. Their results showed that NDVI and PDSI had a good correlation, especially for shrub and grass. The highest value of correlation coefficient was in June when the vegetation is growing and the lower correlation occurred at the end of the growing season for all vegetation types. Zhou et al. (2013) investigated a modified version of Xinanjiang rainfall–runoff model incorporating dynamic remote sensing data for south-east Australian catchments. The results showed that using remote sensing LAI and albedo data in the modified Xinanjiang model led to model performance improvement in wetter bushfire impacted catchments. However, use of vegetation dynamics did not improve runoff time series simulation in a dry catchment. Törnros & Menzel (2014) derived the LAI using the NDVI data for the years 1982–2004. They also made a correlation analysis between precipitation and LAI for several land uses and each month of the year. Based on their results, LAI could be simulated as a function of precipitation. Tesemma et al. (2015a) assessed the effect of observed monthly LAI on hydrological model performance and the simulation of runoff using the variable infiltration capacity (VIC) hydrological model in the Goulburn–Broken catchment of Australia. VIC was calibrated with both observed monthly LAI and long-term mean monthly LAI, which were derived from the global land surface satellite (GLASS) LAI dataset covering the period from 1982 to 2012. Tesemma et al. (2015b) combined a nonlinear model for estimating changes in LAI due to climatic fluctuations with the VIC hydrological model to improve catchment streamflow prediction under a changing climate. The combined model was applied to 13 gauged sub-catchments with different land cover types (crop, pasture, and tree) in the Goulburn–Broken catchment, Australia. Incorporating climate-induced changes in LAI in the VIC model reduced the projected declines in streamflow and confirmed the importance of including the effects of changes in LAI in future projections of streamflow. Zhiqiang et al. (2016) classified the vegetation communities in Poyang Lake wetland using NDVI and linked it with the water regimes by a Gaussian regression model. Tian et al. (2017) examined the effects of revegetation on soil moisture using a biophysically based ecohydrological model. Vegetation covers were investigated based on four types including trees, shrub, grass, and crop. Results showed that trees consume more water than shrub and grass and therefore would result in more soil desiccation.

One of the major effects of urbanization is land cover degradation, therefore, NDVI as a remote sensing index can also be effectively used to monitor urbanization and land cover changes, especially for the arid condition of Iran, where natural watersheds are gradually changing to urbanized areas.

Whereas it is common to consider changes in storm runoff volume due to increases in impervious surfaces (Boyd et al. 1993; Shuster et al. 2005), our investigation of the literature indicates that the relations between temporal changes of NDVI and the conceptual parameters of a rainfall–runoff model are only rarely investigated. In the aforementioned research (e.g., Cheng & Wang 2002; Huang et al. 2008, 2012; Cheng et al. 2010), the effects of urbanization on parameters of the Nash model were investigated using population growth and impervious area data as indices of urbanization. This ignores other causes of increased urbanization such as industrial development. One could mention that impervious surface cover mapping is not straightforward due to heterogeneity of urban areas and possible confusion with bare soils, therefore, vegetation mapping and monitoring is easier. We examine the proposed GUHCR (geomorphological unit hydrograph based on cascade of linear reservoirs) model where the geomorphological properties of the watershed can compensate for a lack of observed event data.

Moreover, while most classical rainfall–runoff models use the historical event data to estimate the model parameters considering the average value for predictions, the effect of climate and anthropogenic changes during the time, which would lead to significant changes in the model parameters, are taken into account for detecting the temporal changes of NDVI using remotely sensed data.

The main aim of this study was to explore the link between the temporal change of NDVI and that of the parameters of the GUHCR model by looking into the consequent impact of land cover changes on flood behavior of the Hervi and Lighvan sub-watersheds in East Azerbaijan, Iran. In this regard, we tried to: (1) compute the NDVIs and determine the land cover changes from 1998 to 2012 using Landsat satellite imagery; (2) propose the rainfall–runoff model and calibrate the model parameter using available storm events within the same period; and (3) relate the change of the model parameter to NDVI values during the study period. These objectives provide the structural sub-headings used in the following Methodology, Results and Discussion sections.

The study area comprises a 77-km² sub-watershed of Lighvan and an approximately 60-km² sub-watershed of Hervi located in north west Iran. Elevation ranges from 1,920 m to 3,453 m above sea level in Hervi and from 2,180 m to 3,453 m in Lighvan. The annual average rainfall in these sub-watersheds is 461 mm, mainly occurring in May to July. The average of mean daily temperature is 12.6 °C and the annual pan evaporation is equal to 1,866 mm. All mentioned values are the long-term mean annual values for the study area. Vegetation in the study area consists of grassland, farmland, and gardens. Although both sub-watersheds contain residential areas, the Hervi sub-watershed is more highly developed area (14.5% versus 10.1% for Lighvan). For the study area, Figure 1 shows an aerial photograph, the location map, and the shuttle radar topography mission (SRTM) digital elevation model (DEM) retrieved via http://data.cgiar-csi.org/srtm/tiles/.

Due to the degree of land use change, this area was selected to address the primary objective of investigating the relationship between NDVI change and model parameters. Table 1 shows the land cover properties of Hervi and Lighvan sub-watersheds in 1976 and 2010. As can be seen in Table 1, dense pasture is the dominant land cover in both Hervi and Lighvan sub-watersheds. The percentages of different land cover have experienced different changes. On the one hand, the percentage of dense pasture, characterized by high NDVI values, showed a considerable decrease from 1976 to 2010. The percentage of sparse pasture, on the other hand, increased considerably in both watersheds.

Discharge measurements at 1-hr intervals at the Lighvan and Hervi gauges are available for 13 events from 1998 until 2012. Rainfall data were collected at Lighvan station by means of a weighing recording rain gauge.

The rainfall of this rain gauge was considered as the input rainfall for both sub-watersheds. The rainfall was assumed to be distributed uniformly over both sub-watersheds because of the small size of the sub-watersheds. The direct runoff hydrograph which can be separated by straight line, fixed based, and variable slope methods, was specified by the fixed gradient method as a simple and commonly used technique (Chow et al. 1988) for each event. It should be noted that the probable bias in this simple method would be offset in the calibration process. Thereafter, the Φ-index method, which applies the continuity equation, was selected among the other methods such as curve number (CN) or SCS to extract the excess hyetograph of each event from observed hyetographs (Chow et al. 1988). In this way, any rainfall prior to the beginning of direct runoff was taken as initial abstraction. Since the study area is small, the rainfall distribution was assumed uniform over the whole watershed. In this study, 13 storm events occurring simultaneously in both sub-watersheds (during 1998–2012) were selected to examine the proposed methodology. Assuming the linear system theory for instantaneous unit hydrograph (IUH) model, the runoff hydrographs for the two sub-watersheds were separated using the concept of linearity. Ten storm events were used for calibrating the model and three events used for verification. Characteristics of the available events are presented in Table 2.

A well-developed global archive of Landsat images with 30 m spatial resolution is available, and is widely used to detect and monitor land cover change (Kepner et al. 2008). Satellite images for the study area (located on path 168, row 34) for the rainy season months May (3), June (2), and July (3) in the period 1998–2012 were downloaded from the USGS website (https://earthexplorer.usgs.gov/login/). Since land cover typically does not change rapidly, one image was assigned to each rainfall–runoff event based on the closest cloud-free image to the date of event. The weather around a storm event's day is usually cloudy, also in some cases there is no image on the exact day of event. Thus we tried to select the cloud-free image of day which is closest to the event date to calculate the NDVI, representing the land cover of the study area for the certain event. Accordingly, seven images were selected for the period 1998–2012, and an extra image from July 1985 was selected to represent land cover prior to the study period. Characteristics of the selected images are presented in Table 3.

The Landsat images were processed to compute the NDVI for each image pixel. Available cloud-free satellite images were used to detect the land cover changes of the study area.

Since the Landsat images are from different times and seasons, the position of the sun, the angle of the terrain, and the atmospheric effects are different. Also, different sensors might not have the same image data. In order to make comparisons between the images obtained by different sensors in different times and seasons, the Landsat images have to be pre-processed by the common methods. The image pre-processing carried out in this study includes radiometric, atmospheric, and topographic corrections. The radiometric correction decreases signal variations uncorrelated to the brightness of the image surface. The purpose of atmospheric correction is to eliminate the atmospheric effects to specify true surface reflectance values. This correction is done by the FLAASH method in ENVI (ENVI 2009). For a DEM, the topographic correction removes the effects resulting from the differences in illumination due to the position of the sun and the angle of the terrain (Lu et al. 2008). This topographic correction has been already applied by the provider (see, http://data.cgiar-csi.org/srtm/tiles/) to the SRTM DEM used in this study. The Landsat images were all geo-referenced to the same co-ordinate system and there was no need to apply geometric correction as this is done by the U.S. Geological Survey using globally available digital topographic data.

Bands 3 (red visible) and 4 (near-infrared) were corrected and the NDVI was calculated for each image. The watershed boundary was then used to subset the NDVI layer via region of interests (ROIs) operation (ENVI 2009).

The main concept behind the NDVI is that for vegetated surface, red and near-infrared wavelengths are characterized by high and low absorptions, respectively (Chen et al. 2003). When sunlight strikes objects, certain wavelengths of this spectrum are absorbed and other wavelengths are reflected. The pigment in plant leaves, chlorophyll, strongly absorbs visible light (from 0.4 to 0.7 μm) for use in photosynthesis. The cell structure of the leaves, on the other hand, strongly reflects near-infrared light (from 0.7 to 1.1 μm). The more leaves a plant has, the more these wavelengths of light are affected, respectively. In general, if there is much more reflected radiation in near-infrared wavelengths than in visible wavelengths, then the vegetation in that pixel is likely to be dense and may contain some type of forest. If there is very little difference in the intensity of visible and near-infrared wavelengths reflected, then the vegetation is probably sparse and may consist of grassland, tundra, or desert. Healthy vegetation absorbs most of the visible light that hits it, and reflects a large proportion of the near-infrared light. Unhealthy or sparse vegetation reflects more visible light and less near-infrared light. Nearly all satellite vegetation indices employ this difference formula to quantify the density of plant growth on the Earth. The NDVI value for a specific pixel always varies from −1 to +1. Water bodies are specified by extreme negative values, surfaces with no vegetation cover result in near to zero NDVIs, and the highest density of green leaves is indicated by NDVI value close to +1 (0.8–0.9). NDVI indicates the vegetation cover level and also acts as a beneficial index for monitoring vegetation variations and land cover changes.

The aforementioned steps were applied to obtain the NDVI maps of the study area for all images. The ‘mean NDVI’ value (NDVI_w) indicates the average of NDVIs in all pixels of land cover image of watershed. It should be noted that the negative values of the pixels (water bodies) reduce the mean values of NDVI. As these negative values exist in all images, they would affect all images during the whole study period and so the mean values would be comparable.

In this study, the geomorphologic rainfall–runoff model of GUHCR with a single model parameter is investigated as a semi-distributed conceptual model. There are some other conceptual models such as Zoch and Nash models (Singh 1988), which could be used here, but this model was selected as a semi-distributed model to analyze the relation between the model parameter – which relates to the geomorphologic properties of the watershed – and NDVI.

The root mean square error (RMSE) measure and a related normalization called the Nash–Sutcliffe efficiency (NSE, Nash & Sutcliffe 1970), are the two model performance criteria most widely used for calibration and evaluation of the hydrological models by comparison with observed data (Gupta et al. 2009). Legates & McCabe (1999) indicated that a hydrological model can be sufficiently evaluated by NSE (Equation (8)) and RMSE (Equation (9)), but owing to the importance of peak discharges in flood control and the volume of runoff in water resources management, three other important criteria, the ratio of error of peak flow (E_P), the ratio of error of time to peak (t_P), and the ratio of error of the hydrograph's volume (E_v) were also used in this study to evaluate the proposed methodology as shown in the appendix (available with the online version of this paper). There are some other criteria that can be used in order to assess the performance of the method. The used criteria are the most widely used measures for calibration and evaluation of the hydrological models by comparison with observed data.

The model parameter was calibrated based on minimizing (maximizing) the RMSE (NSE) and then the performance evaluation was done using the other criteria.

There are several methods for calibration of model parameters such as PEST (Doherty 2007), and meta-heuristic methods like the genetic algorithm (Nourani 2008). In this study, the calibration of the parameter was performed using a semi-automatic trial and error method because there is only one parameter in the model.

Using the calculated values of NDVI and the calibrated parameter of the GUHCR model for each event, a regression equation could be fitted to obtain a relationship between NDVI and the parameter of the model for the study area. A schematic diagram of the proposed methodology is shown in Figure 2.

At the first step, the hydro-climatologic (rainfall–runoff events and geomorphological characteristics of watersheds) data were gathered for both watersheds and the contemporarily cloud-free land cover images were processed to calculate related NDVI values. Then, hydro-geomorphologic data were applied to calibrate and verify the conceptual rainfall–runoff model (GUHCR). Finally, statistical relationship was created between model parameter and NDVI values.

To develop relationships between the model parameter and temporal trends of NDVI due to the land cover change, the 13 storm events occurring simultaneously over both Lighvan and Hervi sub-watersheds were selected and used along with the corresponding Landsat images. All available 13 events data were used for calibration and verification of the hydrological model to reduce the uncertainty of the calibrated model parameter. However, for detecting the land cover changes (via NDVI), since there are not significant changes in land cover in a short period of each year (within 2, 3 months), only one image (most appropriate image) for each year (totalling seven images, see Table 3) was used to find long-term relationship between NDVI of image and the calibrated parameter of the hydrologic model.

Figure 3 shows the NDVI maps, derived from the respective Landsat images. NDVI values smaller than 0.1 correspond to no vegetation or the urbanized part of the watersheds, while the higher NDVI values correspond to higher vegetation coverage. In order to investigate seasonal NDVI dynamics, a representative dynamic and stable pixel for each sub-watershed was selected.

Figure 4 shows the location of pixels, also changes of NDVI in selected pixels for both sub-watersheds. As can be seen in Figure 4, one of the group pixels in both sub-watersheds (Pixel 1) does not change significantly during the seasons while the other pixels which contain dense vegetation cover clearly show a dynamic NDVI behavior. This indicates the seasonal change of land cover. Constant NDVI value for stable pixels demonstrates there are no/limited inconsistencies in derived NDVI maps. Considering both kinds of these pixels, the values of NDVI may have been determined with some degrees of uncertainty, but overall, a decreasing trend can be detected from the NDVI values (see Figure 3 and Table 3). The mean NDVI values of the start and end dates of the study period (1998 and 2012) were respectively computed as 0.283, 0.210 for the Lighvan sub-watershed, and 0.219, 0.172 for the Hervi sub-watershed.

The average of mean NDVI values for all images in Hervi and Lighvan sub-watersheds are respectively 0.17 and 0.21. As can be seen clearly from the images (Figure 3), the NDVI values generally decreased during the 14-year study period (as expected) for both sub-watersheds. However, because of inter-annual changes inconsistent behavior is observed. The changes in the Hervi sub-watershed are more considerable, mostly due to land cover changes, with the northern part experiencing a considerable loss of vegetation between 1998 and 2012 (see Figure 3). In general, the residential portion of the Lighvan sub-watershed is smaller and so the NDVI values are larger, with changes being mainly due to a land-use shift from grassland to farmland. The following sub-sections present and discuss the relationships between the model parameter and temporal evolution of NDVI.

The geomorphological parameters of the GUHCR model (C_i, S_0i, L_i) were first extracted using GIS tools (Table 4). Then, the parameter was determined via Equations (5)–(7), using the sub-watersheds' physical parameters (Table 4) and data of 10 calibration storm events.

Table 5 displays the parameters of the GUHCR models during 1998–2012 calibrated for each rainfall–runoff event. The GUHCR model parameter is related to the watershed's geomorphology and the event's characteristics. As it can be seen from Table 5, parameter decreases from 1998 to 2012 in both sub-watersheds, although it has more fluctuations in Lighvan which has smaller urbanized areas (see Figure 3).

Table 6 summarizes the model performance criteria (NSE, RMSE, Ep, t_P and Ev) for the calibration events which indicate acceptable performance of the GUHCR model. Based on the E_p, t_P and E_v values in Table 6, the GUHCR model simulates the peak flow, time to peak, and the volume of runoff accurately. The results of Table 6 show that the variation of parameter before 2000 is greater than those in the other years for the study period. This may be associated with more rapid increases in urbanization and gardening, and decrease in grassland and farmland. Beginning around the year 2000, the study area was designated as an experimental watershed by East Azerbaijan local water organization, and human activities were restricted.

The previously shown results for the study sub-watersheds demonstrate the effect of land cover changes on the results of the rainfall–runoff model. As mentioned before, the changes of the model parameter are somewhat similar to NDVI changes in the study period.

As can be seen in Tables 3 and 5, both and NDVI values show a decreasing trend from 1998 to 2012. NDVI change in time mainly reflects the impact of human activities, urbanization, and conversion from sparse to dense pasture, as well as seasonal changes. Parameter of the model for both sub-watersheds can therefore be related to the NDVI values through a regression analysis.

Figure 5 illustrates the relationship between the parameter and the decreasing percent of average NDVI values (I_ndvi) during 1998–2012 for Hervi and Lighvan sub-watersheds, respectively. It should be mentioned that the parameter is dependent on each rainfall–runoff event; also, the impact of climate condition is implicitly involved in both rainfall and runoff behaviors and data.

The NDVI values of year 1985 for each sub-watershed (Figure 3, NDVI values are 0.231 and 0.286 for Hervi and Lighvan sub-watersheds, respectively) was considered as the reference NDVI and the percent of NDVI decreasing (I_ndvi) during 1998–2012 was calculated with regard to this reference value. While a power form of the regression equation could be best fitted to the data on parameter using available images of Hervi sub-watershed, for the Lighvan sub-watershed the fitted equation shows somewhat lower correlation between the parameters due to more fluctuations of the parameter. This might be due to the more limited land cover changes in the Lighvan sub-watershed during 1998–2012. One of the most important features of land cover change can be urbanization. This means that the sub-watershed has been mostly covered by vegetation cover with highly temporal fluctuations (from each season and one year to the next) during the study period.

Table 7 presents the verification results of GUHCR model in the Hervi sub-watershed. As is clear in Table 7, the verification results of the GUHCR model show acceptable accuracy in predicting the peak flow, time to peak, and the volume of the hydrograph.

Figure 6 compares the observed and simulated hydrographs of the verification events using the GUHCR model for both sub-watersheds. As can be seen from Table 7 and Figure 6, although the output hydrographs of Hervi sub-watershed are accurately predicted, hydrographs of Lighvan sub-watershed are not well simulated. This may be due to temporal oscillation of the NDVI values of the Lighvan sub-watershed which leads to a weak correlation between the model's parameter and the NDVI values. Hervi sub-watershed shows a clear trend during the study period; however, Lighvan does not show a certain trend in vegetation cover and thus the parameter.

Based on the calibration and verification results, it can be seen that the values of the criteria are reasonable for Hervi sub-watershed, however not so satisfying for the Lighvan sub-watershed. The average values for NSE, RMSE, E_p, t_P, and E_v at the Hervi sub-watershed using the GUHCR model are 0.87, 0.024, 2.59, −15.00, and −12.65, respectively. In contrast, the averages of the aforementioned criteria for Lighvan sub-watershed are 0.55, 0.031, −34.65, −52.78, and −21.46, respectively.

In the current study, it was attempted to examine a proposed geomorphology-based model and remote sensing data to relate the land cover properties of the watershed to the model parameter. This would lead to a more reliable result due to the precise detection of the land use/cover conditions, considering seasonality and permanent land cover changes. Applying this methodology, the model's parameter is determined using the temporal NDVI values. A clear trend of NDVI can be detected in Hervi sub-watershed in the study period. However, in Lighvan sub-watershed, the NDVI values are fluctuate significantly because of a smaller urbanized area and larger natural land cover, and could be more related to climate parameters (seasonality) than anthropogenic effects.

Figure 6 shows a good agreement between simulated and observed hydrographs in the Hervi sub-watershed. It is also worth noting that the GUHCR model is able to simulate both the rising limbs of hydrographs (which are mostly related to storm properties) and the recession limbs (which are usually related to the watershed morphology) appropriately. As indicated by E_p, t_P, and E_v in Table 7, the error of peak flow, time to peak, and the volume of hydrographs show acceptable accuracy for the Hervi sub-watershed. However, the model performance is acceptable in calibration events in both watersheds; the performance criteria are not good in the Lighvan sub-watershed due to the use of regression equation which is poor in Lighvan.

Note that the GUHCR model has an acceptable performance in both sub-watersheds, which is attributable to the use of geomorphologic properties for calculating the model parameter. This model is dependent on physical properties of the watersheds which have less uncertainty in comparison to the historical rainfall–runoff event data (Nourani et al. 2015). Such geomorphological models, which use the physical properties of the watershed to estimate the model parameters, can be employed in ungauged watersheds lacking sufficient data due to less dependency of the model on rainfall–runoff event data. It should be noted that most rainfall–runoff events over the watershed are occurring only in a few months of the year which reduces the uncertainty of the calibrated parameter.

In the GUHCR model, there is only one parameter which is calculated using the geomorphologic properties of the watershed. The model parameter is related to area, slope, and length of the watershed; however, the other geomorphologic parameters can be taken into account in the model. This calibrated parameter is related to land cover variation via a power equation using remote sensing data and can be estimated using the NDVI values. A single power equation has been fitted to create a relation between the NDVI and model parameter; however, having more data for fitting would lead to a more precise model or even to developing an artificial intelligence-based model.

Based on the result of the verification step, higher performance criteria are clearly due to the higher correlation between the NDVI and model parameter in Hervi sub-watershed. The high fluctuation values of NDVI in Lighvan sub-watershed might be more related to climate parameters than human effects. Also, this might be due to less uncertainty in watersheds with urbanized parts where runoff to rainfall ratio is larger than in the natural watershed. For small runoff to rainfall ratios in the natural sub-watershed, the uncertainties associated with rainfall measurement errors and spatial variability can form a dominant part of the overall model prediction uncertainty (high noise to signal ratio). In addition, the parameters' values (i.e., physical and land cover properties) in the urbanized part of a watershed are approximately constant while these values are more variable over the season in natural watersheds. This characteristic would show less uncertainty in the urbanized watershed which leads to a more precise result at Hervi sub-watershed than that of the Lighvan sub-watershed. Here, the purpose was to find an overall trend for the model parameter; however, because of fluctuations in the model parameter around the mean value, the fitted equation could not show an acceptable correlation.

In former studies, the parameters of the Nash model were related to urbanization using only the factors of the impervious area and population density via a power equation (Cheng & Wang 2002; Huang et al. 2008, 2012; Cheng et al. 2010). In the mentioned studies, the population was considered as one of the urbanization factors. In countries with population control programs the urbanization development is not related to the population and therefore this factor would not be an indicator of urbanization in such countries. Table 8 presents a summary of the relevant literature on the effect of land cover changes in rainfall–runoff modeling showing the present study in the last row. As can be seen in Table 8, the effects of land cover changes in conceptual models have been calibrated only based on the event data. The invention of GIS also considering the geomorphologic properties of the watersheds has led to the use of semi-distributed rainfall–runoff models. Although many studies have studied the effect of land cover properties on rainfall–runoff models, there are few studies that have investigated the temporal land cover properties of a watershed using the NDVI and the parameter of a semi-distributed model. As mentioned before, in the present study we tried to relate the changes of NDVI to the changes of the model parameter to address the gap of previous studies (see Table 8). Using five criteria to evaluate the investigated model, it could be demonstrated that the obtained regression equation provides a promising approach for estimating lumped model parameters using remote sensing data for one of the study sub-watersheds.

In this study, the remotely sensed data were used to analyze the NDVI values during 1998–2012. Considering a reference value for the NDVI, the percent of NDVI decreasing (I_ndvi) during the study period was calculated, which was used to determine a relationship between the model parameter and NDVI. Regression relationships were established between temporal variations in remotely sensed NDVI and the parameter of the geomorphologic IUH model. Thereby, a mechanism is provided to monitor and detect the temporal variation of watershed land cover and its consequent impacts on the rainfall–runoff response of the watershed. A single power equation was fitted on the NDVI and the model parameter values. However, having more data would contribute to a more precise model.

The obtained results indicate that decreasing NDVI percentages (e.g., due to urbanization or other land cover/use changes, such as seasonal variation) may show a considerable impact on values of the model parameter. While a power form of the regression equation could be best fitted to the data on parameter using available images of Hervi sub-watershed, the fitted equation for the Lighvan sub-watershed implies slightly lower correlation between the parameters due to more fluctuations of the parameter. As a result, despite the accurate prediction of the output hydrographs of Hervi sub-watershed, hydrographs of Lighvan sub-watershed are not well simulated. This might be the consequence of temporal oscillation of the NDVI values in the Lighvan sub-watershed which leads to a weak correlation between the model's parameter and the NDVI values. To evaluate the performance of the fitted equations, three rainfall–runoff events were examined to verify the rainfall–runoff modeling for both sub-watersheds. The verification results of the GUHCR model show acceptable accuracy in predicting the peak flow, time to peak, and the volume of the hydrograph. The average NSE values for the GUHCR model in the presented study are 0.87 and 0.55 in Hervi and Lighvan sub-watersheds, respectively.

The NDVI analyses of the ‘stable pixels’ of both sub-watersheds show that combination of seasonal and inter-annual changes of NDVI would lead to an overall downwarding NDVI which, in turn, can affect the rainfall–runoff model parameters.

Our results confirmed that the single parameter of the GUHCR model is related to the geomorphologic properties of the watershed, which has a physical meaning. Employing the physical properties of the watershed to calculate the model parameters, such geomorphological models can be efficiently applied in ungauged watersheds lacking sufficient data due to less dependency of the model on rainfall–runoff event data.

As a future research plan, it can be suggested to employ other conceptual IUHs and geomorphologic IUHs or fully distributed models using the proposed methodology. Also, the effect of temporal and seasonal changes of NDVI on models parameters can be analyzed more carefully considering the seasonal variations in type and phenology of vegetation cover. It can be suggested to obtain a combined description by different linear equations which would have better performance. Implementation of other remote sensing indices with different characteristics and abilities (such as LAI, instead of NDVI) to examine the land cover change effects on the performance of the hydrologic model can be considered as a research plan for future studies. The results of different climatic watersheds can be investigated and compared with each other to obtain a global relationship.

The applied modeling and method have some limitations just like other models. The modeling was performed assuming the watershed as a linear system which is nonlinear in reality. This assumption might result in some errors in the model's output. Furthermore, the number of rainfall–runoff events was limited to some seasons, which might be a limitation in determining the parameters' trend; however, the model can be updated to have access to more data in the future.

This paper was supported by a research grant of the University of Tabriz. We thank East Azerbaijan local water organization in Tabriz, Iran, that provided the events data related to the study area (to obtain the data contact [email protected]). The satellite images of the studied area were downloaded from the USGS website (https://earthexplorer.usgs.gov/login/).