Impervious surface area (ISA) has different surface characteristics from the natural land cover and has great influence on watershed hydrology. To assess the urbanization effects on streamflow regimes, the authors analyzed the U.S. Geological Survey (USGS) streamflow data of 16 small watersheds in the White River [Indiana (IN)] basin. Correlation between hydrologic metrics (flow distribution, daily variation in streamflow, and frequency of high-flow events) and ISA was investigated by employing the nonparametric Mann–Kendall method. Results derived from the 16 watersheds show that urban intensity has a significant effect on all three hydrologic metrics. The Variable Infiltration Capacity (VIC) model was modified to represent ISA in urbanized basins using a bulk parameterization approach. The model was then applied to the White River basin to investigate the potential ability to simulate the water and energy cycle response to urbanization. Correlation analysis for individual VIC grid cells indicates that the VIC urban model was able to reproduce the slope magnitude and mean value of the USGS streamflow metrics. The urban model also reproduced the urban heat island (UHI) seen in the Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature products, especially for the grids encompassing the city of Indianapolis, IN. The difference of the hydrologic metrics obtained from the VIC model with and without urban representation indicates that the streamflow regime in the White River has been modified because of urban development. The observed data, together with model analysis, suggested that 3%–5% ISA in a watershed is the detectable threshold, beyond which urbanization effects start to have a statistically significant influence on streamflow regime.
Urbanization, one of the most pervasive land conversions by human activities (Alberti 1999), has a great effect on land cover and landscape characteristics, resulting in a sharp increase in impervious surface area (ISA) on the global landscape. Under natural conditions, precipitation is intercepted by vegetation or stored in the soil body, from which it is taken up by vegetation to evapotranspire or drain slowly into streams (Dunne and Black 1970). In highly urban landscapes, more water flows over impervious surfaces or through drainage networks that can greatly increase the speed of flow to receiving water bodies (Arnold and Gibbons 1996). Increased ISA because of urbanization also leads to greater temperatures in urban areas than in nonurban areas, resulting in urban heat islands (UHIs; Oke 1982). Urban land cover–land use change (LCLUC) can demonstrably reduce infiltration rates on the land surface and runoff response times (Gregory et al. 2006), resulting in higher flow peaks and larger total streamflow volume, shifts in subsurface flow to surface flow, and increases in flood frequency (Konrad and Booth 2005; Tang et al. 2005).
The effects of urban expansion on streamflow patterns have been investigated in various studies as urbanization becomes a global phenomenon. Among all the hydrologic responses, floods in urban environments have received the most attention. This is because with increasing ISA, increasing flood frequency in urban areas poses a greater threat to population and infrastructure than corresponding flood magnitudes in rural areas (Changnon and Demissie 1996; Moscrip and Montgomery 1997; Rose and Peters 2001). Using an upstream versus downstream comparison, Sheeder et al. (2002) found two distinct peaks in outlet discharge, each representing different contributing areas covered by different urban fractions. With the recent realization of the important role of streamflow in fluvial ecosystems and ecological diversity, many researchers have also examined urban effects on biological systems. Poff and Allan (1995) and Clausen and Biggs (1997) identified that hydrologic variability is a significant environmental variable influencing lotic ecosystems. Thus, increasing impervious cover in a watershed has a substantial influence on biological degradation (Booth and Jackson 1997; Wang et al. 2001; Booth et al. 2004).
Much of the previous research on hydrologic response to urbanization relies on empirical studies of daily streamflow observations, such as those obtained from the U.S. Geological Survey (USGS). The problem is that for many streamflow gauge stations, continuous records may not be available for sufficiently long periods (e.g., 50 years) to isolate urban effects, or gauges are not placed to specifically monitor the effect of human regulations. When trend analysis is performed on short-term streamflow data, the results may be more sensitive to climate variability and water resource management decisions than to land use changes. Another challenge is that it is difficult to assess causes of change when analyzing streamflow alone, because hydrologic response is affected not only by LCLUC but also by climate variability (Pielke et al. 2007). Because model simulations can be continuous over long periods, and modeled scenarios can serve as a true control, trends due to land use change will tend to be more statistically significant using simulated streamflow from a physically based hydrologic model (Bowling et al. 2000).
Overall, urban hydrology may have an influence on ecosystem sustainability that has been overlooked in the first-generation studies. At the same time, UHI has a direct effect on regional ecosystems. Predictive tools that can examine both of these major effects are largely lacking, despite the fact that they are linked through the surface water and energy cycles. In this study, a bulk parameterization approach was employed to modify the Variable Infiltration Capacity (VIC) macroscale hydrologic model to make it more suitable for urbanized watersheds, resulting in the VIC urban model, as discussed in section 3. We hypothesize that the streamflow regime is modified by increasing urban intensity, and the VIC urban model can capture these changes in the White River basin as well as the UHI pattern. To test these hypotheses, in section 4 USGS daily streamflow data from 16 small watersheds with different degrees of urbanization in the White River, Indiana (IN), were analyzed to assess differences in hydrologic metrics for watersheds with different amounts of urban land cover. The VIC urban model was then applied in the White River basin to evaluate the simulated hydrologic metrics for 100 individual grid cells with different impervious cover. To facilitate these analyses, the default unit hydrograph used in the VIC streamflow routing model was first improved to better represent the mean characteristics of the White River basin. Simulated and observed hydrologic metrics and simulated and observed surface temperature from the Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature (LST) products were compared in section 5 for the White River basin. Following model evaluation, the VIC urban model was used to assess the critical fraction of urbanization required in the White River watershed to initiate a detectable hydrologic response.
a. Study area
The study area is the 12 000 km2 White River basin above Newberry, IN (Fig. 1). The entire basin is located on a Wisconsinan Age glacial till landscape, representative of large parts of the agricultural midwestern United States. The nature of the till includes low topographic gradients and underlying dense till that restricts vertical water movement, resulting in a functionally homogeneous soil distribution. According to the Koppen climate classification scheme, this basin has a humid continental climate, characterized by variable weather systems and large seasonal variations in temperature. The annual temperature across the study domain exhibits a northeast–southwest increasing gradient, ranging from 11.4° to 13.7°C on average, with some urban regions having a higher relative temperature. The annual precipitation across the domain exhibits the same increasing gradient as temperature, ranging from 1000 to 1200 mm. Precipitation is generally delivered by frontal systems in the winter, and by convective thunderstorms in the summer, resulting in a fairly uniform seasonal distribution of precipitation magnitude but a more seasonal cycle in precipitation intensity.
The White River basin includes the city of Indianapolis, which has experienced a high rate of urbanization during the past 30 years. To investigate how urban expansion affects local hydrologic response, two land use maps with spatial resolutions of 100 m for the years 1983 and 2001 were employed (Figs. 2a,b). The land use–cover classes were derived from an unsupervised classification of Landsat thematic mapper (TM) images and placed in a modified Anderson level-2 classification scheme (Wilson and Lindsey 2005). For this study, they were further aggregated into the following classes: high–low density urban (7.4%, based on 2001 land cover map), road (0.8%), forest (19.5%), grassland (21.6%), cropland (49.4%), and water–wetland (1.2%). Figure 2c shows that most of the increase in urban pixels is around the city of Indianapolis.
b. Data sources
Sixteen small watersheds with different degrees of urbanization were selected within the White River basin for hydrologic analysis based on the availability of USGS daily streamflow data between 1 October 1979 and 30 September 2003 (see Fig. 3). The drainage areas range from 20 to 300 km2. The observed daily streamflow data are obtained from USGS Surface-Water Daily Data for the Nation (available online at http://waterdata.usgs.gov/usa/nwis/sw).
Daily USGS streamflow for the White River at Newberry was also used for model evaluation. Ruddy and Hitt (1990) identified eight reservoirs within the White River basin with a normal capacity of at least 5000 acre–feet (6.17 × 106 m3). The total storage volume of the reservoirs, 0.162 × 109 m3, is small (around 3%) compared to the mean annual volume of the White River at Newberry (4.92 × 109 m3). In addition, Arvin and Spaeth (2007) reported that the urban water usage (e.g., the annual water withdrawal for industry and public supply) from 1986 to 2006 in the White River basin was around 70 billion gallons (0.265 × 109 m3), and the amount of water withdrawal did not change substantially during these 20 years. An estimated 27.2% of water withdrawals in the United States go toward consumptive use (0.072 × 109 m3); the rest will be return flows, in this case to the White River (Ward and Trimble 2003). The fraction of the consumptive use to the mean annual volume at the basin outlet is also small (1.5%). It is therefore assumed that reservoirs, and industrial and municipal water withdrawals do not have a significant effect on streamflow measured at Newberry, and the observed streamflow is used without adjustment in this analysis.
The VIC model is forced with observed surface meteorological data, which includes daily precipitation, and daily maximum and minimum temperature. The meteorological data were gridded at the resolution of the VIC model, a scale of ⅛° latitude–longitude, for the entire midwestern United States using the methodology of Hamlet and Lettenmaier (2005). Daily average surface wind speeds are taken from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis fields (Kalnay et al. 1996). Other meteorological and radiation variables are calculated based on the range in daily temperature as documented in Maurer et al. (2002).
The MODIS surface temperature images used in this study are daytime and nighttime 8-day composite images (MOD11A2.4; Wan et al. 2004). The satellite scans the White River region around 1100 and 2300 LT each day. Both daytime and nighttime 8-day-average products were selected to include a typical summertime period, considered as 28 July in this study for the years 2000–03. They were downloaded from the “Land Processes Distributed Active Archive Center (DAAC) Data Pool” of the MODIS Web site (available online at http://lpdaac.usgs.gov/main.asp), and the spatial resolution is 1 km.
c. Impervious area estimation
There are two urban categories in the Wilson and Lindsey (2005) land classification: high-density and low-density urban. High-density urban areas are defined as locations that have been highly developed as commercial and industrial land uses, whereas low-density urban means areas that have been less intensively developed, such as those dominated by single-family residential development (available online at http://www.epa.gov/mrlc/definitions.html). The two urban categories can be assumed to contain different percentages of ISA and proportional permeable land cover. For this project, turf grass, considered as a vegetation category in urban area for recreational and conservational purposes, is used to represent all pervious urban area. According to the Environmental Protection Agency (EPA)’s definition, 20%–49% of low-density urban development is impervious, and 80%–100% of high-density urban development is impervious (available online at http://www.epa.gov/mrlc/definitions.html). In this study, 90% of high-density urban areas and 35% of low-density urban areas were assumed to be impervious—that is, the mean value of EPA’s standard, and the remainder was represented as turf grass. As shown in Table 1, ISA fraction for each small watershed ranges from 0 to 0.5; more than half of the watersheds have ISA fraction less than 0.05. Some watersheds farther from the city center have almost no increase in ISA between 1983 and 2001. The ISA fraction of the entire White River watershed above Newberry is 0.032 for the 1983 land use map, and 0.041 for the 2001 land use map. Most of the ISA is located around the city of Indianapolis.
d. Statistical analysis
1) Mann–Kendall correlation coefficient
Because the distribution of ISA fraction exhibits a high degree of skewness, the MK method (Kendall 1970; Mann 1945) was employed to calculate the correlation between hydrologic metrics (flow distribution, daily variation in streamflow, and frequency of high flow events) and ISA fraction. It is a rank-based nonparametric procedure; therefore, it is resistant to the effect of extreme values and to deviations from a linear relationship. The statistic S of the MK test is given by
where Xi are the sequential data indexed according to the magnitude of i and n is the dataset length. The sign sgn(θ) is given by
When n ≥ 8, the statistic S is approximately normally distributed with a mean of zero (Kendall 1970). The null hypothesis H0 is that the distribution of hydrologic metrics does not change as a function of ISA fraction.
The trend magnitude β is given by (Hirsch et al. 1982)
2) The Wilcoxon signed-rank test
The nonparametric signed-rank test was used to examine whether the hydrologic metrics calculated from the VIC urban and unmodified VIC model were significantly different. The null hypothesis is that the two samples come from the same population, and the alternative is that they are different in level (means or medians; Wilcoxon 1945).
3) Nash–Sutcliffe coefficient of efficiency
During calibration of the VIC model, the efficiency coefficient E was calculated based on Nash and Sutcliffe (1970) to evaluate model performance, in addition to visually comparing the observed hydrograph and modeled hydrograph; E is given by
For daily streamflow calibration, Qobs is the observed daily streamflow, Qsim is the simulated daily streamflow of the VIC model, and Qmean is the mean observed daily streamflow during the calibration period. Simulation results are considered to be a good fit to observed data for E > 0.75, whereas for E between 0.36 and 0.75, the simulation results are considered to have a satisfactory fit to observed data (Motovilov et al. 1999).
e. Hydrologic metrics
To quantify the influence of urbanization on the streamflow regime, the correlation between averaged hydrologic metrics and ISA fraction was identified for the 16 watersheds. The three hydrologic metrics—flow distribution, daily variation in streamflow and frequency of high-flow events—were calculated annually based on daily streamflow data between the 1980 and 2003 water years (24 years), where the water year starts on 1 October. Although there are many other hydrologic metrics that can be used to measure the streamflow characteristics (Olden and Poff 2003), it has been demonstrated that aquatic ecosystem communities are sensitive to the hydrologic metrics we used (Konrad and Booth 2005). Additionally, these metrics adequately explain the three basic characteristics of hydrologic regime (flow duration, flow variability, and high flow). ISA fractions calculated from the 1983 and 2001 land use maps were compared to the averaged annual metrics from 1980 to 1991 and from 1992 to 2003, respectively. One of the important considerations in urban planning is the relative change of hydrologic metrics with increasing ISA fraction in an urbanized watershed. The analysis of hydrologic metrics therefore focused more on the change magnitude (the slope of the correlation) than on the absolute value of metrics.
1) Flow distribution
The fraction of time that daily streamflow exceeds mean streamflow for each year (Tqmean) was calculated to test flow distribution. Tqmean is based on the duration rather than the volume of streamflow and is given by Konrad and Booth (2002):
where mean streamflow is recalculated each year. In response to increasing ISA, surface runoff (hence storm event flow) is expected to increase, generally at the expense of baseflow. There is often no clear trend in mean annual flow in the absence of significant water withdrawals. The redistribution of water from baseflow to storm flow in urban regions, as well as the increased flashiness of the system, can cause a shorter period of the year when streamflow is greater than mean streamflow, thus decreasing the Tqmean (Konrad and Booth 2005).
2) Daily variation of streamflow
The Richards–Baker flashiness index (R–B index) is used to evaluate the daily streamflow variation. The index is calculated by dividing the sum of the absolute values of day-to-day changes in daily discharge volumes (pathlength) by total discharge volumes for each year. The index, as modified from Baker et al. (2004), is as follows:
where Qi is daily discharge volume for day i and n is the number of days in a year. Daily variation is expected to go up with increasing ISA in urban regions (Konrad and Booth 2005).
3) Frequency of high-flow events
A high-flow event is assumed to be a flow event that exceeds 3 times the median flow. A high-flow event may continue for several days, until the daily value drops below the threshold. The median flow used for this study was calculated based on daily streamflow from 1980 to 2003; therefore, the median flow threshold was constant for 24 years. The frequency of high-flow events in a year is expected to increase with urban development (Konrad and Booth 2005).
3. The Variable Infiltration Capacity model
a. Model description
Among many land surface hydrology models, the VIC model is distinguished by its variable infiltration curve and a nonlinear relationship between baseflow and deep soil moisture (Liang et al. 1994). It is also able to represent multiple vegetation types and soil layers within each grid cell. In recent years, it has been applied in many diverse environments at multiple spatial scales, from ⅛° latitude by longitude to the coarser 2° resolutions used for global applications (e.g., Abdulla et al. 1996; Bowling et al. 2000; Cherkauer and Lettenmaier 1999; Cherkauer et al. 2003; Lohmann et al. 1998a; Nijssen et al. 1997; Nijssen et al. 2003). Abdulla et al. (1996) concluded that the VIC model performs best for humid and subhumid catchments, such as the White River basin. The VIC model was therefore employed to analyze the hydrologic response to varying proportions of urban land cover in the White River basin. The cell size for each grid is ⅛° latitude–longitude resolution. A 1-h time step in full energy mode was used to run the model for streamflow analysis and to obtain the simulated surface temperature.
To evaluate the simulation ability of the VIC model at the basin scale, streamflow in all individual grid cells was routed to the outlet of the White River at Newberry, IN, using the Lohmann et al. (1996) routing model. This model has two components: 1) it routes the raw runoff and baseflow to the outlet of each grid cell according to its unit hydrograph (UH), and 2) it uses the Saint Venant equations to do channel routing (Lohmann et al. 1996, 1998a).
b. Model modification
In general, there are currently two approaches to represent urban areas in land surface models (LSMs). The first approach is to adjust the soil constants (e.g., volumetric heat capacity and thermal conductivity) and parameters (e.g., surface albedo, roughness length, and moisture availability) to make them suitable for urban impervious areas. This method can be considered a bulk parameterization approach (i.e., slab model), which means it treats the urban geometry as a flat surface. The second approach is to couple an urban canopy layer model with an existing LSM (e.g., Taha 1999; Niyogi et al. 2006). Although it is likely to produce better simulation results in small-scale regional studies, this second approach requires detailed information, such as urban geometry, canyon orientation, and differences in physical constants between buildings and roads. It is a potentially overwhelming task to collect all these data for large regions (Lei et al. 2008). If simulations are applied at a mesoscale–macroscale, the detailed information may be lost when averaging data back to coarse model grids (Becker and Braun 1999).
In previous applications of the VIC model, urban land cover has either been classified as bare soil, or the other land use types have been rescaled to eliminate the urban area; therefore, urbanized watersheds were not well represented. In this study, the VIC model was improved by applying a bulk urban parameterization approach. First, the ISA was separated from bare soil; the fraction of each part was determined from the land use map. There is no infiltration on the ISA part, and the precipitation falling on the ISA becomes runoff directly. The “turf grass” category was also added to the vegetation library, with albedo = 0.2, height = 0.1 m, and surface resistance = 50 s m−1. In comparison, Allen et al. (1998) recommended the following values for simulation of a well-watered grass reference crop: albedo = 0.23, height = 0.12 m, and surface resistance = 70 s m−1, which have been used in the FAO-56 Penman–Monteith equation. It was assumed in our simulations that 50% of the turf roots are in the top 10 cm, with a maximum root depth of 30 cm. Monthly leaf area index (LAI) values for turf grass are based on the grassland category in the Myneni et al. (1997) global LAI database, resampled as described by Maurer et al. (2002).
The thermal properties of the ISA were also modified to approximate urban values, rather than using those of bare soil. The bulk parameterization is implemented according to Taha (1999) as follows: 1) the thermal conductivity of the ISA is set larger than 3 W m−1 K−1; 2) the volumetric heat capacity (VHC) of the ISA is set larger than 3 MJ m−3 K−1; and 3) the albedo of the ISA is set to 0.15.
A limited sensitivity analysis of the effect of these parameter values on simulated surface temperature indicates that increasing both thermal conductivity and volumetric heat capacity of impervious surfaces can slightly decrease daytime temperature and increase nighttime temperature; decreasing the albedo can increase daytime temperature without changing nighttime temperature. Within a reasonable variation for these thermal properties (i.e., 10%), the simulated temperatures have less than a 2% change for grids with ISA fraction larger than 0.3, and the changes are much smaller for grids with smaller ISA fraction. The bulk parameterization was therefore implemented using the previously mentioned values. In this way, the VIC model was made more suitable for urban watersheds.
c. Model parameters
The VIC model requires characteristics for the dominant soil type of each grid cell, as well as fractional vegetation coverage and vegetation characteristics. The vegetation fractional coverage was calculated from the Wilson and Lindsey (2005) land use maps. Vegetation characteristics such as architectural resistance, monthly LAI, and roughness length for each vegetation type are taken from Mao and Cherkauer (2009). Most of the soil parameters for each grid were extracted from the continental United States soil characteristic dataset, CONUS-SOIL database (Miller and White 1998), where the State Soil Geographic (STATSGO) soil classifications have been adjusted to account for known problems across state boundaries; these data are distributed as gridded for the conterminous United States.
Preliminary analysis indicated that the default UH used in the VIC streamflow routing model does not represent the timing of water movement in headwater streams within the White River basin appropriately. On the basis of observed daily precipitation and streamflow data, more than 20 storm events in the summer and fall were chosen from small watersheds (less than the size of one VIC grid cell) in the basin with low impervious area. For each storm event, the baseflow was separated using the Web-based hydrograph analysis tool (Lim et al. 2005). The daily direct runoff fractions, or UH ordinates, for a storm event were determined by subtracting estimated baseflow from the observed daily streamflow and dividing by the event total direct runoff. These steps were repeated for all of the selected storm events. An aggregated UH for use in the VIC routing model was derived based on examination of all of the events, with values of 0.65, 0.17, 0.11, 0.05, and 0.02 for lags (days) 0–4, respectively, and 0 for all lags larger than 4. In reality, the UH is likely to vary with drainage area, basin shape, basin slope, and urban land cover. At present, the routing model does not allow representation of spatially varying UH functions; however, we believe that the revised UH better represents the mean characteristics of the White River.
Using results from two different land surface scheme intercomparison projects, Lohmann et al. (1998b) and Nijssen et al. (2003) indicated that calibration of hydrologic models can generally improve the streamflow simulation, by refining parameters that cannot adequately be determined a priori. According to Liang et al. (1994), the six parameters that have the largest effect on the hydrograph, the infiltration parameter—bi, the three baseflow parameters Ds, Dsmax, and Ws, and the second and third soil layer thicknesses (Table 2)—are not directly observed in most cases and are difficult to infer from other parameters. These parameters were therefore adjusted during calibration with a manual iterative approach by visually comparing the routed hydrograph and observed streamflow hydrograph, and by maximizing the daily Nash–Sutcliffe (N–S) index.
a. Model evaluation
The final values of the calibrated parameters are provided in Table 2. The daily Nash–Sutcliffe index for the calibration period (1980–91) is 0.76, and 0.86 for the monthly flows. During the model validation period (1992–03), the daily Nash–Sutcliffe index is 0.77, and 0.87 for the monthly data. The 1983 and 2001 land use maps were used for the calibration period and validation period, respectively. Figure 4 shows the monthly discharge comparison for the calibration period, indicating a general overprediction during the high-flow season. To compare the daily peak flow simulation over the entire period, the peaks-over-threshold (POT) series of simulated versus observed peaks was selected with five matching peaks on average per year. The percentage bias of the POT series was −1.0%, which indicates there is no clear bias for the daily simulated peaks.
To assess model performance separately for different climate conditions (such as for dry and wet years), 78 years (water years 1929–2006) of observed precipitation and temperature data, averaged over the White River basin, were used to partition each year into one of four groups. Each year was assigned to wet–warm, wet–cool, dry–warm, and dry–cool, if the annual value was above–below the 78-yr mean for precipitation and air temperature. The daily N–S index for each year was calculated using the observed and simulated (based on the 2001 land use map) daily White River streamflow data and separated into the previously mentioned four groups. The averaged N–S index for each group was then calculated, with values 0.77, 0.72, 0.60, and 0.68 for those four climate conditions, respectively, which means that the VIC urban model performs best under wet conditions. In addition, the same analysis was done for the unmodified VIC model simulation, resulting in averaged N–S indices of 0.77, 0.69, 0.66, and 0.67 for the four climate conditions, respectively. The results suggested that the VIC urban algorithm improves model performance in wet and/or cool conditions, but model performance decreases for dry and warm conditions. Because poor performance for low-flow periods when streamflow is sustained by baseflow is common for many surface water hydrology models, it is difficult to attribute this to the urban component specifically. The results suggest that the VIC urban model is able to capture the high-runoff response associated with precipitation on impervious and saturated ground, but may underpredict baseflow derived from anthropogenic influences in this mixed use watershed.
b. Hydrologic analysis for small “watersheds”
USGS daily streamflow data from 16 small watersheds in the White River basin (IN) and streamflow for 100 ⅛° grid cells simulated using the VIC urban model were analyzed, examining whether the hydrologic metrics Tqmean, R–B index, and high-flow frequency are sensitive to different degrees of urban intensity within these watersheds. When investigating the correlation between hydrologic metrics and ISA fraction, a large degree of scatter is apparent in the distribution when ISA fraction is small (i.e., the metrics do not seem to exhibit stationarity). For example, Fig. 5 shows that observed Tqmean has no clear correlation with ISA fraction when the latter is less than 0.02. (The correlation for watersheds with ISA fraction between 0.02 and 0.04 is unknown because there are few observed data in this range.) Similar results were derived from the simulated hydrologic metrics; for example, there is no clear correlation when ISA fraction is less than 0.05, implying that a small ISA fraction does not have a significant effect on hydrologic response, and other unknown factors dominate the variability in hydrologic regime. In other words, 0.05 ISA fraction would be a potential threshold for detecting the urbanization effects on hydroclimatic datasets. The correlation between hydrologic metrics and ISA fraction was therefore investigated and compared for watersheds and grid cells with more than 0.05 ISA fraction.
Watersheds with ISA fraction larger than 0.05 were extracted from both 1983 and 2001 land use maps. Four watersheds in period 1 (based on 1983 land use) and six watersheds in period 2 (2001) had ISA fractions larger than 0.05, resulting in a total of 10 comparisons. An examination of the precipitation records indicates that there are more extreme storm events in the later period (1992–2003) than in the earlier time period (1980–91), which could lead to an increase in high-flow frequency and R–B index for the later period. The sensitivity of the averaged three hydrologic metrics during the two multiyear periods to individual extremes appears to be small relative to the influence of increasing impervious fraction.
For the best comparison with the small watershed results, the daily “streamflow” for each of the VIC grid cells was obtained by routing the simulated runoff and baseflow, which is distributed nonuniformly over the area within the grid cell (Nijssen et al. 1997), according to the new UH (without the channel routing). The three hydrologic metrics were calculated for each of the individual grid cells with ISA fractions exceeding 0.05. The correlation analyses for both the observed watersheds and simulated grid cells are summarized in Table 3.
Referring to Figs. 6 –8, when ISA fraction increases from 0.05, simulated Tqmean values exhibit a decreasing trend with a magnitude of −0.04 day day−1 (ISA fraction)−1, which means Tqmean decreases 1.5 days yr−1 with every 0.1 increase in ISA fraction. This slope is less than that of the observed data (4.4 days yr−1), while the mean value is overestimated. The simulated R–B index values exhibit an increasing trend with a magnitude of 0.66, indicating flow variation increased by 7% of mean annual flow (MAF) for every 0.1 increase in ISA fraction, lower than the observed value of 15%, and the mean value is also underestimated. The simulated frequency of high-peak flow events exhibits an increasing trend with a magnitude of 38.7 events (ISA fraction)−1, which means frequency increases by 3.9 events for every 0.1 increase in ISA fraction, close to the observed change rate of 5.2 events yr−1. The mean value also matches with the observed one quite well, except the two observed points with large ISA fraction, which exceeds the range of VIC grid cells (Fig. 8). The Mann–Kendall test shows that all of the simulated and observed correlations between hydrologic metrics and ISA fraction are statistically significant (Table 3).
c. Urban heat island
To further evaluate the ability of the VIC urban model to simulate urban processes such as UHI, MODIS LST image products were compared to the VIC model results with and without the urban representation. To quantify the difference in surface temperature due to urbanization, cropland was used to substitute for the ISA in the vegetation parameter file for the unmodified VIC model. Hourly simulated surface temperatures were averaged from 24 to 31 July for the years 2000–03 between 1100 and 1200 LT for daytime and between 2300 and 0000 LT for nighttime temperature calculation. This adjusts the simulated surface temperatures to better match with the MODIS 8-day composite images for both daytime (1130) and nighttime (2330) image acquisition times over the White River region.
The MODIS images indicated unusually low temperatures for a few pixels as a result of cloud contamination; therefore, any pixel values that were less than 20°C in the daytime and less than 10°C in the nighttime were set to a no-data value during the data preprocessing step. The spatial daytime and nighttime temperature distributions that resulted from the MODIS images, the VIC urban model, and the unmodified VIC model are shown in Figs. 9 and 10. The VIC urban model captured the daytime and nighttime UHI patterns observed in the MODIS images for the grids encompassing the city of Indianapolis, whereas the unmodified VIC model did not. To quantify the simulation difference, grids at the same latitude—for example, 39.8125°N (crossing the center of Indianapolis)—were analyzed for the temperature distribution. Each year contained one MODIS image and two VIC model simulations, and there are four years, for a total of 12 comparisons (Fig. 11). Because of the difference in resolution between the VIC model (∼12 km) and MODIS image (∼1 km), the sample size is different; therefore, the box width is dependent on the dataset size. Figure 12 shows a box plot for the nighttime temperature comparison. The model simulation produced less temperature variation, especially for the unmodified VIC model, whereas simulated higher temperatures for grid cells with larger ISA fractions from the VIC urban model were closer to the MODIS image than that from the unmodified VIC model. On the basis of this quantitative analysis, the VIC urban model yielded a large improvement in simulating the UHI around urbanized regions.
d. Evaluation of urban influence on White River streamflow
We have demonstrated that the VIC urban model can reproduce changes in the hydrologic regime due to urban intensity at the small watershed scale and improves the model performance at the scale of the White River. To evaluate the influence of urbanization for the entire White River basin, averaged hydrologic metrics were calculated for 24 years for both the routed and observed daily streamflow for the whole basin. To quantify the difference in streamflow due to urbanization, cropland was used instead of ISA in the vegetation parameter file for the unmodified VIC model (“undeveloped simulations”), as was done for the UHI simulations.
Compared to the undeveloped simulations, the urbanized simulations produced larger runoff during the summer and fall with little change in the baseflow volume, whereas higher evapotranspiration was produced in the undeveloped scenario. In general, restricting infiltration in urban coverage resulted in increased runoff and decreased evapotranspiration, especially in summer and fall where high temperatures and large storm events are more prevalent. The Tqmean of the urban scenario is close to the observed value (Table 4), whereas it is underestimated for the undeveloped scenario. The R–B index of both simulated scenarios overestimates the observed metrics, with decreased performance for the urban scenario. The high-flow frequency was also improved for the urban scenario, which overall indicates the urban scenario better simulated the hydrologic regime at basin scale with a dominant urban area.
Table 4 shows that all three hydrologic metrics were a little larger in the urbanized scenario than those simulated for the undeveloped case. The high-flow frequency is almost 19% higher in the White River as a result of urbanization, whereas the R–B index is approximately 12% greater. Such increases were expected of the R-B index and high-flow frequency, but it was anticipated that Tqmean would decrease with the representation of impervious cover. By comparing the daily streamflow hydrograph of the two models, it was evident that some daily flow events increased to become larger-than-mean flow with the ISA representation. In addition, for some events there were more days with discharge greater than mean flow during the hydrograph recession period for the urban scenario. Both of these situations lead to an increase in Tqmean of about 2.5% for the urban scenario.
To identify if a threshold exists over which the ISA fraction began to affect the streamflow regime at the basin scale, the hydrologic metrics for the White River basin and seven of its tributaries with different ISA fractions and various drainage areas were investigated using both the urbanized and undeveloped scenarios. The nonparametric signed rank test was employed to assess whether each of the two groups of hydrologic metrics were statistically different from each other. As shown in Table 5, each of the two groups of hydrologic metrics for the White River basin and tributaries 1, 3, 4, and 5 are statistically different at the 0.05 significance level (|Z| value > 1.96), except one value in tributary 5. For tributaries 2, 6, and 7 with small ISA fractions (less than 0.04), the two model scenarios are not statistically different for most of the values, with the exception of two R–B index values for tributaries 2 and 6. Table 5 implies that for the bulk representation of ISA in a basin, ISA fractions larger than 0.03–0.05 will have a statistically significant effect on the streamflow regime. The R–B index appears to be the most sensitive, which means streamflow variation could be first affected by increasing ISA.
By comparing the Mann–Kendall slope and mean values of hydrologic metrics for small watersheds, it was found that overall the VIC urban model is able to capture the general trends of the observed data. The underestimation of the three slopes indicates that the simulated results have less variation with increasing ISA fraction than the observed metrics. The reason might be that some uncertainty factors (e.g., watershed size, human regulations), that may counteract or enhance the effect of ISA on hydrologic regime, are not present in the model simulation.
For both simulated and observed metrics, the R–B index and high-flow frequency of the small watersheds were much larger than those at the river basin scale (Tables 3 and 4). This is largely because at the basin scale, the fast runoff response observed for small areas is offset by neighboring regions with small ISA fractions when routing the streamflow along the drainage network. Tqmean does not vary as much as the R–B index and high-flow frequency between the large and small scales, suggesting that Tqmean is less dependent on watershed scale.
The total ISA fraction for the White River at Newberry is about 0.03–0.04 during the 24-yr simulation period, with most of the ISA located around the city of Indianapolis, leaving the surrounding areas dominated by rural land use–land cover types. The analysis based on seven tributaries (Table 5) implies that at the large watershed scale, a threshold of 0.03–0.05 ISA fraction will start to affect streamflow patterns. This threshold is similar to that found in the observational data for our small watersheds and the simulations of individual VIC grid cells (approximately 5% impervious area). On the basis of the previously mentioned results, 3%–5% ISA is suggested as a threshold over which hydrologic metrics can experience statistically significant changes as a result of impervious cover, for both small and large watersheds. Yeo and Guldmann (2006) found that simulated peak discharge increases significantly with increasing urban land cover when the urban area is greater than 8%, whereas Booth and Jackson (1994) found that peak flows are affected when the impervious area reaches 4.6%. Many studies (e.g., Arnold and Gibbons 1996; Booth and Jackson 1997; Booth et al. 2004; Karr and Chu 2000) have found that biological diversity in a watershed starts to decline when the ISA exceeds 10%. Brabec et al. (2002) concluded that macroinvertebrate health declines above a range of 8%–15% impervious area, and physical measurements such as channel enlargement and sediment deposition are affected with 10%–50% impervious area. This implies that the streamflow regime may be affected first with increasing ISA; subsequently, the stream ecosystem is significantly affected when ISA increases to a certain amount, such as 10%. The changes imposed on the natural stream system by increasing ISA—such as streamflow regime, peak flow discharge, stream ecosystem, and then stream channel appearance—therefore can be considered as a continuum rather than a simple and strict threshold (Booth and Jackson 1997).
In this study, by working with 16 small watersheds located in the White River basin in Indiana, it was found that with a 10% increase in ISA in a watershed, the observed Tqmean decreased 4.4 days yr−1, flow variation increased by 15%, and high-flow frequency increased by 5.2 events yr−1. The VIC model was modified to better represent the ISA in urbanized watersheds, with the development of a slab-type urban land use. Similar correlations—1.5 days yr−1, 7%, and 3.9 events yr−1 for the three metrics, respectively—were observed for the three hydrologic metrics computed from simulation results of the White River basin, using 100 individual grid cells with different ISA fractions.
Results based on the simulated streamflow using both the current urban cover and an undeveloped, agricultural scenario suggest that the hydrologic regime in the White River has been modified because of urban development. Overall, urban development in the White River basin has resulted in a 19% increase in simulated high-flow frequency and a 12% increase in flow variability, as measured by the R–B index, relative to an undeveloped, agricultural scenario. The simulated Tqmean also increased by 2.5%.
Results suggested that 3%–5% ISA in a watershed is the threshold above which ISA starts to have a detectable influence on the local streamflow regime. For watersheds with more than 5% ISA, impervious cover was clearly a dominant factor in the trend of the selected hydrologic metrics—that is, negative correlation for Tqmean, positive correlations for the R–B index, and high-flow frequency with increasing ISA fraction. The influence of ISA fractions of less than 3% on streamflow metrics can hardly be identified, and factors such as climate and other land use types mask the effect of ISA.
For surface temperature simulations, both daytime and nighttime temperature in the urban region increased greatly because of the reduced albedo, increased volumetric heat capacity, and thermal conductivity of the urban land use type defined in the VIC urban model. Modified soil thermal properties hold more heat during daytime and then release the excess heat during the nighttime. By comparing with MODIS LST image products, the ability of the VIC urban model to simulate UHI has been demonstrated.
In general, the entire White River watershed does not have much urban coverage; although the city of Indianapolis has experienced a high degree of urbanization, it has been relatively concentrated around the city center. In the future, evaluation of the bulk parameterization approach for the VIC model will include other midwestern U.S. basins that contain very high urban coverage, such as the Milwaukee River (WI) basin. Also, evaluation will move from watersheds under humid conditions to ones with severe heat island cities, such as Houston or Phoenix.
In addition, the following limitations of our study may be addressed in future studies. Because of the difficulty in determining what percentage of one urban category can be considered as ISA, the mean value according to the EPA definition was employed—that is, ISA of 35% in low urban density and 90% in high urban density. A better and more precise method is needed to quantify this percentage. Further, all of the ISA in this study was assumed to be effective impervious area (EIA), so the precipitation falling on it will become runoff directly, which may result in higher values of the simulated average R–B index and high-flow frequency than observed data at basin scale. Despite these limitations, the study results provide a significant and detectable effect of urbanization on hydrologic regime using a combination of observations and model results, which also demonstrate that the bulk parameterization approach of VIC model can yield significant improvements in the representation of urban areas.
We gratefully acknowledge the support of the NASA Land Cover Land Use Change program for this research.
Corresponding author address: Guoxiang Yang, Department of Agronomy, Purdue University, Lilly Hall, 915 W. State Street, West Lafayette, IN 47907. Email: email@example.com