Assessing the Difference between Soil and Water Assessment Tool (SWAT) Simulated Pre-Development and Observed Developed Loading Regimes

The purpose of this research was to assess the difference between Soil and Water Assessment Tool (SWAT) simulated pre-development and contemporary developed loading regimes in a mixed-land-use watershed of the central United States (US). Native land cover based on soil characteristics was used to simulate pre-development loading regimes using The Soil and Water Assessment Tool (SWAT). Loading targets were calculated for each major element of a pre-development loading regime. Simulated pre-development conditions were associated with increased retention and decreased export of sediment and nutrients when compared to observed developed conditions. Differences between simulated pre-development and observed developed maximum daily yields (loads per unit area) of suspended sediment (SS), total phosphorus (TP), and total inorganic nitrogen (TIN) ranged from 35.7 to 59.6 Mg km−2 (SS); 23.3 to 52.5 kg km−2 (TP); and, 113.2 to 200.8 kg km−2 (TIN), respectively. Average annual maximum daily load was less during simulated pre-development conditions when compared to observed developed conditions by ranges of 1307 to 6452 Mg day−1 (SS), 0.8 to 5.4 kg day−1 (TP), and 4.9 to 26.9 kg day−1 (TIN), respectively. Hydrologic modeling results indicated that the differences in annual maximum daily load were causally linked to land use and land cover influence on sediment and nutrient loading. The differences between SWAT simulated pre-development and observed contemporary loading regimes from this study point to a need for practical loading targets that support contemporary management and integrated flow and pollutant loading regimes.


Introduction
Rapid agricultural and urban expansion have caused wide-spread impairment to receiving waters [1][2][3][4]. Currently, over 54% of assessed stream and river miles are threatened or impaired for aquatic life harvesting (e.g., fishing and fish consumption) in the United States (US). [5]. Excessive sediment and nutrient inputs from agricultural and urban development are listed among the leading sources of impairment. The associated impairments due to excessive sediment and nutrient inputs are coupled to socioeconomic losses that exceed billions of dollars annually [6]. In addition to socioeconomic losses, studies have reported environmental problems that are not limited to the loss of physical habitat, water toxicity, alterations to trophic state, eutrophication and subsequent hypoxic conditions, reductions in biodiversity, and increases in invasive species [7]. Given the trajectory Department of Agriculture National Resource Conservation Service (USDA-NRCS) show native vegetation that is dependent on climate and soils properties. While the SWAT model can be used to simulate the daily timeseries flow, sediment, and nutrient data needed to quantify loading regimes, there is a need to assess SWAT simulated pre-development conditions of sediment and nutrient regimes that are important for defining a pre-development loading regime.
The primary objective of the current work was to assess the difference in loading regimes between SWAT simulated pre-development conditions and observed developed conditions to inform policy makers and managers of what it might take to reach pre-development conditions. Subobjectives were to quantify: (i) a sub-set of statistically significant indices to represent each key element of a pre-development loading regime; (ii) the extent of alterations to sediment and nutrient loading regimes; and, (iii) loading targets based on pre-development loading estimates using four years of observed suspended sediment, total inorganic nitrogen, and total phosphorus data collected at five nested gauging sites in a representative mixed-land-use watershed of the central US.

Study Site
Hinkson Creek Watershed (HCW) is a mixed-land-use (i.e., 32% forested, 39% agricultural, 26% urban, and 3% wetland/water) catchment, which was located in the Lower Missouri River Basin ( Figure 1). The catchment has a drainage area of approximately 228 km 2 and elevations that range from 274 m above mean sea level (AMSL) in the headwaters to 177 m AMSL at the outlet. The headwaters of HCW are primarily rural agricultural lands and forested areas; and, the lower reaches are mostly urbanized (Table 1). Approximately 60% of the city of Columbia, Missouri (population 113,225) is located in the lower elevations [35].    (26) 88.9 (39) 73.0 (32) 6.8 (3) Soils in HCW are comprised of poorly drained to well drained prairie-forest transitional soils [36]. Mexico-Leonard series soils in the headwaters are clay loam till soils that were underlain by Pennsylvanian sandstone. Weller-Bardley-Clinkenbeard series soils in the lower reaches are composed of cherty clay solution residuum underlain by Mississippian limestone [36,37]. Soils in the headwaters are associated with increased surface runoff potential (hydrologic groups C and D) with infiltration rates that were as low as 0.001 cm h −1 [38] due to the presence of a well-developed claypan of smectitic minerals found in the argillic (Bt) horizon. Floodplain alluvial soils in the lower reaches have infiltration rates ranging from 0.1 cm h −1 at agricultural sites to 126.0 cm h −1 at bottomland hardwood forest sites [39].
Climate in Missouri is dominated by maritime and continental tropical air masses in the summer and continental polar air masses in the winter (http://climate.missouri.edu/). A recent 18-year climate record (2000-2017) showed average air temperature and annual total precipitation of 13.5 • C and 967 mm, respectively, in HCW. Air temperatures ranged from −23.1 • C to 41.3 • C. A well-defined wet season occurs during spring (April to June), and a relatively sporadic wet season occurs during early fall (September to October) [40].
The main channel of HCW, Hinkson Creek, is approximately 56 km in length (including meandering) ( Table 1); flowing southwesterly from the headwaters in Hallsville, Missouri through Columbia, Missouri to its confluence with Perche Creek. After confluence with Perche Creek, Hinkson Creek flows approximately 18 km (including meandering) through Perche Creek to its confluence with Missouri River. During winter 2008, Hinkson Creek was instrumented with a nested-scale experimental watershed study design [41]. Five hydroclimate gauging sites were positioned along the main channel to divide HCW into five sub-basins. (Table 1, Figure 1). Site #5 is located near the watershed outlet, and Sites #1 to #4 are nested within (Figure 1). At Site #4, the stage of Hinkson Creek has been intermittently monitored since 1967 using a USGS gauging station (#06910230).

Data Collection
Precipitation and stage were monitored at each site (n = 5) during the study period (2010-2014) ( Figure 1). Rainfall was sensed using TE525WS tipping bucket rain gauges with an accuracy of +1% at 2.5 cm h −1 rainfall rate to −3.5% at 5.1 to 7.6 cm h −1 rainfall rate. Stage was sensed using Sutron Accubar ® constant flow bubblers with accuracy of 0.02% at 0-7.6 m to 0.05% at 7.6-15.4 m. Precipitation and stage data were stored on Campbell Scientific CR-1000 data loggers. Methods that were used to estimate discharge included standards that were published by Turnipseed and Saur (2010) [42]. Stream velocity was measured using FLO-MATE™ Marsh McBirney flow meters and wading rods when stage was less than 1-m deep. Storm flows were measured using a USGS Bridge Board™. Stage-discharge rating curves were generated using the incremental cross section method to estimate stream flow [43].

Pre-Development SWAT Modeling
The SWAT model was used to simulate hydrologic conditions of a loading regime during pre-development conditions (circa 1800) at each gauging site in HCW. The model was previously calibrated and validated at each gauging site concurrently from the headwaters to the watershed outlet using four years of observed daily streamflow, sediment, phosphorus, and nitrogen data in the study catchment [44]. The process included auto-calibration using SWAT-cup as per methods that were published by Arnold et al., 2012 [34]. Calibration parameters were carefully adjusted during auto-calibration to keep the parameters within physically realistic bounds [34,44]. This and additional information regarding the input datasets and model forcing's can be found in Zeiger and Hubbart (2016) [44].
Land use and land cover was converted from a calibrated SWAT model [44], as per recommendations that were published by Bressiani et al., 2015, who suggested using a calibrated SWAT model prior to land cover conversion to simulate pre-European development land use and stream management conditions [45]. The effects of land use and land cover changes on loading regimes were the focus of this research, therefore, current climate and elevation input data were used as inputs in the SWAT model [45]. National Land Cover Dataset (NLCD 2011) input data were converted to native land cover dependent on soil characteristics using NRCS-USDA ESD and OSD. For example, all the areas with 'Mexico' soils were assumed native grassland land cover pre-development, and areas with 'Weller' soils were assumed forested land cover pre-development. Management operations were changed to reflect native vegetation. For example, native grassland land cover was attributed 'tall fescue' and forested land cover was attributed deciduous tree cover. Soil Conservation Service curve numbers (CN2; a parameter that controls surface runoff) and overland roughness (OV_N; a parameter that controls surface roughness) parameters were updated to reflect native land cover under 'good' conditions using values that were provided in the USDA Technical Release 55 [46]. Urban Harvester soils where converted to native soils that are located in the sub-basin (e.g., Weller soils) to simulate pre-development conditions as Harvester soils are described in an USDA OSD as "disturbed materials in areas where loess deposits that have been graded and reshaped for urban and suburban development". All of the other model parameters were left at default excepting parameters adjusted during model calibration outlined in Zeiger and Hubbart (2016) [44]. In summary, the general work flow of pre-development SWAT modeling was completed in the following six steps: (1) SWAT was set-up using contemporary elevation, soils, land use and land cover (LULC), climate and management operations data in HCW [44], (2) SWAT was calibrated and validated at each sub-basin concurrently from the headwaters to the watershed outlet using four years of observed stream flow, sediment, nitrogen, and phosphorus data collected at five nested gauging sites in HCW [44], (3) LULC was updated to reflect native land cover prior to European settlement (circa 1800), the SWAT parameter database (SWAT2012.mdb) was updated to reflect native land cover prior to European settlement (circa 1800), (4) SWAT input files were re-written, (5) SWAT was run, and (6) model output data were exported for comparison of simulated pre-development conditions and observed developed conditions.

Statistical Analysis
A USGS R package EflowStats was downloaded (https://github.com/USGS-R/EflowStats) and was used to calculate 171 'ecologically relevant' hydrological indices [11] at each gauging site in HCW. A detailed description of the 171 indices that were included in EflowStats can be found in "Users' Manual for the Hydroecological Integrity Assessment Process Software" [33]. Ultimately, 171 hydrological indices were calculated from four-year daily sediment and nutrient loading records during baseline conditions (i.e., SWAT simulated pre-development sediment and nutrient loading regimes) and observed current developed conditions at each gauging site. Principal components analysis was performed to examine intercorrelation and to minimize multicollinearity among indices [11,13]. One primary index was selected for each major element of a loading regime, including magnitude, frequency, duration, timing, and rate of change across low, average and high flows, similar to methods used by Olden and Poff (2003) to define a natural flow regime [13]. A total of eight flow elements were considered in this study, including low, average, and high magnitude; low and high duration, frequency, timing, and rate of change. Olden and Poff (2003) divided flow magnitude indices (n = 94) into average, high, and low flow sub-categories, while flow duration (n = 44) and flow frequency indices (n = 14) were divided into high and low flow sub-categories. Timing (n = 10) and rate of change indices (n = 9) were not divided into sub-categories for a grand total of nine sub-categories. In the current work, the sub-categories were considered excepting low frequency (n = 3) and high flow frequency indices (n = 11) were lumped into one category that was termed 'frequency', similar to how Olden and Poff (2003) lumped average (n = 3), high (n = 4), and low flow (n = 3) timing variables into one category termed 'timing'.
The PCA analyses were conducted using an R software package FactoMineR that was published by Lê, Sébastien et al., 2008 [47]. Broken-stick models were used to select the number of statistically significant principal components where the number of observed eigenvalues that exceeded simulated random eigenvalues showed the appropriate number of principal component axes for use in final PCA [48]. Results from PCA, in part, showed extracted eigenvectors (i.e., loadings) to each hydrologic index for each principal component axis. Absolute 'loadings' were sorted in ascending order to isolate indices that explained the most variance along each principal component axes [11,13]. Indices were ranked by loadings. Ranks were summed across all loading regimes (i.e., sediment, TIN, and TP) to select a subset of eight indices by greatest overall rank. The subset of eight indices was used to calculate loading targets corresponding to each key component of a loading regime at each gauging site.

Pre-Development SWAT Modeling
Climate inputs included wet (precipitation = 1323 mm), average (precipitation = 984 mm), and dry years (precipitation = 668 mm), in compliance with model calibration guidelines that were suggested by Arnold et al., 2012 [34], who suggested that the model calibration period should include wet, average, and dry years ( Table 2). Air temperatures ranged from −28.7 • C during winter months to 42.8 • C during summer months with an annual mean air temperature of 14.1 • C, and total solar ranged from 0.5 to 29.5 MJ m −2 , with a mean of 13.5 MJ m −2 ( Table 2). When considering that current climate and topography inputs were used in pre-development SWAT model simulations, observed differences in discharge and loading regimes were assumed due to land use and land cover changes. Differences in mean daily discharge between observed developed and simulated pre-development ranged from −0.6 m 3 s −1 (rural Site #2) to 0.5 m 3 s −1 (urban Site #5), showing that discharge generally decreased at rural sites and increased at urban sites following land use and land cover changes (Table 3). Similarly, average annual total discharge generally decreased in the rural headwaters and increased in the lower urbanizing reaches. Differences in median daily discharge between observed developed and simulated pre-development ranged from −0.4 m 3 s −1 (rural Site #1) to −0.7 m 3 s −1 (urban Site #4), indicating decreased median discharge following development at each gauging site. Differences in maximum daily discharge between observed developed and simulated pre-development ranged from 12.1 m 3 s −1 (rural Site #1) to 166.0 m 3 s −1 (urban Site #5) showing increased maximum daily discharge at each gauging site. Variability in daily discharge increased by 0.9 m 3 s −1 (rural Site #2) to 6.9 m 3 s −1 (urban Site #5). These results from SWAT simulations correlated land use and land cover alterations to the magnitude of discharge in HCW. These results indicated that land use and land cover alterations have increased surface runoff and the subsequent streamflow response at rural and urban sites of the study catchment. While median discharge decreased at rural sites, maximum discharge and human impacts on sediment and nutrient inputs increased sediment and nutrient yields in HCW. However, identifying the specific hydrologic mechanisms that decreased base flow at rural sites and increased the base flow at urban sites was beyond the scope of the current work. Mean daily sediment, TIN, and TP loading generally increased by an order of magnitude as average annual total loading increased from one to three orders of magnitude following the development in HCW. For example, observed developed mean daily sediment loads ranged from 13.2 to 68.1 Mg day −1 , while simulated pre-development mean daily sediment loads ranged from 0.2 to 13.4 Mg day −1 . Observed developed average annual total sediment loads ranged from 4812.2 to 24,877.2 Mg day −1 , while simulated pre-development average annual total sediment loads ranged from 62.1 to 4907.9 Mg day −1 . Results showed that, on average, annual total nutrient loads increased from 842 to 25,823 kg day −1 (TIN), and 322 to 5045 kg day −1 (TP). Thus, while simulated average annual total discharge decreased in the headwaters, the results showed consistent increases in sediment and nutrient export due to land use and land cover changes in HCW.

Statistical Analysis
Broken-stick models indicated the number of statistically significant principal component axes to consider ranged from 1 to 4 (PC1-PC4) ( Table 4, Figure 2). The principal component axes were orthogonal, and thus, independent. The first principal component (PC1) spanned the direction of the most variation in loading indices, the second principal component (PC2) spanned the direction of the second most variation in loading indices, and so on. Results showed PC1 explained 38 to 81% of the variation in loading indices with eigenvalues that ranged from 4.2 to 28.3. The second principal component explained an additional 7 to 35% of the cumulative variance in loading indices with eigenvalues that ranged from 1.1 to 13.9. There was no clear difference of intercorrelation between sediment, TIN and TP loading indices.   [11].

ma1
Average daily load of the entire study period. ma4 Standard deviation of the percentiles divided by the mean of percentiles (every 5th percentile is included). ma9 Spread in daily loading is the ratio of the difference between the 90th and 10th percentiles to median. ma12-23 Monthly median load of each month (January = ma12; December = ma23). ma24-35 Monthly coefficient of variation for each month (January = ma24; December = ma35). ma45 The mean load minus the median of the annual means divided by the median of the annual means. ml1-12 Monthly median minimum load for each month (January = ml1; December = ml12). ml13 Coefficient of variation of minimum monthly loads. ml19 The ratios of the minimum annual load to mean annual load for each year. mh1-12 Monthly median maximum load for each month (January = mh1; December = mh12). mh13 Coefficient of variation of maximum monthly loads. mh21 The average load for events above median load. mh23 The average load for events above seven times the median. mh26 The average peak load value for events above seven times the median. dl2 Annual minimum of 3-day moving average load. dl5 Annual minimum of 90-day moving average load. dl6 Variability of annual minimum daily average load. dl7 Variability of annual minimum of 3-day moving average load. dl9 Variability of annual minimum of 30-day moving average load. dl12 Annual minimum of 7-day moving average load divided by the median. dh1 Annual maximum of 1-day moving average daily load. dh14 The 95th percentile value divided by the mean of the monthly means. dh18 Average duration of events with loads above three times the median. dh20 Average duration of events with loads above the 75th percentile. fl1 Average number of events with loads below the 25th percentile. fh1 Average number of events with loads above the 75th percentile. fh4 Average days per year that the loading is above seven times the median load. fh5 Average number of loading events with loads above the median load. ta1 Loading constancy. ta2 Loading predictability. ta3 Seasonal predictability of loading. th2 Variability in Julian date of annual maxima. th3 Seasonal predictability of non-loading. tl1 Julian date of annual minimum. tl3 Seasonal predictability of low loading. ra2 Variability in rise rate. ra3 Fall rate. ra4 Variability in fall rate. ra5 Number of day rises. ra6 Median positive change of loading.

Loading Targets
Results from PCA indicated that a total of 20 unique loading statistics explained the most variance among 171 loading indices of PC1 across all of the loading regimes (i.e., sediment, TIN, and TP) in the study catchment. Ranking indices by "loadings" of PC1 across all of the loading regimes reduced the 171 loading indices to a subset of eight indices appropriate for loading targets in the study catchment. The subset of indices with greatest rank included: (1) median load during January (ma12), (2) median maximum load during January (ml1), (3) the average load for events above seven times the median (mh23), (4) annual minimum of 90-day moving average load (dl15), (5) annual maximum of one-day moving average daily load (dh1), (6) average number of events with loads above the 75th percentile (fh1), (7) loading constancy (ta1), and (8) variability in rise rate (ra2).
Average and low magnitude loading indices (ma12 and ml1) were generally greater during observed developed conditions with the largest differences in TIN loading. For example, observed ma12 (TIN) ranged from 3.4 to 16.9 kg day −1 , while the simulated pre-development ma12 (TIN) ranged from 0.5 to 1.8 kg day −1 ( Table 6). Average daily maximum load (dh1) was consistently greater during the observed developed conditions, especially for sediment loading. Observed dh1 (sediment) ranged from 1320 to 6930 Mg day −1 , and simulated dh1 (sediment) ranged from 11 to 477 Mg day −1 . Loading frequency was generally lesser at urban sites (Sites #3, #4, and #5) during developed conditions when compared to pre-development conditions. For example, fh1 sediment loading values were 44.5 days year −1 during pre-development conditions and 14.5 days year −1 during developed conditions at Site #3. Rate of change was consistently greater during developed conditions when compared to baseline conditions indicating increased variability in loading rise rates at all of the sites due to land use and land cover influences. For example, ra2 ranged from 437 to 736% during developed conditions and 134 to 319% during pre-development conditions. These results reflect a land use influence to loading magnitude, frequency, duration, and rate, especially during high flows.

Pre-Development SWAT Modeling
The advantages and disadvantages of continuous time physically based watershed scale models, like SWAT, have been well documented [49][50][51][52], and thus, were not the focus of the current work. However, model accuracy of pre-development simulations is often unknown due to a lack of historical observed data, hence the impetus for this investigation. The methodology that is used in the current work of calibrating the SWAT model to the constituents of interest (in this case flow, sediment, nitrogen, and phosphorus) prior to land use and land cover conversion has been accepted as an appropriate method to simulate pre-European development land use and stream management conditions [45]. Thus, while it was not possible to quantify model accuracy due to a lack of historical data, calibrating the model prior to land cover conversion helped to ensure more accurate simulations of pre-development flow, sediment, and nutrient regimes.
Historical land cover reconstructions were important for understanding pre-settlement loading regimes in this research. Native land cover was estimated based on physiographic site characteristics (e.g., climate, geology, landscape position, and soils). Natural resource management agencies use such methods regularly to provide the best estimates of ecological reference site conditions that are important for management of natural resources. Native vegetation coverage of many US watersheds can be found in Ecological Site Descriptions (ESD) (https://esis.sc.egov.usda.gov/) and Official Soils Series Descriptions (OSD) (https://soilseries.sc.egov.usda.gov/) that were developed by United States Department of Agriculture National Resource Conservation Service (USDA-NRCS). Once land cover was converted and SWAT simulations were run.
Results showed agricultural and urban development increased maximum flows, which translated to increased maximum pollutant loading in HCW. Differences in maximum streamflow ranged from 12 m 3 s −1 (Site #1) to 166 m 3 s −1 (Site #5). Differences in maximum daily sediment loading were relatively greater, ranging from 36 kg km −2 (Site #2) to 60 kg km −2 (Site #1). Differences in maximum daily TP loading were also large ranging from 23 kg km −2 (Site #2) to 53 kg km −2 (Site #5). Differences in maximum daily TIN loading were the greatest, ranging from 113 kg km −2 (Site #2) to 201 kg km −2 (Site #5). Increased pollutant loading was expected considering approximately 50% of total watershed area that was native forested land cover was cleared for agricultural and urban development in the study catchment. In addition, approximately 12% of the land cleared was replace with impervious surfaces associated with increased volume and velocity of surface runoff, and subsequent increases in peak flow in HCW [53]. However, the results that are indicated differences in streamflow alone did not account for the observed pollutant loading in HCW. For example, flow and load duration curves show increased median streamflow did not translate to increased median pollutant loading ( Figure 2). Model simulations showed that increased pollutant export (especially TIN loading) was also due to differences in management operations (i.e., fertilizer applications, grazing, and lawn fertilizer applications).
The upper Midwest of the Missouri River Basin (MRB) has been targeted as a critical source area of sediment, N and P loading [54][55][56][57]. The greatest nitrogen fertilizer application rates of more than 2.5 t km −2 year −1 [54], with nine central states, including Missouri, which contribute more than 75% of the total N and P terrestrial loads to the Gulf of Mexico [56]. In the current work, the maximum daily loading of TIN was four orders of magnitude greater during developed conditions when compared to pre-development conditions. These results highlight substantial changes to N loading, in part, due to land use alterations to streamflow regimes, but primarily due to fertilizer applications in agricultural and urban areas of the study catchment. Thus, the restoration of pre-development stream flow regime may achieve limited results without source reductions of sediment and nutrients in HCW and similar mixed-land-use watersheds.
Results from this work are in agreement with Bernhardt et al., 2008 and others [58][59][60], in that a small fraction of N was exported during pre-development conditions while excessive nutrient enrichment was observed during developed conditions. Major anthropic sources of N included crop and lawn fertilization, soybean cultivation, and grazing. Ultimately, results from the current work in combination with other regional studies show anthropogenic sources of increased N production have substantially altered N loading regimes [58,61]. subsequent increases in peak flow in HCW [53]. However, the results that are indicated differences in streamflow alone did not account for the observed pollutant loading in HCW. For example, flow and load duration curves show increased median streamflow did not translate to increased median pollutant loading (Figure 2). Model simulations showed that increased pollutant export (especially TIN loading) was also due to differences in management operations (i.e., fertilizer applications, grazing, and lawn fertilizer applications).

Statistical Analysis
Similar assumptions and limitations that are associated with environmental flows assessment using PCA were also applicable to the current work. For example, indices dependent on drainage area were area-normalized for better site comparison [11], which reduced variability that is important in the selection of loading targets using PCA. Additionally, PCA results are a product of statistical tests that are not physically-based. Selecting loading indices based on statistical results can be problematic when selected targets are not casually linked to loading alteration. Therefore, results should be carefully examined and validated against process-based understanding to ensure that indices are physically meaningful.
When the deviation between simulated pre-development and observed developed loading indices were examined using bi-plots, and the results showed a clear separation of pre-development and developed loading indices along the first principle component axes (PC1) (Figure 3). The second principal component axes (PC2) generally represented the effects of spatial scale on indices between sites as the drainage area increased from Site #1 in the headwaters to Site #5 near the watershed outlet ( Figure 3). The loading indices that were ranked by absolute loading values for each major element of the loading regime and the results are also shown in bi-plots (Figure 3). Indices with the greatest absolute distance from the origin in the bi-plot explained the most variability in the group of 171 indices. Land use and land cover conversion (associated with PC1) had greater influence on the loading regime than the effects of scale (associated with PC2) in the study catchment with the exception of pre-development sediment simulations where sediment loads increased between Sites #4 and #5 (Figure 3).  Table 5. Loading indices that overlapped were slightly offset to improve readability. Total inorganic nitrogen is TIN, and total phosphorus is TP.
Results showed land use induced loading alterations were evident across every major element of the sediment and nutrient loading regimes in the study catchment, especially high load magnitude and duration (Figure 4). For example, annual maximum of one-day moving average daily load (dh1) was associated with differences between pre-development and developed conditions that ranged from 1307 to 6452 Mg day −1 (sediment), 5 to 27 kg day −1 (TIN), and 1 to 5 kg day −1 (TP). Results indicated the increase in dh1 was attributed to (1): land use and land cover alterations to the volume and velocity of surface runoff, and (2) management operations (e.g., tillage, lawn and crop fertilization, grazing, and soybean cultivation) that increased the export of sediment and nutrients. Previous research in the study catchment showed nearly all suspended sediment and nutrient loads were exported during high flow interval (e.g., 99% of suspended sediments, 92% of TIN, and 95% of TP loads) [40]. Results from the current work show that development has increased maximum annual daily load of sediment and nutrient export by 3 to 4 orders of magnitude in HCW.  Table 5. Loading indices that overlapped were slightly offset to improve readability. Total inorganic nitrogen is TIN, and total phosphorus is TP.

Loading Targets
Flow, sediment, and nutrient loading regimes shape aquatic communities [10,14,58], as long-term patterns of water, sediment, and nutrient availability facilitate the spatial arrangement of soils, microclimate, and hillslope hydraulic flow paths [58]. In the current work, pre-development conditions were associated with increased retention and decreased export of flow, sediment, and nutrients when compared to developed watersheds. These results, in combination with previous literature [58], highlight potential ecological implications of land use induced loading alteration observed in Hinkson Creek and point to the need for loading targets to support ecological health.
Results showed land use induced loading alterations were evident across every major element of the sediment and nutrient loading regimes in the study catchment, especially high load magnitude and duration (Figure 4). For example, annual maximum of one-day moving average daily load (dh1) was associated with differences between pre-development and developed conditions that ranged from 1307 to 6452 Mg day −1 (sediment), 5 to 27 kg day −1 (TIN), and 1 to 5 kg day −1 (TP). Results indicated the increase in dh1 was attributed to (1): land use and land cover alterations to the volume and velocity of surface runoff, and (2) management operations (e.g., tillage, lawn and crop fertilization, grazing, and soybean cultivation) that increased the export of sediment and nutrients. Previous research in the study catchment showed nearly all suspended sediment and nutrient loads were exported during high flow interval (e.g., 99% of suspended sediments, 92% of TIN, and 95% of TP loads) [40]. Results from the current work show that development has increased maximum annual daily load of sediment and nutrient export by 3 to 4 orders of magnitude in HCW.   Table 5. Total inorganic nitrogen is TIN, and total phosphorus is TP.
Previous studies have reported land use alterations to the frequency of flow and non-point source pollutant delivery to receiving waters, especially in urban areas [62,63]. However, high load frequency (fh1) showed reductions to fh1 at urban Sites #3, #4, and #5 near the watershed outlet. The index fh1 was defined as the average number of events with loads that were above a threshold equal to the 75th percentile. Differences in fh1 were attributed to differences in the 75th percentile of loading between pre-development and developed conditions. For example, the 75th percentile of sediment loading threshold was 18.7 kg day −1 during pre-development conditions when compared to 566.6 kg day −1 during developed conditions at Site #3. These results indicated that management efforts focused on restoring a pre-development sediment regime should increase fh1 of sediment delivery at Site #3. However, it is important to note that increasing fh1 would only be appropriate coupled to a 30% reduction to the 75th percentile of sediment loading. Additional research is needed to validate the efficacy of loading targets based on a pre-development loading regime in contemporary mixed-landuse watersheds that have been extensively altered by past and present land uses.
Restoration of a natural loading regime may not be practical in developed watersheds. There are best management practices (BMPs) that can increase the residence time and the retention of pollutant  Table 5. Total inorganic nitrogen is TIN, and total phosphorus is TP.
Previous studies have reported land use alterations to the frequency of flow and non-point source pollutant delivery to receiving waters, especially in urban areas [62,63]. However, high load frequency (fh1) showed reductions to fh1 at urban Sites #3, #4, and #5 near the watershed outlet. The index fh1 was defined as the average number of events with loads that were above a threshold equal to the 75th percentile. Differences in fh1 were attributed to differences in the 75th percentile of loading between pre-development and developed conditions. For example, the 75th percentile of sediment loading threshold was 18.7 kg day −1 during pre-development conditions when compared to 566.6 kg day −1 during developed conditions at Site #3. These results indicated that management efforts focused on restoring a pre-development sediment regime should increase fh1 of sediment delivery at Site #3. However, it is important to note that increasing fh1 would only be appropriate coupled to a 30% reduction to the 75th percentile of sediment loading. Additional research is needed to validate the efficacy of loading targets based on a pre-development loading regime in contemporary mixed-land-use watersheds that have been extensively altered by past and present land uses.
Restoration of a natural loading regime may not be practical in developed watersheds. There are best management practices (BMPs) that can increase the residence time and the retention of pollutant loads by diverting runoff to naturally occurring near-surface valley bottom storage (e.g., bottomlands, wetlands, bank soils, and alluvial fills) [62,64]. However, the megatons of N and P fertilizers applied during March in the Mississippi River Basin is not a natural magnitude and timing of nutrient inputs for undeveloped watersheds. Such nutrient enrichment has caused extensive and wide-spread impairment to receiving waters [54,56,57]. Thus, source contributions will also need to be managed effectively to achieve loading targets that are based on pre-development loading regimes, especially nutrients in the study catchment.

Study Implications and Future Work
The progression of management efforts from sustaining minimum flows and the reduction of point source pollutant loading, to the restoration of a natural flow regime and the reduction of non-point source pollutant loading has helped to improve water quality and ecological health in the US. A focus on reductions to point and non-point source loading have been justified while considering that the magnitude of total loading is usually the most altered aspect of loading regimes currently. However, alterations to the duration, frequency, timing, and rate of change in loading have also been reported [14,62]. Building on work from Poff et al., 1997, Wohl et al., 2015 discussed how flow management plans that do not consider sediment inputs might achieve limited success considering aquatic ecosystem health depends on adequate flow and sediment regimes [10,14]. Results from this research are in agreement with Wohl et al., 2015 [14], there is a need to broaden the environmental flow assessment and load reduction efforts to account for alterations to all major aspects of loading regimes.
Building on work from Wohl et al., 2015 [14], the results from this study show that nutrient regimes can be more extensively altered than flow and sediment regimes in the central US. About 7% of assessed stream miles have been deemed threatened or impaired with flow alteration, while over 20% of assessed miles are threatened or impaired with sediment and/or nutrient alterations in the US. [5], with extensive nutrient alterations being reported in the Mississippi River Basin [54,56,57]. Such wide-spread impairment drives the impetus for future work that is focused on development of regional flow and loading targets to restore pre-development flow and loading remiges. Results from this assessment show the same methodologies outlined in regional environmental flows assessment frameworks (e.g., ELOHA and HIP) could be expanded to include loading targets as being appropriate for a more holistic approach to regional stream and river management.

Conclusions
Human land use, land-change related impacts have caused wide-spread alterations to flow and pollutant loading regimes, and the current pace of impairment exceeds mitigation efforts globally. The current work assessed the efficacy of loading targets based on SWAT simulated pre-development loading regimes in a contemporary mixed-land-use catchment. This modeling application was possible due to the long-term monitoring of precipitation, stage, and water quality constituents at nested gauging sites in the study catchment.
Results showed land use impacts to each major element of sediment and nutrient loading regimes in HCW. Differences in maximum daily loading were 55 Mg km −2 (sediment), 201 kg km −2 (TIN), and 53 kg km −2 (TP) at urban Site #5 near the watershed outlet. Land use alterations to high load duration (dh1) were also apparent with differences that ranged from 1307 to 6452 Mg day −1 (sediment), 5 to 27 kg day −1 (TIN), and 1 to 5 kg · day −1 (TP) between pre-development and developed conditions. Variability in the rise rate of loading was two to three times greater during developed conditions. These results indicated that loading targets dependent on pre-development loading regimes may not be socioeconomically practical in watersheds already extensively altered by land use. Ultimately, there is a need for practical loading targets that support balanced flow, sediment, and nutrient regimes to guide regional stream and river management.