Identifying critical source areas of nonpoint source pollution with SWAT and GWLF

Ecological Modelling 268 (2013) 123–133

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Identifying critical source areas of nonpoint source pollution with SWAT and GWLF Rewati Niraula a , Latif Kalin b,∗ , Puneet Srivastava c , Christopher J. Anderson b a b c

Hydrology and Water Resources, University of Arizona, Tucson, AZ 85721, United States School of Forestry and Wildlife Sciences, Auburn University, Auburn, AL 36849, United States Department of Biosystems Engineering, Auburn University, Auburn, AL 36849, United States

a r t i c l e

i n f o

Article history: Received 15 June 2013 Received in revised form 13 August 2013 Accepted 14 August 2013 Available online 14 September 2013 Keywords: SWAT GWLF Critical source area Water quality modelling Nonpoint source pollution

a b s t r a c t Identification of critical source areas (CSAs) (areas contributing most of the pollutants in a watershed) is important for cost-effective implementation of best management practices. Identification of such areas is often done through watershed modeling. Various watershed models are available for this purpose, however it is not clear if the choice (and complexity) of a model would lead to differences in locations of CSAs. The objective of this study was to use two models of different complexity for identifying CSAs. The relatively complex Soil and Water Assessment Tool (SWAT) and the simpler Generalized Watershed Loading Function (GWLF) were used to identify CSAs of sediment and nutrients in the Saugahatchee Creek watershed in east central Alabama. Models were calibrated and validated for streamflow, sediment, total nitrogen (TN) and total phosphorus (TP) at a monthly time scale. While both models performed well for streamflow, SWAT performed slightly better than GWLF for sediment, TN and TP. Sub-watersheds dominated by urban land use were among those producing the highest amount of sediment, TN and TP loads, and thus identified as CSAs. Sub-watersheds with some amount of agricultural crops were also identified as CSAs of TP and TN. A few hay/pasture dominated sub-watersheds were identified as CSAs of TN. The identified land use source areas were also supported by field collected water quality data. A combined index was used to identify the sub-watersheds (CSAs) that need to be targeted for overall reduction of sediment, TN and TP. While many CSAs identified by SWAT and GWLF were the same, some CSAs were different. Therefore, this study concludes that model choice will affect the location of some CSAs. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Approximately 67% of lakes, reservoirs and ponds, and 53% of rivers and streams in the U.S. are classified as impaired, needing immediate attention (USEPA, 2013). Impairment of water bodies due to elevated levels of nutrients and sediments originating from upland areas (i.e. watersheds) is a serious problem around the world. High level of nutrients can cause problems such as toxic algal blooms, oxygen deficiency, fish kills, and loss of biodiversity. These problems can also make the water unsuitable for drinking, industrial, agricultural and recreational use (Carpenter et al., 1998). Watershed management offers a strong basis for developing and implementing effective management strategies (such as

∗ Corresponding author at: 602 Duncan Dr., Auburn, AL 36849, USA. Tel.: +1 334 844 4671; fax: +1 334 844 1084. E-mail addresses: [email protected] (R. Niraula), [email protected], [email protected] (L. Kalin), [email protected] (P. Srivastava), [email protected] (C.J. Anderson). 0304-3800/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ecolmodel.2013.08.007

riparian zones, vegetation strips, retention ponds, etc.) to protect water resources (USEPA, 2003). Past efforts in reducing pollutant loads from watersheds have mainly focused on point sources and have failed to adequately address the impact of nonpoint sources. If nonpoint sources of pollutants are not addressed, water bodies can continue to be impaired (USEPA, 2003). However, nonpoint sources of nutrients and sediments are difficult to identify and control because they originate from spatially and temporally varying areas (Carpenter et al., 1998). The level of sediment and nutrient contribution from different parts of a watershed can vary substantially. Some typically small and well defined areas contribute much of the sediment, P, and N into the watershed outflow (Walter et al., 2000; Pionke et al., 2000) and over relatively short periods (Dillon and Molot, 1997; Heathwaite et al., 2005). But in many situations source areas are not well defined but diffused. Certain areas with a particular type of soil, land use/cover and slope are more vulnerable than the others in terms of nutrient and sediment loss. These areas are known as critical source areas (CSAs). It is extremely important to identify these sources of pollutants for cost-effective management practices.

124

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

Identifying nutrient and sediment loss prone areas in a watershed and concentrating management efforts in those areas have been recommended in numerous studies (e.g., Nonpoint Source Task Force, 1984; Tim et al., 1992; Zhou and Goa, 2008). Such areas can be identified through sub-watershed level water monitoring, simulation modeling, or both (Sharpley et al., 2002). Direct water monitoring and field studies are usually costly and labor intensive, and require a number of years of monitoring to sufficiently account for climatic fluctuations. The use of watershed models, such as Soil and Water Assessment Tool (SWAT) (Arnold et al., 1998) and Generalized Watershed Loading Function (GWLF) (Evans et al., 2002), can avoid most limitations associated with field studies and can help in identifying and prioritizing sub-watersheds for costeffective implementation of management practices (Tripathi et al., 2005; Ouyang et al., 2008; Georgas et al., 2009). GWLF has widely been used for estimating streamflow, nitrogen (N) and phosphorus (P) loadings (Swaney et al., 1996; Lee et al., 2000), hydrochemistry (Schneiderman et al., 2002) and also for assessing changes in streamflow (Chang, 2003; Wu et al., 2007) and water quality (Tu, 2009) under different land use scenarios. GWLF has also been used for identification of CSAs at sub-watershed level (Markel et al., 2006; Georgas et al., 2009). Similarly, SWAT has been used around the world for predicting streamflow, and sediment and nutrient loads from watersheds (Spruill et al., 2000; Kirsh et al., 2002; Veith et al., 2005; Srivastava et al., 2006; Jha et al., 2007; Niraula et al., 2012a,b). SWAT has also been used in several studies for identification and prioritization of CSAs of sediments and nutrients (Tripathi et al., 2005; Ouyang et al., 2008; White et al., 2009; Ghebremichael et al., 2009; Panagopoulos et al., 2011; Shang et al., 2012; Niraula et al., 2012a). Wanger et al. (2007) studied the impact of alternative water quality models, by comparing GWLF and SWAT, on pollutant loading for Total Maximum Daily Loadings (TMDL) development. The use of alternative water quality models resulted in differences in required sediment reduction. The SWAT model load estimates were consistently larger than loads from GWLF. However, it is not clear if the use of different models would lead to different CSAs that are significant enough for practical implementation of Best Management Practices (BMPs). In a related study, Niraula et al. (2012b) found that calibration of the SWAT model had very little effect on locations of nutrients and sediment CSAs. Therefore, the objective of this study was to assess the effect of model choice on CSA locations. Two models of different complexity, SWAT and GWLF, were used for this purpose. Both models were utilized to identify sediment and nutrients CSAs in the Saugahatchee Creek watershed in east central Alabama.

2. Methodology 2.1. Study area The 570 km2 Saugahatchee Creek watershed (Fig. 1), selected for this study, is a sub-watershed of the Lower Tallapoosa sub-basin in east central Alabama. The watershed, as determined using National Land Cover Data (NLCD, 2001), was comprised of 67.8% forest, 10.0% grassland, 11.7% agricultural land (hay/pasture and row crops) and 8.4% urban area (Fig. 1). Although most of the watershed lies in the Piedmont physiographic province, a small portion lies in the Coastal Plains. The Piedmont covers a transitional area between the mostly mountainous Appalachians in the northeast and the relatively flat Coastal Plains in the southeast Alabama. While the soils in the Piedmont are dominated by loam and sandy loam, soils in remaining coastal plains are sandy loam based on the STATSGO soil database. Elevation in the watershed varies from 103 to 255 m. The study area is characterized by hot summers and mild winters

with average temperatures of 26 ◦ C and 7 ◦ C, respectively. The long term annual average rainfall in the watershed is 1336 mm. Alabama Department of Environmental Management (ADEM) has identified two segments within the Saugahatchee Creek watershed as being impaired for nutrients and organic enrichment/dissolved oxygen (ADEM, 2008). The nutrient of concern in both of the tributaries is phosphorus. ADEM also recommended development of TMDLs for addressing water quality problems in this watershed. 2.2. Watershed models 2.2.1. Soil and Water Assessment Tool (SWAT) The SWAT is a semi-distributed model that was primarily developed to predict the impact of land management practices on water, sediment, and agricultural chemical yields in large complex watersheds over long periods of time (Neitsch et al., 2005). The model inputs consist of topography, soil properties, land use/cover type, weather/climate data, and land management practices. The watershed is sub-divided into sub-watersheds and each sub-watershed is further divided into hydrological response units (HRU) based on topography, land use, and soil (Neitsch et al., 2005). Surface runoff in each HRU was estimated using a modification of the Soil Conservation Service Curve Number (SCS-CN) method (USDA, 1972). In the curve number method, daily precipitation is partitioned between surface runoff and initial and continued abstractions as a function of antecedent soil moisture condition. The total sub-watershed discharge computed by SWAT includes runoff from its HRUs and subsurface flow including lateral flow and return flow. Flow in SWAT is routed through channels using either Muskingum routing method or variable storage coefficient method (Neitsch et al., 2005). The latter was used in this study. Erosion and sediment yield from each HRU are estimated based on the Modified Universal Soil Loss Equation (MUSLE) (Williams, 1975). Sediment is routed through channels using a modification of Bagnold’s sediment transport equation (Bagnold, 1977). This equation estimates sediment transport capacity as a function of flow velocity. The model either deposits or erodes sediment, depending on the sediment load entering the channel and the capacity of the flow. SWAT models nitrogen and phosphorus cycles through five different pools of nitrogen (two inorganic forms: NH4 + and NO3 − ; three organic forms: fresh, stable and active) and six different pools of phosphorus (three inorganic forms: solution, active and stable; three organic forms: fresh, stable and active) in soil (Neitsch et al., 2005). Mineralization, decomposition, and immobilization are important processes in both N and P cycles. Organic N and P transport with sediment is estimated using a loading function developed by McElroy et al. (1976) and later modified by Williams and Hann (1978). Daily organic N and P runoff losses are calculated by loading functions based on the concentrations of these elements in the top soil layer, the sediment yield, and an enrichment ratio. Nitrate concentration in mobile water is calculated and multiplied with mobile water volume to estimate total nitrate lost from the soil layer. Mobile water is the sum of runoff, lateral flow and percolation. The soluble P removed in runoff is estimated using the P concentration in the top soil layer, runoff volume and a P soil partitioning coefficient. Further details can be found in Neitsch et al. (2005). 2.2.2. Generalized Watershed Loading Function (GWLF) The GWLF model is a combined distributed/lumped parameter, continuous watershed model (Evans et al., 2002), which has the ability to simulate runoff, sediment, and nutrient (N and P) loads from various source areas, each of which is considered uniform with respect to soil and cover. GWLF uses land use, soil, and daily weather data for calculation of water balance. For estimation of sediment and nutrient loads, monthly calculations are made based

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

125

Fig. 1. Saugahatchee Creek watershed in east central Alabama. Also shown are land use/cover and USGS flow and NOAA weather gaging stations.

on the daily water balance aggregated to monthly values. It is considered a distributed model for surface loading because it allows multiple land uses within an area, but for each area other parameters are considered to be uniform by the model. This model does not consider spatial location of source areas, but simply adds the loads from different source areas within each sub-watershed. The model works as a lumped parameter model using a water balance approach for modeling sub-surface loading (Haith and Shoemaker, 1987; Haith et al., 1992). In GWLF, surface runoff is simulated using the SCS-CN method. Erosion and sediment yield is modeled using the Universal Soil

Loss Equation (USLE). GWLF simulates soil erosion by considering i) soil detachment by rainfall and ii) runoff transport relationships developed by Meyer and Wischeier (1969). A sediment delivery ratio, based on watershed size, and a transport capacity, based on average daily runoff, is then used to estimate sediment yield from each source area. In USLE equation while soil erodibility factor (K) depends on soil properties, other factors such as cover factor (C) and practice factor (P) depend on land use type. The daily runoff volume which transport sediment is calculated based on CN which is also a function of soil and land use/cover.

126

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

Table 1 Calibrated parameters of SWAT model with their default and calibrated values. Parameter

Descriptions

CN2 ESCO GWQMN SPEXP PRF SURLAG ADJ PKR PPERCO PHOSKD P-UPDIS PSP SOL LABP NPERCO SOL NO3

Initial SCS Runoff curve number for moisture condition II Soil evaporation compensation factor Threshold depth of water in shallow aquifer required for the return flow to occur Exponent parameter for calculating sediment reentrained in channel sediment routing Peak rate adjustment factor for sediment routing in the main channel Surface runoff lag time Peak rate adjustment factor for sediment routing in the subbasins(tributary) Phosphorus percolation coefficient Phosphorus soil partitioning coefficient Phosphorus uptake distribution factor Phosphorus sorption coefficient Initial soluble P concentration in surface soil layer (mg/kg) Nitrogen percolation coefficient Initial NO3 concentration in the soil(mg/kg)

a b

Range 25–92 0.01–1 0–5000 1–1.5 0–2 1–24 0.5–2 10–17.5 100–200 0–100 0.01–0.7 0–100 0–1 0–100

Default value

Calibrated value

Variesa 0.95 0 1 1 4 0.5 10 175 20 0.4 0 0.2 0

−5b 0.2 1200 1.5 1.97 1 2 17.5 200 80 0.7 11 1 21

Varies with land use and soil type. Value reduced by 5.

Nutrient loads from rural areas are transported in runoff water and eroded soil from numerous source areas. Dissolved loads are calculated by multiplying dissolved N and P coefficients with runoff. Solid phase nutrient loads are calculated by multiplying average sediment nutrient concentrations with monthly sediment yield. All N and P inputs from urban areas are assumed to be in solid phase. The model uses an exponential accumulation and wash-off function for estimating urban loadings. Sub-surface losses are calculated using dissolved N and P concentrations (which are model inputs) considering a single lumped-parameter contributing area (Evans et al., 2002). 2.3. Model inputs ArcSWAT 2.1 (Winchell et al., 2008) and AVGWLF 7.1 (Evans et al., 2008) were used to set up and develop models for the Saugahatchee Creek watershed. Data required in this study included Digital Elevation Model (DEM), soil properties, land use/cover, climate data such as precipitation, and minimum/maximum temperature, and point sources. A 10-m resolution DEM was used to delineate the watershed and sub-watershed boundaries, which were used in both models. State Soil Geographic (STATSGO) database was used to derive soil parameters. Land use/cover data were obtained from the National Land Cover Dataset (NLCD) for the year 2001 (Fig. 1). Point source discharges from three point sources (North Auburn and Opelika waste water treatment plants, and West Point Stevens) were included. Daily precipitation and minimum and maximum air temperature data between January 1995 and December 2008 were collected from three NOAA weather stations (Fig. 1). Flow data for the period from 2000 to 2008 were obtained from a USGS gauging station located within the watershed (Fig. 1). Water quality data for the same location (USGS station) was obtained from Alabama Department of Environmental Management (ADEM, 2002) and Department of Fisheries and Allied Aquaculture at Auburn University (DFAA, 2001). At least 3 data points per month for total suspended solids (TSS), nitrate (NO3 ), nitrite (NO2 ), Ammonia (NH3 + NH4 ), total nitrogen (TN), soluble reactive phosphorus (SRP) and total phosphorus (TP) including measurements immediately after some storm events, were available for the period 2000–2003. Organic P and N were calculated by subtracting the sum of mineral components from the TP and TN, respectively. 2.4. Calibration and validation of models Both SWAT and GWLF models were run from 1995 to 2008. The periods 2000–2004 and 2005–2008 were selected as the

calibration and validation periods, respectively, for flow. The first five years (1995–1999) were used as a warm up period to minimize uncertain initial conditions (e.g., soil moisture, groundwater level, ground residue, nutrient pool, etc.). A manual calibration technique was adopted. To identify the most sensitive parameters for calibration, a manual one-at-a-time (OAT) relative sensitivity analysis was carried out for the list of parameters synthetized from review of the current literature. Models were first calibrated for streamflow using data from the USGS gauge station (Fig. 1). For this purpose, observed streamflow was separated into surface runoff and baseflow components with a baseflow separation filter program (Arnold et al., 1995; Arnold and Allen, 1999). In addition to matching hydrographs of model simulated baseflow and excess runoff with observed counterparts, quantitative measures (percent bias, correlation coefficient, and Nash-Sutcliffe efficiency) were also used during calibration. Once the models were calibrated for flow, they were subsequently calibrated for sediment, TN and TP at the same USGS station. Due to lack of sufficient water quality data, monthly sediment was calibrated for year 2000 and validated for year 2002. High quality data were available only for those years. TN and TP were calibrated for the year 2000 (to be consistent with sediment calibration, which affects TN and TP) and validated for the period 2001–2002. The Load Estimator (LOADEST) software (USGS, 2013) developed by USGS was used to estimate monthly loads from spontaneous grab samples of sediment and nutrients, instantaneous streamflow, and daily average streamflow. The degree of fit (R2 ) between the observed instantaneous flow data and TSS, TN and TP used for estimating loads were 0.77, 0.72 and 0.76, respectively. Various hydrologic and water quality parameters (Tables 1 and 2) were changed within their ranges to get the best fit with the observed data. Four evaluation criteria: percent bias (PBIAS; Moriasi et al., 2007), Nash–Sutcliffe efficiency (NSE; Moriasi et al., 2007), coefficient of determination (R2 ; Niraula et al., 2012b) and line graphs were used to assess streamflow, sediment, TN and TP loads simulated by SWAT and GWLF. It should be noted that the calibrated model parameters are likely away from their true values because of the uncertainties coming from measurement errors in water quality data and the additional uncertainty exerted by the LOADEST model.

2.5. Identification of critical source areas (CSAs) CSAs were identified at the sub-watershed level. The sediment and nutrient yields from each sub-watershed were analyzed based on loads per unit area to identify the CSAs. Maps were created based on these loadings to depict the CSAs separately for sediment, TN

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133 Table 2 Calibrated parameters of GWLF model with their default and calibrated values. Parameter

Default value

Calibrated value

Curve number Recess coefficient Seepage coefficient Nitrogen in sediments (mg/kg) Phosphorus in sediments (mg/kg) Sediment A factor Erosivity coefficient January–February March April–December C factor: low/high intensity dev P factor: low/high intensity dev Nitrogen runoff coefficient Low intensity dev (kg/ha/day) High density dev (kg/ha/day) Phosphorus runoff coefficient Cropland(mg/l) Forest(mg/l) Low intensity dev (kg/ha/day) High density dev (kg/ha/day) N (mg/l) in groundwater P (mg/l) in groundwater

Variesa 0.1 0 3000 750 1.0997 × 10−4

+5b 0.01 0.006 7500 1500 1.0997 × 10−3

a b

0.18 0.28 0.18 0.02 0.02

0.18 0.4 0.06 0.4 0.6

0.012 0.101

0.02 0.15

0.079 0.006 0.002 0.011 1 0.01

0.6 0.036 0.007 0.015 0.4 0.02

Varies with land use and soil type. Value increased by 5.

and TP using simulation results from both SWAT and GWLF models. In order to do this, sub-watersheds were ranked in descending order based on yields. The sub-watershed with the highest yield was ranked first and the one with the lowest yield was ranked last. Moving from the highest ranking to the lowest, sub-watersheds that collectively contribute to 20% of the sediment, TP or TN were considered as CSAs. Different CSAs for sediment, TN and TP were obtained. The 20% threshold was chosen to provide results for comparison across two models. Actual threshold is a function of cost of implementing management practices. For instance, if budget is limited for implementation of management practices a lower threshold is needed. A combined index was also defined to identify the subwatersheds that can be considered as CSAs and need to be targeted for overall reduction of sediment, TN and TP. Although TMDL programs target reduction of individual stressors and, therefore identification of CSAs of individual stressors seems more reasonable, sometimes a water body can be on the 303(d) list for more than one pollutant. Since many BMPs can simultaneously reduce N, P and TSS loads (e.g. vegetative filter strips), with a single BMP combined load reduction can be achieved. The combined index can help identify those areas that are critical for multiple stressors. Having BMPs in those areas will be more economical. This index is given by Ij =



(ωi Yi,j )

(1)

where, Ij is the combined index for sub-watershed j, Yi,j is an index for sediment (i = 1), TN (i = 2) and TP (i = 3) for sub-watershed j, and is defined as Yi,j =

Ri,min − Ri,j Ri,min − Ri,max

(2)

where, Ri,j is the rank of watershed j for constituent i, and Ri,min and Ri,max are respectively the lowest and highest ranks for constituent i. Note that Ri,max = 1, Ri,min = total number of sub-watersheds and 0 ≤ Yi,j ≤ 1. For our study watershed Ri,min and Ri,max are 106 and 1, respectively. In Eq. (1), ωi is a subjectively chosen weight given to ωi = 1. For the purpose each Yi based on their importance, where of our study we assigned equal weights (1/3) to all three components (sediment, TN and TP). Theoretically Ij varies between 0 and

127

1. If the same sub-watershed is ranked last in all three constituents then Ij = 0. Similarly, if the same sub-watershed is ranked highest in all three, then Ij = 1. 3. Results and discussion 3.1. Calibration/validation of SWAT and GWLF (monthly time scale) Both SWAT and GWLF models were able to predict the monthly streamflow with good accuracy (Fig. 2a and Table 3). According to the performance statistics, SWAT and GWLF performed equally well with respect to all three measures. GWLF had slightly better NSE values and SWAT had marginally better PBIAS values, but for practical purposes the differences are insignificant. Both models were able to capture the months with high- and low-flows (Fig. 2a). Although both SWAT and GWLF models were able to predict the monthly sediment loadings with sufficient accuracy (Fig. 2b and Table 3), SWAT performed better than GWLF during both calibration and validation periods on the basis of NSE and R2 . However, SWAT and GWLF performed equally well based on PBIAS during both calibration and validation periods. SWAT overestimated sediment by 0.9% and 3.7% during the calibration and validation periods, respectively, while GWLF underestimated the monthly sediment loadings by 0.5% during the calibration period and overestimated by 4.0% during the validation period. Both models predicted the monthly TN loadings with sufficient accuracy (Fig. 2c and Table 3). SWAT performed slightly better than GWLF on the basis of NSE and R2 during both calibration and validation periods. However, SWAT performed markedly better based on PBIAS. SWAT underestimated TN loading by 2.4% during the calibration period and underestimated by 0.3% during the validation period. GWLF, on the other hand, overestimated TN loading by 7.8% during the calibration period and underestimated by 9.0% during the validation period. Even though both SWAT and GWLF models predicted monthly TP loadings well (Fig. 2d and Table 3), SWAT performed better than GWLF on the basis of performance statistics (NSE, R2 and PBIAS) during both calibration and validation periods. SWAT overestimated phosphorus loading by 7.0% during the calibration period and underestimated by only 0.1% during the validation period. GWLF overestimated TP by 10.2% and underestimated by 7.8% during the calibration and validation periods, respectively. 3.2. Critical source areas (CSAs) The calibrated SWAT and GWLF models were used to identify and compare the CSAs in the Saugahatchee Creek watershed at the sub-watershed level. Different CSAs were identified with respect to sediment, TN and TP loadings because the factors driving each are likely to be different, but not mutually exclusive. To eliminate the effect of differences in areas of the sub-watersheds, average of annual loadings per unit area (i.e. yield) for the 2000–2008 period were used to identify the CSAs. 3.2.1. CSAs of sediment Both models produced somewhat similar areas as CSAs of sediment (Fig. 3a and d). Based on 20% contribution, 5 sub-watersheds, covering 4.2% of the watershed, were identified as CSAs by GWLF. Similarly, 6 sub-watersheds, covering 4.6% of the watershed, were identified as CSAs by SWAT. Although the ranks are not exactly in the same order, sub-watersheds 25, 53, 56, 69 and 73 were identified as CSAs by both SWAT and GWLF models. Sub-watershed 23 was also identified as a CSA of sediment by SWAT. While average sediment yield ranged from as high as 9.77 tons/ha/yr (t/ha/y) to as low as 0.06 t/ha/y based on the SWAT results (Fig. 3a), it ranged

128

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

Fig. 2. Observed, SWAT-simulated, and GWLF-simulated monthly (a) streamflow, (b) sediment load, (c) TN load, and (d) TP load, for the calibration and validation periods.

between 0.07 t/ha/y and 6.15 t/ha/y based on GWLF results (Fig. 3d). This showed that sediment yields obtained by SWAT showed larger variation than GWLF generated sediment yields. The average basin sediment yield was found to be 0.94 t/ha/y with GWLF compared to 1.44 t/ha/y with SWAT. 3.2.2. CSAs of TP In spite of the fact that most of the CSAs of TP identified by each model overlap, several sub-watersheds were identified as CSAs by one model and not by the other. While SWAT identified 7 sub-watersheds (5.8% of total watershed area) as CSAs for TP, GWLF identified 11 sub-watersheds (9.6% of total watershed area) as CSAs. Although the ranks are not exactly in the same order, subwatersheds 25, 51, 53, 56, 59, 69, and 73 were identified as CSAs of TP by both the SWAT and the GWLF models. Most of these subwatersheds were also identified as sediment CSAs, indicating that TP load is highly correlated with sediment load. Sub-watersheds

64, 79, 88 and 97 were also identified as CSAs of TP by GWLF, but not by SWAT. TP yields ranged from 0.02 kg/ha/y to 0.87 kg/ha/y according to SWAT results (Fig. 3b), and between 0.06 kg/ha/y and 0.71 kg/ha/y according to GWLF results (Fig. 3e). This again showed that SWAT generated outputs showed slightly wider range than GWLF generated outputs. The basin-average TP yield was found to be 0.22 kg/ha/y with GWLF, which is slightly higher than the SWAT generated average of 0.19 kg/ha/y. Although the average TP load from SWAT was slightly smaller than the GWLF average, the fact that highest TP load from all the sub-watersheds was obtained with SWAT points to the strong correlation between P and sediment. 3.2.3. CSAs of TN In the case of TN, while 10 sub-watersheds (10.5% of total watershed area) were identified as CSAs by SWAT, 13 sub-watersheds (13% of total watershed area) were considered as CSAs by GWLF. Both models identified sub-watersheds 4, 56, 69, 79, 88, 97 and

Table 3 Performance measures of the SWAT and the GWLF during calibration and validation periods for flow, sediment, TN and TP. Parameter

Model

Calibration

Validation 2

NSE

R

PBIAS (%)

NSE

R2

PBIAS (%)

Streamflow

SWAT GWLF

0.90 0.91

0.90 0.92

+3.1 −6.7

0.74 0.79

0.78 0.82

−1.4 +3.5

Sediment

SWAT GWLF

0.72 0.68

0.72 0.68

+0.9 −0.5

0.86 0.78

0.87 0.81

+3.7 +4.0

TN

SWAT GWLF

0.75 0.70

0.81 0.78

−2.4 +7.8

0.90 0.87

0.92 0.91

−0.3 −9.0

TP

SWAT GWLF

0.85 0.65

0.89 0.79

+7.0 +10.2

0.88 0.87

0.91 0.88

−0.1 −7.8

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

129

Fig. 3. Sediment, TP and TN yields from each sub-watershed as estimated by SWAT (a, b, c, respectively) and GWLF (d, e, f, respectively). Numbered sub-watersheds are identified as CSAs.

103 as CSAs. Sub-watersheds 25, 73 and 102 were also identified as CSAs based on SWAT, but not by GWLF. Similarly, sub-watersheds 1, 11, 14, 53, 59 and 85 were identified as CSAs by GWLF, but not by SWAT. While TN yields ranged between 0.57 kg/ha/y and 5.31 kg/ha/y based on SWAT results (Fig. 3c), it ranged between 0.85 kg/ha/y and 4.73 kg/ha/y according to GWLF (Fig. 3f). Again, SWAT had slightly more variation in TN yields. The average basin TN yield was found to be 1.93 kg/ha/y with GWLF compared to 2.1 kg/ha/y with SWAT. The differences in values obtained by SWAT and GWLF can be related to the way NO3 pool is handled. NO3 is mostly transported through subsurface flow in watersheds as was the case with SWAT in this study too. However, with GWLF, a dissolved N coefficient is simply used to simulate subsurface N, which might not be sufficient to capture the NO3 -N loads from the watershed.

3.3. Combined CSA Index Combined CSA index Ij was determined for each sub-watershed j by assigning the same weights to sediment, TN, and TP (ωi = 1/3 in Eq. (2)). Sub-watersheds with Ij ≥ 0.9 were subjectively identified as CSAs collectively for all parameters (sediment, TN and TP). Again, the purpose here was to present the idea. The threshold 0.9 was chosen so that we deal with a manageable number of CSAs. A map (Fig. 4) was produced to show these CSAs based on results from both SWAT and GWLF models to visualize the differences. Sub-watersheds 25, 53, 56, 59, 69 and 73 had Ij ≥ 0.9 based on both models and thus were identified as CSAs by both models. However, sub-watersheds 23 and 51 had Ij ≥ 0.9 only based on SWAT results and thus were also identified as CSAs based on SWAT. Likewise, sub-watershed 88 had Ij ≥ 0.9 only based on GWLF results

130

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

Fig. 4. Critical source areas based on the combined index (Eq. (1)). The numbered sub-watersheds should be first targeted for management practices for overall reduction of sediment, TP and TN loads.

and thus was identified as CSA by the GWLF model. CSAs identified by SWAT covered only 6.5% of the watershed area and contributed 26.5% of the sediment, 23.1% of TP and 13.9% of TN loadings. Similarly, CSAs determined by GWLF covered 5.6% of the watershed and were responsible for contributing 23.1% of sediment, 16.5% of TP and 12.7% of TN. While sub-watersheds 23, 25, 51, 53, 56, 59, 69 and 73 were all dominated by urban land use with considerable amount of pasture land, sub-watershed 88 was comprised of some cropland and lay predominantly on the coastal plain soils. It was clear from this study that, sub-watersheds dominated by urban area were among those producing highest amount of sediment and nutrient loads and thus were identified as CSAs. Also sub-watersheds with some amount of agricultural crops were identified as CSAs of TP and TN. A few hay/pasture dominated sub-watersheds were identified as CSAs of TN (e.g. sub-watershed 4, which is 46% hay/pasture). Table 4 provides the land use composition of the combined CSAs only.

3.4. Differences in model performances and identified CSAs As mentioned earlier, both SWAT and GWLF performed equally well for streamflow. This could be because both models use the same method (SCS-CN) for estimating runoff. However, in the case of sediment, TN and TP, SWAT performed slightly better than GWLF. In the case of sediment, there was also a difference in one of the CSAs identified by the two models. While GWLF uses the conventional USLE method for estimating soil erosion, SWAT uses the MUSLE equation. USLE uses rainfall intensity as the erosive energy, whereas MUSLE uses runoff volume and peak flow rate to simulate sediment erosion and yield. This improves the prediction accuracy and also eliminates the need for delivery ratio. In GWLF, nutrient loads from rural areas are calculated based on dissolved N and P coefficients in surface runoff, and a sediment coefficient in sediment load. Nutrient loads form urban areas are estimated through exponential accumulation and wash-off function. On the other hand, SWAT models N and P cycles comprehensively. Processes like mineralization, decomposition, and immobilization are allowed to take place in the soil in each HRU. Thus, SWAT provides a more mechanistic and process-based approach than GWLF. As a result, SWAT predicted sediment, TN and TP loads better than GWLF did (see Table 3). Further, because SWAT and GWLF conceptualizes sediment, TN and TP processes differently, there were some variations in the locations of identified CSAs. Also, since SWAT divides the subwatersheds into smaller computational units, i.e. HRUs, there is a greater chance that it can distinguish sub-watersheds with higher

loadings than the others. The wide range of output values from different sub-watersheds for sediment, TN and TP obtained with SWAT further support this.

3.5. Implications of not choosing the right model Models are simplified mathematical representations of real watershed systems, and thus are never perfect (Gupta et al., 1998). Model structural adequacy can vary with “engineering” view point (focused on functional adequacy and usually take a decision making perspective); “physical science” viewpoint (consistency with the physical system) and “system science” viewpoint (data-based hybrid of the first two which stresses on both consistency and principle of parsimony) (Gupta et al., 2012). Thus, choice of the right model depends on number of factors. Water quality models are typically chosen for a specific application based on the model’s capability to simulate the dominant processes and desired output, professional and personal preference, software availability, prior experience in using the model, and computational cost that include both time and financial resources (Wanger et al., 2007). Results from this study showed that not choosing the right model may have important implications. For instance, GWLF identified 4 extra sub-watersheds as CSAs for TP compared to SWAT. They had 7 sub-watersheds in common as CSAs. So, if GWLF is used as the base in deciding where to implement BMPs, about 65% more area should be targeted compared to SWAT. This might have extremely important economic implications. The total areas targeted for BMPs implementation in conservation programs are generally the limiting factor due to limited funding. The differences in locations of CSAs were more evident with TN. Out of the 13 CSAs identified by GWLF only 7 were also recognized as CSA by the SWAT model. That means GWLF has 6 sub-watersheds not identified as CSA by SWAT. Note that we are not promoting one model over another, rather pointing out the differences in CSAs due to the use of two different models. Since no measured data were available at sub-watershed level (we had 106 sub-watersheds), we cannot establish which model was identifying the locations of CSAs more properly, which is a limitation of this study. Although we did not have the water quality data to verify the CSAs identified by the models, some unpublished water quality data exists from several first and second order creeks around the City of Auburn (Fig. 1) to support most of the land use source areas identified by the models. Table 5 summarizes these water quality

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

131

Table 4 Land use composition of combined CSA’s in the Saugahatchee Creek watershed. Sub. #

Water (%)

Urban (%)

Forest (%)

Shrubland (%)

Hay/pasture (%)

Cropland (%)

Identified bya

23 25 51 53 56 59 69 73 88

0.4 0.7 1.1 1.1 0.0 0.2 0.9 1.6 1.8

43.4 76.1 54.3 64.1 82.0 70.4 69.1 54.7 2.8

43.1 20.6 41.8 31.0 13.2 20.8 25.3 34.2 60.0

4.4 0.2 2.2 1.1 1.4 2.2 1.9 3.5 25.8

8.9 2.4 0.6 2.7 3.4 6.4 0.3 5.2 3.3

0.0 0.0 0.1 0.1 0.0 0.0 2.5 0.8 6.3

S S, G S S, G S, G S, G S, G S, G G

a

S: SWAT G: GWLF

Table 5 Mean (±SE) total suspended solid (TSS), TN, and TP concentrations (mg L−1 ) for urban and reference creeks in the Saugahatchee Creek basin sampled during 2009–2010 (unpublished data). Creek

TSS Ashton Camden Shadowwood Tuscanny Webster Urban Avg. Rattlesnake Clara Ref. Avg. TN Ashton Camden Shadowwood Tuscanny Webster Urban Avg. Rattlesnake Clara Ref. Avg. TP Ashton Camden Shadowwood Tuscanny Webster Urban Avg. Rattlesnake Clara Ref. Avg.

No. of samples

BF

BF

FL

RL

9 9 9 9 9

7 3 3 5 6

3 6 7 4 4

9 9

4 4

5 6

8 8 8 8 8

6 2 1 4 4

2 6 6 4 4

8 8

3 3

6 5

8 8 8 8 8

6 2 1 4 4

2 6 6 4 4

8 8

3 3

6 6

12.2 (5.6) 4.4 (0.5) 56.1 (42.1) 9.5 (2.0) 8.4 (3.0) 18.1 (9.6) 7.4 (2.3) 3.8 (1.2) 5.6 (1.8)

Storm flow FL

RL

Average (FL + RL)/2

14.6 (6.5) 14.5 (10.4) 22.5 (8.6) 27.8 (7.8) 22.0 (4.7) 20.3 (2.6) 5.9 (4.4) 2.7 (1.9) 4.3 (1.6)

60.9 (31.2) 51.5 (18.9) 199.7 (34.4) 41.8 (19.8) 25.5 (4.3) 75.9 (31.5) 152.5 (15.5) 184.8 (15.1) 168.7 (16.1)

37.8 33.0 111.1 34.8 23.8 48.1 (15.9) 79.2 93.4 86.5 (7.1)

0.79 (0.07) 0.35 (0.03) 1.07 (0.27) 0.39 (0.09) 0.87 (0.16) 0.69 (0.14) 0.27 (0.06) 0.21 (0.04) 0.24 (0.03)

0.88 (0.07) 0.98 (0.32) 1.02 0.43 (0.15) 0.73 (0.08) 0.81 (0.11) 0.40 (0.12) 0.27 (0.12) 0.34 (0.07)

1.63 (0.53) 0.78 (0.18) 1.25 (0.31) 0.64 (0.23) 0.73 (0.07) 1.01 (0.19) 0.55 (0.16) 0.41 (0.18) 0.48 (0.07)

1.26 0.88 1.13 0.53 0.73 0.91 (0.13) 0.48 0.34 0.41 (0.07)

0.053 (0.016) 0.017 (0.005) 0.066 (0.030) 0.026 (0.005) 0.119 (0.080) 0.056 (0.018) 0.027 (0.007) 0.018 (0.003) 0.023 (0.005)

0.070 (0.015) 0.093 (0.008) 0.059 0.044 (0.019) 0.051 (0.018) 0.063 (0.009) 0.038 (0.026) 0.029 (0.019) 0.034 (0.004)

0.114 (0.041) 0.067 (0.018) 0.040 (0.009) 0.065 (0.039) 0.034 (0.015) 0.064 (0.014) 0.039 (0.015) 0.043 (0.013) 0.041 (0.002)

0.092 0.080 0.050 0.055 0.043 0.064 (0.009) 0.039 0.036 0.037 (0.001)

SE: standard error, BF: baseflow, FL: falling limb, RL: rising limb.

data and generally supports that urban areas represented CSAs of nutrients and sediment into the Saugahatchee Creek. Water quality from urban areas around the City of Auburn showed mixed results for TSS but consistently higher concentrations of TN and TP compared to reference streams (streams draining forested watersheds). Surprisingly, some of the highest TSS concentrations are found in the reference streams but only during the rising limb of stormflow. During the baseflow and falling limb, the average urban creek had 3–5 times greater TSS concentration than the reference creeks. Nutrient concentrations (N and P) were consistently higher in urban creeks compared to the reference creeks (often near twice the concentration) for baseflow, falling limb, and rising limb. These data supported indications provided by the two models that urban land generated higher levels of nutrients and possibly sediment. 4. Summary and conclusions Two watershed models differing in their conceptual structures and complexities, SWAT and GWLF, were set up, calibrated, and

validated in an east-central Alabama watershed. The models were then utilized to identify critical source areas (CSAs) of sediment, TN and TP for implementation of cost effective management practices in the watershed. Based on the overall model performance statistics, it was observed that SWAT performed slightly better than GWLF. Although performance of GWLF was similar to SWAT for streamflow during both calibration and validation periods, SWAT performed better for sediment, TN, and TP. The calibrated and validated models were used to identify the CSAs in the study watershed at sub-watershed level based on loadings per unit area (yield). In general, sub-watersheds dominated by urban areas were among those producing the highest amount of sediment, TN and TP yields, and thus were identified as CSAs. However, sub-watersheds with some amount of agricultural crops were also identified as CSAs of TP and TN. A few hay/pasture dominated sub-watersheds were identified as CSAs of TN. This study revealed that CSAs can vary based on the parameters of interest (sediment, TN or TP). Based on a combined index, while 8

132

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133

sub-watersheds were identified as CSAs by the SWAT model, 7 subwatersheds were identified as CSAs by the GWLF, among which 6 sub-watersheds were identical. Therefore, although there were similarities in many of the identified CSAs, not all the CSAs identified by the models were mutual. GWLF did not recognize two sub-watersheds, which were identified as CSAs by the SWAT model. Similarly, SWAT did not detect one of the sub-watersheds identified as CSA by GWLF. Dissimilarities in CSAs are attributed to the differences in model conceptualizations implemented in the SWAT and GWLF models. Although SWAT performed slightly better than GWLF for predicting flow, sediment, TN, and TP, we cannot conclude that SWAT was better for identifying and locating the CSAs, since there was no measured data to verify the CSAs identified by each model. Unpublished water quality data from several first and seconds order tributaries generally agreed with the modeled indications of CSAs in the Saugahatchee Creek watershed. We found that urban creeks in the basin normally had higher concentrations of sediment, N, and P although considerable variability (particularly regarding sediment) was noted. Targeted and more intensive monitoring of creek sediment may shed light on the dynamics and variability found among creeks with different land uses. This study was conducted in a forest dominated watershed with significant urban portion concentrated at the upstream headwaters. It would be interesting to carry out similar studies in future in intensely cultivated watersheds and/or more heterogeneous watersheds representative of the majority of NPS pollution sources. References Arnold, J.G., Allen, P.M., Muttiah, R., Bernhardt, G., 1995. Automated base flow separation and recession analysis techniques. Ground Water 33 (6), 1010–1018. Arnold, J.G., Srinivasan, R., Muttiah, R.S., Williams, J.R., 1998. Large area hydrologic modeling and assessment, Part I: model development. J. Am. Water Resour. Assoc. 34 (1), 73–89. Arnold, J.G., Allen, P.M., 1999. Automated methods for estimating baseflow and groundwater recharge from streamflow records. J. Am. Water Resour. Assoc. 35 (2), 411–424. ADEM, 2002. Surface water Quality Screening Assessment of the Tallapoosa River Basin-2000. Aquatic Assessment Unit, ADEM. September 2002. ADEM, 2008. Final Saugahatchee Creek Watershed Total Maximum Daily Load (TMDL) – Nutrients & OE/DO. Alabama Department of Environmental Management, Water Quality Branch, EPA Region 4. Bagnold, R.A., 1977. Bedload transport in natural rivers. Water Resour. Res. 13, 303–312. Carpenter, S.R., Caraco, N.F., Correll, D.L., Howarth, R.W., Sharpley, A.N., Smith, V.H., 1998. Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecol. Appl. 8 (3), 559–568. Chang, H., 2003. Basin hydrological response to changes in climate and landuse: the Constoga River Basin, Pennsylvania. Phys. Geogr. 24 (3), 222–247. DFAA, 2001. Water quality, loading and bio-assessment of Saughatchee Creek, 2000–2001. Department of Fisheries and Allied Aquaculture, Auburn University, December 2001. Dillon, P.J., Molot, L.A., 1997. Effect of landscape form on export of dissolved organic carbon, iron, and phosphorus from forested stream catchments. Water Resour. Res. 33 (11), 2591–2600. Evans, B.M., Lehning, D.W., Corradini, K.J., Petersen, G.W., Nizeyimana, E., Hamlett, J.M., Robillard, P.D., Day, R.L., 2002. A comprehensive GIS-based modeling approach for predicting nutrient loads in watersheds. J. Spat. Hydrol. 2 (2), 1–18. Evans, B.M., Lehning, D.W., Corradini, K.J., 2008. AVGWLF Version 7.1. User’s Guide. Penn State Institutes of Energy and the Environment, The Pennsylvania State University, University Park, PA. Ghebremichael, L.T., Watzin, M.T., Veith, T.L., 2009. SWAT modeling of critical source areas of runoff and phosphorus loss: Lake Champlain Basin, VT. In: 2009 International SWAT Conference Proceedings, Texas Water Resources Institute Technical Report No. 356, pp. 385–393. Georgas, N., Rangarajan, S., Farley, K.J., Jugupilla, S.R.K., 2009. AVGWLF-based estimation of the nonpoint sources nitrogen loads generated within long island sound sub-watersheds. J. Am. Water Resour. Assoc. 45 (3), 715–733. Gupta, H.V., Sorooshian, S., Yapo, P.O., 1998. Toward improved calibration of hydrologic models: multiple and non-commensurable measures of information. Water Resour. Res. 34, 751–763. Gupta, G.V., Clack, M.P., Vrugrt, J.A., Abramowitz, G., Ye, M., 2012. Towards a comprehensive assessment of model structural accuracy. Water Resour. Res. 48, W08301, http://dx.doi.org/10.1029/2011WR011044. Haith, D.A., Shoemaker, L.L., 1987. Generalized watershed loading functions for stream-flow nutrients. Water Resour. Bull. 23 (3), 471–478.

Haith, D.R., Mande, l, R., Wu, R.S., 1992. GWLF: Generalized Watershed Loading Functions User’s Manual, Version 2.0. Cornell University, Ithaca, NY. Heathwaite, A.L., Quinn, P.F., Hewett, C.J.M., 2005. Modeling and managing critical source areas of diffuse pollution from agricultural land using flow connectivity simulation. J. Hydrol. 304 (1–4), 446–461. Jha, M.K., Gassman, P.W., Arnold, J.G., 2007. Water quality modeling for the Raccoon River watershed using SWAT. Trans. ASAE 50 (2), 479–493. Kirsh, J., Kirsh, A., Arnold, J.G., 2002. Predicting sediment and phosphorus loads in the Rock River basin using SWAT. Trans. ASAE 45 (6), 1757–1769. Lee, K.Y., Fisher, T.R., Jordan, T.E., Correll, D.L., Weller, D.E., 2000. Modeling the hydrochemistry of the Choptank River basin using GWLF and GIS. Biogeochemistry 49, 143–173. Markel, D., Somma, F., Evans, B.M., Tyson, J., 2006. Using a GIS transfer model to evaluate pollutant loads in Lake Kinneret watershed, Israel. Water Sci. Technol. 53 (10), 75–82. McElroy, A.D., Chui, S.Y., Nebgen, J.W., Aleti, A., Bennet, F.W., 1976. Loading Functions for Assessment of Water Pollution from Non Point Sources. EPA Document EPA 600/2-76-151. USEPA, Athens, USA. Meyer, L.D., Wischeier, W.H., 1969. Mathematical simulation of the process of soil erosion by water. Trans. ASAE 12, 754–762. Moriasi, D.N., Arnold, J.G., Van Liew, M.W., Bingner, R.L., Harmel, R.D., Veith, T.L., 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50 (3), 885–900. Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., King, K.W., 2005. Soil and Water Assessment Tool: Theoretical Documentation, version 2005 (available at http://www.brc.tamus.edu/swat/). Niraula, R., Norman, L.M., Meixner, T., Callegary, J., 2012a. Multi-gauge calibration for modeling the Semi-Arid Santa Cruz watershed in Arizona–Mexico border area using SWAT. Air Soil Water Res. 2012, 41–57, http://dx.doi.org/10.4137/ASWR.S9410. Niraula, R., Kalin, L., Wang, R., Srivastava, P., 2012b. Determining nutrient and sediment critical source areas with SWAT model: effect of lumped calibration. Trans. ASABE 55 (1), 137–147. NLCD, 2001. National Land Cover Data 2001. Downloaded from alabamaview.org on September 2008. Nonpoint Source Task Force, 1984. National Nonpoint Source Policy. USEPA, Washington, DC. Ouyang, W., Hao, F.H., Wang, X.L., 2008. Regional point sources organic pollution modeling and critical areas identification for watershed best environmental management. Water Air Soil Pollut. 187, 251–261. Panagopoulos, Y., Makropoulos, C., Baltas, E., Mimkou, M., 2011. SWAT parameterization for the identification of critical diffuse pollution source areas under data limitations. Ecol. Model. 222 (19), 3500–3512. Pionke, H.B., Gburek, W.J., Sharpley, A.N., 2000. Critical source area controls on water quality in an agricultural watershed located in the Chesapeake basin. Ecol. Eng. 14 (4), 325–335. Schneiderman, E.M., Pierson, D.C., Lounsbury, D.G., Zion, M.S., 2002. Modeling the hydro-chemistry of the Cannonsville watershed with Generalized Watershed Loading Functions (GWLF). J. Am. Water Resour. Assoc. 38 (5), 1323–1347. Sharpley, A.N., Kleinman, P.J.A., McDowell, R.W., Gitau, M., Bryant, R.B., 2002. Modeling phosphorus transport in agricultural watersheds: processes and possibilities. J. Soil Water Conserv. 57 (6), 425–439. Spruill, C.A., Workman, S.R., Taraba, J.L., 2000. Simulation of daily and monthly stream discharge from small watersheds using the SWAT model. Trans. ASAE 43 (6), 1431–1439. Srivastava, P., McNair, J.N., Johnson, T.E., 2006. Comparison of process-based and Artificial Neural Network approaches for streamflow modeling in an agricultural watershed. J. Am. Water Resour. Assoc. 42 (3), 545–563. Swaney, D.P., Sherman, D., Howarth, R.W., 1996. Modeling water, sediment and organic carbon discharges in the Hudson–Mohawk basin. Estuaries 19 (4), 833–847. Shang, X., Wang, X., Zhang, D., Chen, W., Chen, X., Kong, H., 2012. An improved SWAT-based computational freamework for indtifying critical source areas for agriclutiurla pollution of the lake basin scale. Ecol. Model. 226, 1–10. Tim, U.S., Mostaghimi, S., Shanholtz, V.O., 1992. Identification of critical nonpoint source areas using Geographic Information System and water quality modeling. Water Resour. Bull. 28 (5), 877–887. Tripathi, M.P., Panda, R.K., Raghuwanshi, N.S., 2005. Development of effective management plan for critical sub-watersheds using SWAT model. Hydrol. Process. 19, 809–826. Tu, J., 2009. Combined impact of climate and land use changes on streamflow and water quality in eastern Massachusetts, USA. J. Hydrol. 379 (3–4), 268–283. USDA Soil Conservation Service, 1972. National Engineering Handbook. U.S. Government Printing Office, Washington, DC, Hydrology Section 4 (chapters 4–10). USEPA, 2003. National management measures for the control of nonpoint pollution from agriculture. EPA-841-B-03-004. Available at: http://water.epa.gov/polwaste/nps/agriculture/agmm index.cfm USEPA, 2013, National summary of water quality assessments of each waterbody type in US. Available at: http://ofmpub.epa.gov/waters10/attains nation cy. control#prob surv states USGS, 2013. Load Estimator (LOADEST): A Program for Estimating Constituent Loads in Streams and Rivers. http://water.usgs.gov/software/loadest Veith, T.L., Sharpley, A.N., Weld, J.L., Gburek, W.J., 2005. Comparison of measured and stimulated phosphorus with indexed site vulnerability. Trans. ASAE 48 (2), 557–565.

R. Niraula et al. / Ecological Modelling 268 (2013) 123–133 Walter, M.T., Walter, M.F., Brooks, E.S., Steenhuis, T.S., Boll, J., Weiler, K., 2000. Hydrologically sensitive areas: variable source area hydrology implications for water quality risk assessment. J. Soil Water Conserv. 55 (3), 277–284. Wanger, R.C., Dillaha, T.A., Yagow, G., 2007. An assessment of the reference watershed approach for TMDL with biological impairments. Water Air Soil Pollut. 181, 341–354. White, M.J., Storm, D.E., Busteed, P.R., Stoodley, S.H., Phillips, S.J., 2009. Evaluating nonpoint sources critical areas contributions at the watershed level. J. Environ. Qual. 38, 1654–1663. Williams, J.R., 1975. Sediment routing for agricultural watersheds. Water Resour. Bull. 11 (5), 965–974.

133

Williams, J.R., Hann, R.W., 1978. Optimal Operation of Large Agricultural Watersheds with Water Quality Constraints. Technical Report No. 96. Texas Water Resources Institute, Texas A&M University. Winchell, M., Srinivasan, R., Di Luzio, M., Arnold, J.G., 2005. 2008. ArcSWAT 2.1 interface for SWAT. User’s Guide. Grassland, Soil and Water Research Laboratory. USDA Agricultural Research Sevice Temple, Texas. Wu, W., Hall, C.A.S., Scetena, F.N., 2007. Modeling the impact of recent landuse changes on the stream flows in northeastern Puerto Rico. Hydrol. Process. 21, 2944–2956. Zhou, H.P., Goa, C., 2008. Identifying critical source areas for non-point phosphorus loss in Chaohu watershed. Huan Jing Ke Xue 29 (10), 2696–2702.