Long-term effects of land-use change on water resources in urbanizing watersheds

The changes in energy balance resulting from land-use change may significantly affect the amount and timing of water loss to the atmosphere as evapotranspiration (ET). Also, these will impact water fluxes in the watershed system, influencing runoff rate, flow volume, intensity, and frequency of floods. During the past century, land-use change in the SuAsCo (Sud-bury-Assabet and Concord) watershed has altered basin hydrology, sediment, and nutrient load that is detrimental to water resources in SuAsCo. This study uses an integrated physically-based model Hydrological Simulation Program-FORTRAN (HSPF), along with Land Transformation Model (LTM), to assess predicted temporal and spatial changes in water, nutrient, and sediment yields for future land-use scenarios of 2035, 2065, and 2100. Results showed that a 75% increase in effective impervious area and a 50% decrease in forest area in 2100 (from 2005 baseline levels) are projected to cause a 3% increase in annual stream-flow and a 69% increase in total yearly mean surface runoff. The average annual total suspended solid (TSS) yield at the watershed outlet is estimated to increase by 54% in 2100. An increase of 12% and 13% concentrations of average annual total phosphorus (TP) and total nitrogen (TN) are predicted by 2100 due to urban expansion and increased runoff volume. This integrated modeling approach will inform watershed managers and landowners about critical areas of the SuAsCo watershed to apply best management practices (BMPs) to mitigate the effects of land-use land cover (LULC) change.


Introduction
The global expansion of agricultural and urban areas, along with significant increases in energy, water, and fertilizer consumption, led to changes in hydrological processes and tremendous losses to biodiversity [1][2][3][4][5].Land-use change directly influences hydrological processes, such as ET [6,7], infiltration [8], and runoff [9,10].For instance, land-use conversion from forest to agricultural or urban would typically be associated with increased runoff with constant precipitation.This is because transpiration rates and ET in farmlands are generally lower than in the forest [11,12].Moreover, infiltration through impervious urban land is lower [13][14][15].Studies also suggested that land-use change is directly responsible for a 0.08 mm/year increase in global runoff [16,17].on the State's List of Impaired Waters (303d) under the Clean Water Act [58].Nutrients enter these rivers from nonpoint sources carried by stormwater runoff and point sources such as discharge from wastewater treatment plants.In addition, during the past century, land use changes in the SuAsCo watershed have altered basin hydrology, sediment, and nutrient load that is detrimental to the ecology and societal significance of water resources in SuAsCo [26,59].
This study is unique in applying a dynamic and long-term approach to understanding how water quantity and quality change spatially and temporally under future land use change, which is dominated by urbanizing processes in the SuAsCo watershed.An integrated system of a macroscale, physically based model HSPF, and machine learning ANN-based LTM model is used.HSPF is calibrated and validated to simulate SuAsCo watershed processes, and the LTM is used to generate future land-use scenarios.The statistical functional relationship built by ANN captured the pattern of the LULC change.This study also demonstrates how GIS and state-of-the-art ANN tools can be applied to predict LULC change.The purpose of this study is to assess the impacts of land-use change on water resources in an urbanizing watershed system, and the objectives include 1) to quantify the effects of land use on runoff, 2) to evaluate the impacts of future land use on total suspended sediments (TSS), and nutrients such as total nitrogen (TN), and (TP).In addition, this dynamic hydrologic modeling approach at a regional scale will provide an essential methodology for further research on urban expansion and its impacts on hydrology.
Null hypotheses of this study are i) Impacts of future land-use change on water flows are insignificant means <10% change from baseline, ii) Impacts of future land-use change on the concentration of TP, TN, and TSS in water are insignificant means <10% change from baseline.Specific alternative hypotheses tested in this study are (i) Impacts of future land-use change on water flows are significant; Ha: ΔW Quan /ΔC = Significant; where H a is an alternate hypothesis, ΔW Quan = change in water quantity such as a change in surface runoff water, storms or low flows, ΔC = land-use change (i.e., conversion of one land category to another category), significant means >10% change from baseline; and (ii) Impacts of future LULC on water quality are significant; H a : ΔW Qual /ΔC = Significant); where Ha is an alternate hypothesis, ΔW Qual = change in water quality such as a change in concentration of TP, TN, and TSS in water, significant means >10% change from baseline.

Study area
The SuAsCo (Sudbury, Assabet, and Concord Rivers) is a small, semi-urban watershed in eastern Massachusetts, USA.It is roughly 1.61 miles west of the Boston metropolitan area and is one of the 27 major watersheds in Massachusetts.The total drainage area of SuAsCo is 1012.69km 2 .The SuAsCo watershed encompasses partially or wholly 36 Massachusetts towns [59].The Lower Concord River Basin drains directly into the Concord River, formed at the convergence of the Sudbury and Assabet Rivers.Sudbury River Basin covers 419.6 km 2 , about 44% of the total SuAsCo basin, the Assabet River Basin 458.4 km 2 is about 41% of the SuAsCo Basin, and the Lower Concord River Basin covers 155.4 km 2 (Fig 1).The mean annual streamflow from the watershed outlet at NWIS gaging station Concord River below Meadow Brook (station no.1099500) is about 18.4 m 3 /s.
The SuAsCo watershed has a humid continental climate with snowy winters and warm summers.The annual average precipitation in SuAsCo is 1210 mm, with annual temperature and evapotranspiration at 9.15 ˚C and 650 mm, respectively, for the period 1973-2008.The forest was the predominant land use in the watershed in 2005.About 43% of the watershed was forested, 5% was under agriculture, pasture, and brushland, and about 35% was urban land based on land-use data from 2005.Wetlands constituted about 13% of the watershed.The operation of reservoirs in SuAsCo for recreation may affect streamflow, even though there are no withdrawals from these reservoirs.The Aquatic Life Use, an indicator for suitable habitat for flora and fauna, including water quality, was assessed to be impaired for 58% of the total acreage of the SuAsCo watershed [60].Impairment of Aquatic Life Use had resulted from many stressors.Those stressors include excessive total phosphorus and algal growth, low dissolved oxygen/saturation, flow regime alterations, and baseflow depletion from groundwater withdrawals and on-site septic systems.In addition, the slow-moving water of the rivers in this watershed and extensive wetlands depress dissolved oxygen levels and elevate phosphorus levels that contribute to eutrophication.

Water balance.
In HSPF, the watershed is subdivided into subbasins representing areas with different elevations and slopes.Each subbasin is divided into hydrologic response units (HRUs) based on different soil and land use combinations.The HRUs are further divided into impervious and pervious areas [61].Pervious landforms also represent wetlands.Each pervious HRU considers the following physical processes: interception, evapotranspiration, surface detention, infiltration, surface runoff, shallow subsurface flow (interflow), deep percolation, and baseflow.However, surface detention and surface runoff are the only components simulated in impervious areas.
Water Budget Impervious (IWATER) module was used to simulate the retention, routing, and evaporation of water from an impervious land segment.Unevaporated water produces surface runoff from impervious areas.In watershed segmentation, the area of impervious land segments represented by IMPLND is called the "effective" impervious area, or EIA, rather than the total or mapped impervious area [62].An EIA is any impervious area directly connected to a drainage system, such as streams, rivers, lakes, storm drains or rooftops that drain directly into storm drains.Runoff that drains into pervious areas before entering impervious areas was not included in the impervious simulation.EIAs are usually lower than impervious areas in watersheds, particularly in less densely populated areas.However, in highly urbanized areas, the EIA and total impervious areas may be very similar.
Slope (SLUR) for the overland flow plane, Chezy-Manning equation (NSUR), and average values of the surface roughness length (LSUR) of each HRU are used to generate overland flow [63].These capitalized abbreviations are for coefficients that undergo parameter calibration.Stream channels, lakes, and reservoirs are represented by reaches (RCHRESs), and these reaches receive flows from pervious and impervious areas via channel networks.A conceptual framework for hydrological processes in HSPF is shown in Fig 2 .An essential component of water balance is actual evapotranspiration (ET).Based on land surface vegetation, interception is simulated, assuming an interception storage capacity for rainfall abstractions.Water is lost initially in the form of evapotranspiration from the riparian vegetation (wetlands), then interception storage, upper-zone storage, active groundwater (AGWETP), lower-zone storage (LZETP), and baseflow (BASETP) [64].
Three conceptual parameters are used in HSPF to separate precipitation and snowmelt into two fractions that undergo 1) infiltration or 2) runoff.Those three conceptual parameters include a surface storage capacity value (UZSN), an infiltration-capacity index (INFILT), and an interflow-inflow index (INTFW) [62].The fraction of infiltrated water entering groundwater storage is considered active groundwater outflow and is available for discharge to surface channels.The fraction of infiltrated water that enters deep aquifers is known as inactive groundwater and is controlled by DEEPFR.Three parameters are used to calculate the outflow of groundwater from active groundwater storage: 1) active groundwater storage (AGWS); 2) active groundwater recession coefficient (AGWRC), 3) and an active groundwater outflow modifier (KVARY) [65].A list of parameters used for streamflow calibration is given in Table 3.

Sediments (TSS) and Nutrients (TN, TP)
. HSPF models inputs of sediments such as sand, silt, and clay in the form of loads from neighboring lands and upstream reaches [24].The model considers the exchange of sediments between the suspension and bed and the outflows of sediments to the downstream reaches.The erosion rate of sediments is based on the Universal Soil Loss Equation (USLE) [66].The critical shear stress theory was used to model scour, deposition, and transport of sand silt and clay particles [63,67].The critical bed shear stress parameter (TAUCS) was used for deposition.
Fig 3 demonstrates the schematic representation of fluxes and storages used in sediment simulations.During SEDMNT, two sediment fluxes are directly added to the detached sediment storage variable DETS, while the other fluxes are calculated in subordinate routines [62].The net additions or removals of sediment caused by human activity or wind are computed using NVSI.Sediment removal by water is simulated by scouring the matrix soil (SCRSD) and washing of detached sediment (WSSD).During the washoff process, sediment is detached from or attached to the soil matrix and transported.Rainfall results in detachment (DET).When there is no rainfall, attachment occurs; AFFIX specifies the attachment rate.Overland flow is the method of transport for detached sediment.Scouring involves both picking up and transporting the matrix soil by overland flow.
A list of parameters used for sediment calibration is given in Table 3.The following steps were taken for sediment calibration.1) Target sediment loading rates were estimated from the landscape as a function of land use, topography, and management practices.2) The parameters of scouring (KGER, JGER), deposition (TAUCD, TAUCS), and transport parameters (Table 3) were used for the stream channel to simulate the physical dynamics of the streams/ waterbodies.3) Overall sediment budgets for stream contributions and land use were evaluated.4) Observed and simulated sediment concentrations, including particle size distribution and load information, were compared.Land use is one factor affecting sediment loading, so a range of parameters was used to represent sediment loading based on land use (Table 3).For example, the AFFIX parameter measures the soil's tendency to compact.Wetlands and forest soils resist compaction more effectively because more organic matter is present in those soils.Therefore, forests and wetlands soil have lower values for the AFFIX parameter's typical range (0.01-0.1).
Phosphorus loading from the landscape is usually particle related.The sediment potency factor (pounds of phosphorus per ton of sediment) was used for phosphorus simulation [68].Hence sediment calibration has a significant influence on phosphorus calibration.Phosphorus was simulated in interflow and groundwater flow.Nitrogen loading from the landscape was simulated as ammonia and nitrate-nitrite.Nitrate-nitrite was represented as a buildup-wash off parameter on the land surface and is associated with groundwater and interflow.Ammonia was a minor constituent of the total nitrogen loading for most land uses and was modeled in the same fashion as nitrate-nitrite.Nitrate-nitrite loading comes from groundwater and interflow sources.Unlike phosphorus, nitrate-nitrite was not simulated as a sediment-associated (sediment-dependent) variable.

Land Transformation Model (LTM)
In regression based LULC change models, variables that predict spatial location are correlated with locations of predicted change [69].The methodologies include logistic regression [70], hedonic price models [71], and artificial neural network models [50].GIS has been used a lot more recently to house LULC change models with spatial analytical applications.Models like LTM that predict the spatial and temporal patterns of land use provide more information about the impacts of change than do approaches that focus on estimating aggregate land use amounts within geographic units, such as counties [72].
GIS-based LTM models are becoming more popular as reliable ANN-based models to forecast land-use changes.Using spatially explicit digital maps, these models simulate transitions and produce graphical outputs.A GIS preprocessing process helps manage and control historical spatial data layers, while ANNs help to identify input patterns (historical land use dynamics) and socioeconomic, political, and environmental drivers.Using ArcGIS grid files, model outputs can be directly incorporated into a hydrologic model [73].Therefore, using LTM, changes in land use can be linked to ecohydrological models for groundwater flow, solute transport, and forest cover [74].Studies have demonstrated that the LTM can accurately predict land use data from the national land cover data (NLCD) [75][76][77].
The error level of the LTM model used in this study was stabilized over 250,000 cycles, and the model output validity and reliability were assessed using kappa coefficients and percent correct metrics (PCMs).While reduced goodness of fit for our models could originate from challenges associated with the lack of a standard land-use/cover classification system across all datasets, one of the reported strengths of a neural network is its ability to perform well on datasets that are not part of the training dataset [51,78].Model outputs in LTM are in the form of binary results of "change (1) and "no change (0)" in land use.A type of ANN, a multi-layer perceptron (MLP), is used to understand the interaction between these factors or input drivers for changes.This information informs the LTM model to predict locations of change or no change with the help of the software Stuttgart Neural Network Simulator (SNNS).The process of LTM modeling involves the following steps.1) The first step of processing spatial data involves processing inputs to the LTM model in a series of base layers (Fig 5) [72].2) In step 2, the spatial effects of driver cells on land-use transitions are quantified with the help of GIS [50].3) In the next step, predictor variables are integrated with the help of MLP and SNNS.Finally, the number of cells that transitioned to urban land is calculated using the principle index driver (PID).
U t is the amount of newly developed urban land required in the time interval t, and A for the given time interval t is the per capita requirements for urban land.dP/dt in a given time interval in a particular area is the number of newly arrived people.The details about LTM GIS parametrization and ANN parametrization is described in [52].Finally, a percent correct match (PCM) metric is used to calculate the model's accuracy.

PCM ¼ # cells correctly predicted to change cells transitioning * 100 ð2Þ
An additional Kappa statistic calculates the percentage of LTM's success relative to chance.To identify the critical factors for land-use changes and predict the present, future, and past land-use change, the LTM model, has been applied and validated in various locations worldwide [79][80][81].In this study, the already calibrated LTM model [82] was used for the SuAsCo watershed to predict the present, future, and past land-use change and its impacts on runoff, sediments, total nitrogen, and phosphorus.
Three LTM scenarios: 2035, 2065, and 2100 were used to see the impacts of LULC on hydrology, sediments, and nutrients.Since the objectives of this study include studying the impacts of land-use change on water quantity and quality, meteorological input data was not changed for a future scenario.Instead, meteorological data from the baseline period  was used for future land-use scenarios to determine how water quality and quantity will be affected by changing land use.

Model inputs and calibration & validation data set
Meteorological input data were obtained from three weather stations in the SuAsCo watershed: Worcester WSO AP (Station no.MA 199923), Walpole 2 (Station no.198757), and Bedford (Station no.MA 190535) for time 1973-2008 (Table 1).Data from meteorological stations include hourly precipitation, air temperature, and solar radiation (Table 2).Missing data was filled with next available value.In this study, single-site calibration at the outlet is considered appropriate for a semi-distributed model.The performance of multisite calibrations in a small watershed with similar land-use types and meteorological conditions (uniform rainfall throughout the watershed) is almost identical regardless of how many parameters are calibrated [83].However, a multisite calibration scheme can increase parameter complexity and uncertainty in streamflow projections along with parameter equifinality.Studies in catchments of different sizes under various climates have shown that temporal parameter transfer outperforms spatial and spatiotemporal transfer modes for semi-distributed models [84][85][86][87].
Stream gage located at Concord River below Meadow Brook (USGS #01099500, RCHRES 157) was used for HSPF model calibration for 36 years from January 1, 1973, to December 31, 2008.Also, one stream gage at Sudbury River (USGS #01098530, RCHRES 140) and two stream gauges at Assabet River (USGS #01097300, RCHRES 99, USGS #01097000, RCHRES 142) were used for hydrology validation.In addition, observed data for sediments and nutrients were obtained from Massachusetts Department of Environmental Protection [88].Sediment and nutrient calibration were done after the hydrology was deemed acceptable because water quality modeling is dependent on the streamflow and hydrology of the modeled system.

HSPF multi-objective calibration.
In this study, temporal calibration is performed for time 1973-2008 to assure annual and seasonal variability, and spatial validation helped ensure local water balance at the sub-watershed level.While internal and outlet streamflow runoff can be correlated, streamflow runoff is not the only criterion used for calibration/validation.The multi-objective function of HSPF is used during the automatic calibration of hydrological processes utilizing an expert system to calibrate HSPF (HSPEXP+) that includes several criteria describing different fit features between model outputs and observed data.This is because a single objective function is insufficient to properly evaluate the simulation of all the essential physical processes of a hydrologic system [46,[89][90][91].Criteria taken into account in this study include squared errors of daily flows, storm/runoff volumes, storm peaks, seasonal flows, monthly flows, 50% lowest flow exceedance, and 10% highest flow exceedance.In addition, interflow, ET, and base-flow recession rates were also considered as criteria for calibration.Three internal stream gauges were used for model validation.Model performance is also tested for other internal variables, such as sediments and nutrients throughout the watershed, which indicates the credibility of the model prediction.
HSPF as a continuous semi-distributed model requires a spin-up time (about two to three years) to reach a dynamic agreement and steady state with the initial water budget conditions of the watershed.In this research, simulations started on July 1, 1970, and starting the model with a wet season helped the optimization algorithm prevent response bias for initial conditions.These initial conditions affect direct runoff or overland flow, soil water storages, interflow, and baseflow.Then calibration is performed for the period 1/1/1973-12/31/2008.Seventy-two storms were selected at each gaging station.Storm events were selected based on the recommendations in BASINS Technical Note 5 [92].The starting date for each storm was the day precipitation began, and the ending date was the day when observed flow returned to prestorm (or almost prestorm) conditions.Based on the recommendations in BASINS Technical Note 5, storm events were selected.The day precipitation began was the starting date for each storm, and the end date was when observed flow returned to prestorm conditions.The values of parameters that drive the respective hydrological processes were adjusted to attain agreement between observed and simulated results.

HSPF model parameters.
A list of parameters used for streamflow calibration is given in Table 3.In addition, 16 parameters for sediments and five parameters for nutrients were used for adjustment during the calibration period, and optimized values for those parameters are shown in Table 3.Values for initial parameters were obtained from studies in watersheds with similar meteorological conditions [61], topography [32], land use, soil properties, and management practices [93].Based on the literature on HSPF application in a watershed, a total of 17 parameters were used for runoff (Table 3).Monthly values of hydrology-related parameters LZETP, CEPSC, and UZSN were used, and those values varied depending on land use.The values of other parameters, such as DEEPFR, INTFW, and IRC, were not changed based on land use and remained constant throughout the year.The range of calibrated model parameters for watersheds with similar land use, such as Ipswich, Blackstone, and Taunton River [32,61,93] watershed was used as guidance during the calibration process.
With the help of the HSPEXP+ program, calibration was performed.Based on whether the errors are within the limits of the criteria (10-20% error) (Table 4), HSPEXP+ follows a hierarchical order for analysis of statistics and advice generation for changing related parameters [23].In addition, optimum model parameter values (Table 3) were derived using the calibration process with the HSPEXP+ program [94] for daily, monthly, and yearly flows (Table 5).HSPEXP+ interactive interface allows users to edit the User Control Input (UCI) file of the HSPF model, run a simulation, plot HSPF hydrology outputs against observed values, compute error statistics on the simulation, and provide expert advice regarding how to improve calibration parameters.

Statistical test for model calibration and validation.
The statistical tests of model results were performed to compare simulated flow, sediment, TN, and TP loads with the observations.Those statistical tests include (1) regression coefficient: R 2 (2) percent stream flow difference [calculated as (total model simulated flow-total observed flow)/total observed flow], (3) the Nash-Sutcliffe efficiency (NSE) (Table 5) and Pearson Correlation Coefficient (CC (Q).These statistics are used to evaluate the predictive power of the calibrated model.
The statistical tests of model results were performed to compare simulated flow, sediment, TN, and TP loads with the observed (field measurements) flow, sediment, TN, and TP loads.Those statistical tests are (1) regression coefficient: R 2 (2) percent stream flow difference [calculated as (total model simulated flow-total observed flow)/total observed flow], (3) the Nash-Sutcliffe efficiency (NSE) (Table 5) and Pearson Correlation Coefficient (CC (Q).These statistics are used to evaluate the predictive power of the calibrated model.
where Y sim and Y obs are the simulated and observed streamflow at time step i, and Y avg the average observed streamflow for the simulation period, n is the number of time steps.3 Results

Calibration and validation of the HSPF model
The model is calibrated for the SuAsCo watershed for 36 years , and a summary of calibration results is given in Table 5. HSPEXP model performance criteria were satisfied in terms of errors percent in total flow (m 3 /s), storm or runoff volumes (mm), storm peaks (mm), storm peaks (mm), and seasonal flow (mm).The observed streamflow was generally in reasonable agreement with the simulated flow for calibration and validation with Pearson correlation � > 0.8 (Fig 6), R 2 � 0.7 or above (Table 5).The observed data were split for calibration and validation.Based on evaluation metrics in Tables 4 and 5, calibrated parameters performed almost equally well for the validation period, and no overfitting occurred.The model under-simulated the streamflow volume by 8.9%.On the other hand, low flows (Table 4) were overstimulated.
Although observed water quality data is not available at every modeled reach or for every land use category, scatter plots (Fig 7) help the modeler quickly visualize sediment and nutrient behavior in the entire watershed model.For example, a scatter plot for simulated TSS, TN, and TP mean daily TSS concentration (mg/l) (Fig 7 ) shows that simulated data is in reasonable agreement with observed data (Table 6) with R 2 > 0.65 and a Pearson correlation coefficient of about 0.8.

LTM scenarios and impacts on hydrology
Based on land use data from 2005, about 43% of the watershed was forested, and 35% was urban land.In addition, wetlands constituted about 13% of the watershed, and about 5% of the area was under pasture, brushland, and agriculture.Land-use changes from 2005 to the LTMprojected conditions for 2035, 2065, and 2100 are illustrated in Fig 8 for the simplified landuse categories used in the hydrologic model.Based on LTM scenarios, most land-use conversions were from forest to low-density residential development in SuAsCo.Agriculture/pasture is projected to decrease by 30% in LTM 2100 (Table 7).The commercial/industrial and highdensity areas are expected to increase by 72% and 62%, respectively.Medium-density and lowdensity residential areas will increase by 83% and 93%, respectively.Forested areas and wetlands decreased by half (50% decrease) and 45%, respectively, while open water will remain unchanged by 2100 (Fig 8).LTM model scenarios showed a 75% increase in effective impervious areas from 2005 to 2100 because of urbanization, which caused a 50% reduction in forested land for the 2100 scenario (Table 7, Fig 8).Consequently, a 2.7% increase in total annual runoff (streamflow), a 69.4% increase in yearly surface runoff, and a 6.4% reduction in evapotranspiration led to a 4.2% decline in interflow compared to the baseline scenario of 2005 (Fig 9).In addition, the range of annual average surface runoff from subbasins in the watershed was 0.44-17.46inches      the 2100 scenario (Table 8).Thus land-use change has the effect of rising flood peaks during storm periods and decreasing baseflows between storms.For example, Fig 11 shows the impact of land-use change scenarios on four individual storm events (one storm in each decade of the simulation period) where peak storm volume was highest and low flows (recession limb) were smallest under the 2100 scenario.Without changing metrological data for future land use scenarios, the annual average increase in winter (months 1-3) streamflow volume was 5.1% under the 2100 scenario (Table 8).But the increase in summer (months 7-9) streamflow volume (about 10%) is more significant than winter volume.There was a substantial increase of 90.6% in summer storm volume and a reduction in winter storm volume (about 22.6%) under the 2100 scenario.Precipitation is distributed relatively evenly over most of the year, at about 89 mm per month for most months, increasing to 101 mm in March and November.Nearly a 6% increase in winter storm 9 (Table 8) was predicted by 2100.Low flows (Q90 = flows exceeded 90 percent of the time) (Fig 12) were increased under the 2100 scenario (Table 8) because of a decrease in evapotranspiration caused by future deforestation.Additionally, a 22% increase is projected for an annual average of 10% low flows (inches) for 2100 scenarios (Table 8).A slight decrease of 1.1% in 10% high flows (Table 8) was also projected.Overall, the annual average storm (total discharge) peak volume (m 3 /s) increased by 4.1% (Table 8) for the 2100 land use scenario.Low flows were increased under the 2100 scenario because of a decrease in evapotranspiration caused by future deforestation.

Assessment of changes in land use on sediments (TSS) and nutrients (TP, TN) in SuAsCo watershed
The HSPF model was run under the LTM scenario to study the effects of land-use changes on TSS.The annual average TSS yield at watershed outlets increased by 21% in 2035, 35% by 2065 (or about mid-century), and 54% by 2100 (Fig 13).The sediment yield at the outlet was

HSPF calibration & validation and impacts of LTM scenarios
The HSPF-based hydrological model was developed for the SuAsCo watershed to study the impacts of LULC change on hydrology and water quality.The observed streamflow was generally agreed with the simulated flow for calibration and validation.However, precipitation variability over the basin can be the reason for some of the model errors.Additionally, the Sudbury River at Saxonville has an upstream pond that undergoes occasional regulation and is not simulated in the model, resulting in a low R 2 .Also, the stage-discharge relation is distorted by beaver dams resulting in the backwater that may have caused discrepancies between observed and simulated streamflow for the calibration period (Fig 6).
Overestimation of low flows could be because when wetlands are simulated as pervious lands (PERLNDs) in HSPF, ET water loss equals precipitation water in wetlands minus outflows from wetlands.Since wetlands receive surface and subsurface lateral flows from surrounding upgradient areas, water fluxes in wetlands are not limited to input by precipitation, ET loss, and outflows from wetlands.Outflows are dominated by groundwater seepage and evapotranspiration (ET) in wetlands present in Massachusetts [95].Hence, in reality, wetlands have more available moisture for ET loss than moisture from precipitation simulated in HSPF.

Assessment of changes in land use on hydrology in SuAsCo watershed
Based on LTM scenarios, an increase in flood peaks during storm periods and a decrease in baseflows between storms is caused by land-use change.A double-peaked storm hydrograph (Fig 11b) could result from high rainfall and high initial moisture on steep slopes over the impervious surface.Rapid surface runoff from adjacent areas to the stream channel and impervious areas such as roads and trackways contribute to the initial hydrograph peak.The lower bound of the runoff range is produced from the subbasin in forested land, and the upper bound of the runoff range is associated with subbasins with urban land use.Similar results were obtained in the literature [96,97].
Similar to the streamflow at the outlet of the SuAsCo watershed, a slight increase in total yearly streamflow has also been estimated in other studies with urbanizing watersheds in Germany and China [98,99].Under the 2100 scenarios, the increase in summer stream flow volume is partly because of higher potential evapotranspiration during summer in forested and agricultural areas compared to urban areas.Consequently, clearing vegetation and canopy cover can increase runoff.In addition, soils are drier in summer than in winter, so the rate of infiltration is higher during summer [100].As a result, nearly a 6% increase in winter storms  7).
Low flows were increased under the 2100 scenario because of a decrease in evapotranspiration caused by future deforestation.Rainfall with high intensities on land with dispersed and scattered vegetation (or a low canopy covering) can cause siltation.Also, the hydraulic conductivity of the soil can decrease due to disconnected macrospores triggered by intense rainfall [101,102].This increase in sedimentation, crusting, and surface soil compaction because of land cover change can lead to a reduction in infiltration that contributes to the rise in storm peak volume.While results show that land-use change impacts low flows more conspicuously than high flows, it could be due to hydrologic model characteristics.
In a semi-distributed model such as HSPF, the reliability of the HSPF simulation results is higher for low flow conditions due to its lower temporal variability.On the other hand, high flow magnitude shows changes hourly rather than daily, monthly, or annually.In other words, while the time scale of the HSPF model is sufficient to capture daily flows, it is not precise enough to capture the changes in peak flows (10% high flows) where daily average flows cannot reflect flashy storms that last for minutes.Although studies by Middelkoop [103] in Germany and Borah & Bera [104] in Midwest, USA, found that semi-distributed hydrologic models are efficient for long-term hydrologic simulation, they cannot capture intense single events.Beighley et al. [105] also reported that small peaks are more sensitive to increased impervious land than large peaks.Urbanization can increase annual runoff depths more than maximum yearly discharges.
The slight decrease of 1.1% in 10% high flows (Table 8) may also be because most of the high flows occurred during winter in SuAsCo, and land-use change reduced the winter stream storms.Despite that, a 1.2% increase in 50% high flow (mm) was predicted under the 2100 scenario.Verbunt et al. [106] also observed this slight increase in high flows.This marginal change in high flows at the SuAsCo watershed outlet indicates that the conversion of vegetation to urban areas impacts small peaks more significantly than higher peaks.Additionally, studies conducted in urbanized watersheds in Switzerland, Washington DC, USA, and California, USA Moglen & Beighle, [107] showed high flows and runoff from an upstream subbasin decline at the basin outlet.

Assessment of changes in land use on sediments (TSS) and nutrients (TP, TN) in SuAsCo watershed
With the decrease in forested area and wetlands by 50% and 45%, respectively, sediment yield is predicted to increase by 2100 in the SuASCo watershed.This result is because of decreased sediment deposition as wetlands are converted into developed areas under the future land-use scenario.In addition, wetlands would typically capture and retain high quantities of dislodged sediments during rainfall events [108], so a decrease in wetland area would increase the TSS loading downstream of the Concord basin.
The sediment load from the Concord basin was more than the Assabet and Sudbury basins because of projected future urbanization (Fig 8) in the Concord basin downstream of the SuAsCo watershed.Even in the baseline scenario of 2005, the Assabet basin is less urban, with scattered patches of high-density urban areas, than the Sudbury and Concord basins.On the other hand, high-density urban areas in the Sudbury and Concord basins are more clustered.In addition, urbanization can increase peak flow (Fig 8) which could produce significant amounts of sediment from channel enlargement.These results are consistent with a study done in the eastern part of Alabama, USA [109], where areas with elevated levels of TSS also have a high magnitude of surface runoff.Similarly, in nine urbanizing watersheds in the Louisiana area [110], the land-use change caused an increase in TSS concentration.
The Assabet River provides a more conducive environment for nitrogen fixation because of duckweed (Lemma), a host for nitrogen-fixing bacteria [88].Therefore, more tributaries were affected by nitrogen loads in the Assabet River than in the Sudbury River (Fig 13).The increase in TN yield is also due to an increase in urban areas and runoff volumes.The highest increase in TN load was observed in the downstream Concord basin because of more clustered urban areas than in the Assabet and Sudbury basins (Fig 8).A study in the urbanized tropical region of Brazil also showed increases in nutrient concentrations with increases in impervious land [111].Moreover, urban streams can differ significantly from forested ones, with low dissolved oxygen concentrations, high respiration rates, and dissolved inorganic nitrogen [112].The nitrogen load transported by these urban fast-flow pathways to streams predominantly comprises ammonium-nitrogen (NH 4 -N), particulates, and dissolved organic nitrogen rather than nitrate (NO 3 -N).
Lower water levels caused by dams/reservoirs and slow-moving water increase phosphorus loads in SuAsCo.The water in SuAsCo moves mainly along the surface of the land.However, water also flows along the pathway of near-surface lateral flow.So, phosphorus transport (PO 4 -P) into surface waters via sediments is favored during storm events.Field measurements in the Assabet, Concord, and Sudbury rivers showed that sediment oxygen demand is constant with time.Other than that, high levels of DO and supersaturation are caused by nutrients rather than sediments [88].The SuAsCo watershed encompasses many impoundments, and anaerobic conditions in deep water provide a conducive environment for phosphorus to unbind.This dissolved phosphorus is then taken up by aquatic plants [57], and these floating plants further reduce dissolved oxygen and create eutrophic conditions in the water.Changes in phosphorus loading over the LTM scenarios are more strongly related to increasing phosphorus outputs from urban rather than agricultural sources (Table 7).Tasdighi & Osmond [113] also found strong and significant relationships between urban land, TP, and TN pollution using regression models.
Sediment yield was predicted to increase by 2100 in the SuAsCo watershed because of the decrease in the forested area and wetlands by 50% and 45%.This is also because of decreased sediment deposition as wetlands are converted into developed areas under the future landuse scenario.In addition, wetlands would typically capture and retain high quantities of dislodged sediments during rainfall events [108], so a decrease in wetland area would increase the TSS loading downstream of the Concord basin.Sediments act as sinks and sources of nitrogen and phosphorus.Through physical and biochemical reactions, nutrients that enter water tend to transfer from water to sediments [114].In addition, nutrients accumulated in the sediments can be released into the overlying water through hydrologic disturbances that lead to eutrophication [115].Therefore, ignoring the inputs of nutrients released from the sediments and only controlling anthropogenic inputs of TN and TP cannot effectively restore eutrophic water.
A limitation of this study is that LTM did not address urban density differences or densification processes.The LTM does not consider intensity levels, which play an integral role in defining the quality of cities, in comparison to focusing only on transitions (from non-urban to urban areas) [53,116].Future studies can focus on incorporating such features into the LTM to simulate multiple urban densities.

Best management practices (BMPs)
BMPs can minimize the impacts of LULC change on hydrology and water quality.For example, adding rain gardens to neighborhood streets, installing depressed parking lot islands, and using permeable pavement in large parking lots are critical site-level planning and design strategies to achieve stormwater management and water quality goals.In addition, conservation practices such as filter strips and contour farming could be adopted in the critical sub-watersheds of the study area [117].Additionally, bioswales are site-specific BPMs, where grasses, shrubs, and wetland plants are planted with infiltration systems to retain and infiltrate stormwater.Infiltration trenches with engineered soil and gravel and the installation of bioswales provide additional stormwater storage and facilitate water infiltration into the surrounding soil and groundwater.Native vegetation can also be planted in vegetated swales to enhance their water quality benefits.In addition to structural BMPs, non-structural BMPs can reduce TSS, TP, and TN loading through practices like street sweeping and catch basin cleaning, as well as illegal discharge detection and elimination (IDDE).
The BMPs commonly implemented in agricultural areas include split fertilization, alternative tillage, and mulching [118,119].Sediment reduction can be achieved by stabilizing channel banks and grassed waterways [120].In minimizing sediment yields at the watershed level, structural BMPs are better than agricultural BMPs [121].While baseflow, percolation, and aquifer recharge can be increased, surface runoff can be reduced.In the SuAsCo watershed, rain gardens and low-impact development technologies were evaluated for their effectiveness in improving stormwater runoff phosphorus concentrations.It was shown that rain gardens could improve runoff water quality in the Nashoba Brook watershed (a tributary watershed for SuAsCo) when they are installed at different percentages of total residential land use [122].Further, a parking lot in Acton, Massachusetts, was designed to manage stormwater more efficiently by incorporating infiltration basins and rain gardens.The quality of water was also improved as a result.
Public education measures should be implemented to address potential pollution sources from residential and commercial activities, such as lawn care, septic system maintenance, car washing, and household chemical use and disposal.SuAsCo watershed structural BMPs include hydrodynamic separators and deep sump catch basins for pre-treatment of TSS and catch basins for reducing the direct discharge of highway runoff [123].In addition, redesigning the drainage ditch to create a vegetated retention area can also be used.Non-structural BMPs include reducing sand and salt application and frequently cleaning catch basins [124].To effectively mitigate runoff, nutrients, and sediment yields in the study area, the concerned decision-makers can adopt a suitable combination of BMPs.

Conclusion
LULC change caused by population and economic growth has adverse effects on rivers and streams.This study uses a hydrologic model HSPF with LTM to assess the LULC impacts on runoff, TSS, TN, and TP in the SuAsCo watershed.Overall, HSPF and LTM models reasonably simulate watershed processes and can be used to assess runoff volume, sediment, and nutrient yield at a watershed scale.This study has shown that storm runoff, and the quantity of sediments and nutrients entering freshwaters is increased by decreasing forests and expanding the urban area.An increase in effective impervious area and a reduction in forest area are projected to cause a 69% increase in surface runoff.In addition, the reduction in evapotranspiration due to reduced forested land also caused a substantial rise in summer and winter storm runoff under the 2100 scenario.In addition to changes in water balance, TSS yield at the watershed outlet increased by 54% in 2100.The rate of TP and TN loads in the watershed outlet are projected to increase by 13% and 12%, respectively, by 2100.The sediment and nutrient yields in the downstream region were generally higher in the downstream Concord River basin than in the upstream Assabet and Sudbury basins.
Land-use change should be given more attention on the watershed scale.Therefore, there is a need for sustainable and site-specific conservation measures to mitigate the disruptive consequences of LULC change dynamics on runoff, sediments, and nutrient yield.Land managers can mitigate the magnitude and frequency of floods on a regional scale through urban land use planning.This research showed that areas with high surface runoff should be targeted for BMPs, such as rain gardens, permeable pavement in large parking lots, and vegetated swales.Future research can focus on temporal and spatial variability in TSS, TP, and TN and uncertainty analysis for hydrological and LTM models.In addition, the effective combination of land use can reduce flood peaks, sediments, and nutrient loads at basin outlets or in different subbasins.Hence using BMPs to improve water quality is especially relevant in areas where cost-efficient solutions to water quality impairment are required.

Fig 1 .
Fig 1.The SuAsCo watershed: The watershed is shown along with three main rivers.The length of the rivers is also presented.https://doi.org/10.1371/journal.pwat.0000083.g001

Fig 2 .
Fig 2. Conceptual framework of LULC impacts on watershed system.The grey dotted rectangular box represents the watershed boundary.Framework of processes simulated by HSPF is shown.The EIA in future landuse change scenarios is also included.https://doi.org/10.1371/journal.pwat.0000083.g002 Fig 4 represents the conceptual framework of PQUAL (water quality) module of HSPF.A detailed description of equations and parameters necessary to simulate processes related to sediments and nutrients is given in BASINS Technical Note 8 (Sediment Parameter and Calibration Guidance for HSPF) [68].

( 11 .
18-443.5 mm) during 2005, and runoff increased for 2100 scenarios with a new range of 0.46-20.95inches (11.7-532.13mm) (Fig 10).The change of forested land to urban areas also increased the average storm peak volume.For example, the average storm peak volume increased from 57.9 m 3 /s in 2005 to 60.3 m 3 /s for

Fig 9 .Fig 10 .
Fig 9. Percentage change in annual average water balance components such as surface runoff, interflow, evapotranspiration, and total runoff under LTM scenarios is presented.For each water balance component, landuse scenarios of 2023, 2065 and 2100 are shown from left to right.https://doi.org/10.1371/journal.pwat.0000083.g009

Fig 13 .
Fig 13. a) Annual Average TSS Yield from the watershed outlet, b) Percentage change in annual average TSS Yield discharge from the watershed outlet with Future Land use Scenarios.https://doi.org/10.1371/journal.pwat.0000083.g013

Fig 15 .
Fig 15. a) Change in the annual average concentration of TN under land-use Scenario and b) percentage change in the annual average concentration of TN c) Change in concentration of TP and d) percentage change in the annual average concentration of TP under land-use scenarios.https://doi.org/10.1371/journal.pwat.0000083.g015

Table 1 . Three weather stations were used for climate data in SuAsCo.
Annual average value for precipitation, temperature and potential evapotranspiration is presented.
https://doi.org/10.1371/journal.pwat.0000083.t001Forvalidation,twogauging stations at Assabet River (Acton stream gage # 01097300 and Maynard stream gage # 01097000) and a stream gauge at Sudbury River (Saxonville stream gage #01098530) were used.Eight hundred-eight observations were used for sediment calibration (Fig 7).Out of 808 observations, 509 samples were collected from the Assabet River, 158 samples from the Sudbury River, and 141 samples were collected from the Concord River.Therefore, 919 observations were analyzed for total nitrogen calibration, including 617 samples collected from Assabet River, 205 samples from Sudbury River, and 97 samples from Concord River.

Table 2 . Model data inputs and sources and resolution is presented
. BASINS interface allows to import input data in HSPF based on watershed boundary. https://doi.org/10.1371/journal.pwat.0000083.t002

Table 3 . List of adjusted parameters for calibration of hydrology, sediments, and nutrients in the HSPF model.
Optimized values are shown for calibrated parameters.The range indicates that values may differ depending on soil and land use. https://doi.org/10.1371/journal.pwat.0000083.t003

Table 5 .
Calibration and validation statistics.R 2 , NSE and Pearson correlation (PC) is shown to evaluate model performance for daily, monthly, and yearly flows.

Table 4 . Summary of calibration statistics for SuAsCo watershed.
Percentage error (%) between observed and simulated flows and volume is shown.Criteria are also shown for the limits of the error.