Nutrient Inputs to the Laurentian Great Lakes by Source and Watershed Estimated Using SPARROW Watershed Models

Abstract Nutrient input to the Laurentian Great Lakes continues to cause problems with eutrophication. To reduce the extent and severity of these problems, target nutrient loads were established and Total Maximum Daily Loads are being developed for many tributaries. Without detailed loading information it is difficult to determine if the targets are being met and how to prioritize rehabilitation efforts. To help address these issues, SPAtially Referenced Regressions On Watershed attributes (SPARROW) models were developed for estimating loads and sources of phosphorus (P) and nitrogen (N) from the United States (U.S.) portion of the Great Lakes, Upper Mississippi, Ohio, and Red River Basins. Results indicated that recent U.S. loadings to Lakes Michigan and Ontario are similar to those in the 1980s, whereas loadings to Lakes Superior, Huron, and Erie decreased. Highest loads were from tributaries with the largest watersheds, whereas highest yields were from areas with intense agriculture and large point sources of nutrients. Tributaries were ranked based on their relative loads and yields to each lake. Input from agricultural areas was a significant source of nutrients, contributing ∼33-44% of the P and ∼33-58% of the N, except for areas around Superior with little agriculture. Point sources were also significant, contributing ∼14-44% of the P and 13-34% of the N. Watersheds around Lake Erie contributed nutrients at the highest rate (similar to intensively farmed areas in the Midwest) because they have the largest nutrient inputs and highest delivery ratio.


S1. Description of the SPARROW model
SPARROW is a GIS-based watershed model that uses a mass-balance approach to estimate nutrient sources, transport, and transformation in terrestrial and aquatic ecosystems of watersheds under long-term steady-state conditions (Smith et al., 1997;Alexander et al., 2008).
SPARROW includes non-conservative transport, mass-balance constraints, and water flowpaths defined by topography, streams, and reservoirs, based on a stream-reach network with delineated reach catchments. The model-estimated flux leaving each reach (i) in SPARROW, * i F , is given (S1) The first summation term represents the flux ( * j F ) from all upstream confluent reaches ( ) i J that are delivered downstream to reach i. i δ is the fraction of upstream flux delivered to reach i; i δ generally equals 1 unless the upstream end of reach i is the location of a diversion. ( ) ⋅ T is the stream transport function representing attenuation processes acting on the flux as it travels along the reach pathway (instream loss). This function defines the fraction of the flux entering reach i at the upstream node that is delivered to the reach's downstream node. The factor is a function of measured stream and reservoir characteristics, denoted by the vectors S Z and R Z , with corresponding coefficient vectors S θ and R θ .
The second summation term represents the incremental flux (that introduced to the stream network in reach i). This term is composed of the flux originating from specific sources, indexed by S N n , , 1  = . Associated with each source is a source variable, denoted n S . Depending on the nature of the source, this variable could represent the mass of the source available for transport to streams, or it could be the area of a particular land use. The variable n α is a source-specific coefficient that converts source variable units to flux units. The function ( ) ⋅ n D represents the land-to-water delivery factor. For sources associated with the landscape, this function along with the source-specific coefficient determines the amount of a constituent delivered to streams. The land-to-water delivery factor is a source-specific function of a vector of delivery variables, denoted by D i Z , and an associated vector of coefficients D θ . The last term in the equation, function ( ) ⋅ ′ T , represents the fraction of flux originating in and delivered to downstream end of reach i. This term is similar in form to the stream delivery factor defined in the first summation term but is used to transport flux only from the midpoint of the reach to the outlet.
Parameter coefficients associated with the sources, land-to-water delivery factors, and instream-and reservoir-loss terms were statistically estimated using weighted nonlinear least squares regression, based on calibrations with long-term mean annual normalized loads at each monitored station throughout the study area. A more indepth description of the SPARROW model and its calibration is given in Schwarz et al. (2006).

S3. Summary of measured (with Fluxmaster) and estimated (with SPARROW) loads for each monitored site
In Tables S2_TP and S2_TN (attached to the end of the Supporting information), we provide a summary of measured (with Fluxmaster) and estimated (with SPARROW) loads for each monitoring site used in the calibration of the SPARROW models. In this table we also provide a summary of the location of each of the sites and the concentration data used to estimate the loads with Fluxmaster.
For each of the sites, there is a description of its location and the concentration data used in Fluxmaster. For each site, midmonthly mean and median concentrations for each site were calculated using a subset of the data. The subset included only samples collected closest to the middle of the month for sites at which more than one sample per month was collected. This was done to reduce possible bias associated with the frequency of sampling at different sites (such as weekly or storm sampling events compared to the more common monthly or less than monthly sampling).
The long-term mean annual nutrient loads for each monitored site were computed with the rating curve/regression procedure in the Fluxmaster computer program (Schwarz et al., 2006). This procedure combines water-quality data at a monitoring station with daily flow values to provide more accurate load estimates than can be obtained by using individual water-quality measurements alone. TP and TN loads were determined with log-linear water-quality regression models that related the logarithm of constituent concentration to the logarithm of daily flow, decimal time (to compensate for trends), and season of the year (expressed using trigonometric functions of the fraction of the year). Regression models were fit to data from each potential load site (sites with > 25 samples and corresponding long-term flow data, see Saad et al., 2011 [this issue] for a more complete description). 3479, and 4961. The point source load estimates used in the model calibration were based on estimates for 2002. If the 2002 data were missing and 1997 data were available, then data from 1997 were used. If both 2002 and 1997 data were missing and 1992 were available, then 1992 data were used. The distribution of point source inputs of P and N were quite variable, but largest inputs were concentrated around major urban areas ( Figure S2

Nutrients from confined and unconfined manure
Manure inputs are based on county-level estimates of animal wastes from Ruddy et al. (2006). Nutrients associated with livestock wastes reflect contributions from the excreted wastes of unconfined animals on farms, pastures, and rangelands and from the excreted wastes of confined animals, including those in concentrated animal feeding operations. Confined animal wastes include recoverable manure that may be applied to nearby farmlands as well as unrecoverable manure that is lost during the collection, storage, and treatment of the waste.
Manure inputs were derived from 2002 county livestock population data from the U.S. Census of Agriculture using species specific rates. The county data were then allocated to each SPARROW stream catchment by the fraction of the catchment's agricultural land and grasslands. Highest confined manure input rates for P occurred in southern Minnesota, Iowa, southwestern Wisconsin, and northern Indiana, whereas highest unconfined manure input rates for P occurred in eastern Iowa, western Kentucky, and Tennessee ( Figure S4). The distribution of confined manure input rates for N is similar to that for P, highest inputs occurred in southern Minnesota, Iowa, southwestern Wisconsin, and northern Indiana ( Figure S5). A.

Inputs from specific land-use/land-cover types
Land-use/land-cover related inputs (additional agricultural inputs from cultivated agricultural areas, urban inputs from urban and open areas, and natural inputs from forested areas) were based on the respective amount of area in each of these land types based on 2001 National Land Cover Data (USGS, 2000). The area consists primarily of forested areas in north and southeast, and agricultural areas in the west and central areas. Several major metropolitan areas are in the basin, such as Chicago, IL, Detroit MI, and Cleveland, OH ( Figure S6).

Areas of tile drainage
Tile drain information was derived from the 1997 National Resources Inventory dataset compiled by the National Resource Conservation Service. The distribution of tile drains is inversely related to the permeability of the soil. The highest percentage of land with tile drains is found in the central part of MRB3 (southern Minnesota, Iowa, southeastern Wisconsin, central Illinois, southern Michigan northern Indiana, and northern Ohio), where the soils are least permeable ( Figure S8). Figure S8. Percentage of the area drained by tiles for MRB3, by SPARROW catchment. Intervals represent five equal quantiles (0-20%, 20-40%, 40-60%, 60-80%, 80-100%).

Drainage density
The drainage density of each catchment was calculated as the total length of all stream reaches in a catchment divided by its area. Stream reaches were obtained from the enhanced stream-reach file 1 (RF1; 1:500,000 scale) (Brakebill et al., 2011;this issue). The highest drainage density is in the southeast part of MRB3, which was partly caused by the increased precipitation in this area and the increased resolution in the enhanced RF1 coverage. Lowest drainage densities are in the northwest part of MRB3. Figure S9. Drainage density (summed length of all stream reaches in each catchment divided by its area) for MRB3, by SPARROW catchment. Intervals represent five equal quantiles (0-20%, 20-40%, 40-60%, 60-80%, 80-100%).

Air temperature
Air temperature data, averaged from minimum and maximum daily temperature values over 30 years (1971 to 2000), were obtained from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) digital data network (http://www.prism.oregonstate.edu/).
These data show a strong north to south gradient of increasing mean air temperature ( Figure   S11). Mean air temperatures range from 2.3 to 14.9 C.   Hydraulic loading was used to describe nutrient removal in reservoirs (Schwarz et al. 2006). Hydraulic loading for each reservoir was calculated as average flow divided by reservoir surface area based on information from the National Inventory of Dams (USACE, 2007). The reservoirs used in the SPARROW models are shown in Figure S15. Nutrient losses in reservoirs that are not included in this coverage were assumed to occur as part of instream loss.

SPARROW model calibration.
In the model calibration process, the studentized residuals (difference in the logarithmically transformed observed and simulated loads) were evaluated for evidence of local or regional biases based on visual inspections of residual maps ( Figure S16). In the calibration process, the model used the monitored transport from each upstream basin rather than the modelestimated transport to represent the load leaving the upstream basin in an attempt to eliminate prediction errors from cascading downstream to other monitored sites. In Figure 2 of the main manuscript, the full model predictability is evaluated, because errors are allowed to cascade downstream to other monitored sites. In Figure 2 of the main manuscript, full prediction errors are reported; however, in Figure S16, the residuals in the calibration process are using studentized residuals are reported to enable comparisons in the calibration of other MRB SPARROW models. Overall, the models showed little signs of regional biases. A. B.

Output from the MRB3 SPARROW TP and TN models for each Great
Lake and nearby large river basin (Tables S3 and S4), each HUC8 (Table S5) and each tributary with a drainage basin greater than 150 km 2 (Table S6).
The total TP and TN load and yield to each Great Lake and nearby large river basin is provided in Tables S3 and S4. The percent contribution of each source of nutrients is given in each 10.8 a Delivery ratio is computed as the delivered load divided by the total non-decayed load. Table S3. Total annual phosphorus loads and yields for each Great Lake and major river basin in Major River Basin 3, with the percent contribution by source.

Percent Contribution by Source
All loads and yields from individual tributaries were adjusted to remove known spatial biases by only predicting areas that are not monitored.
[kg, kilogram; km 2 , square kilometer; Yield is load per unit area of the watershed]  Robertson and Saad, 2011, Journal of the American Water Resources Association, Nutrient Inputs to the Laurentian Great Lakes By Source and River Basin Estimated Using SPARROW Watershed Models. Table S4. Total annual nitrogen loads and yields for each Great Lake and major river basin in Major River Basin number 3 (MRB#3), with the percent contribution by source.

Percent Contribution by Source
All loads and yields from individual tributaries were adjusted to remove known spatial biases by only predicting areas that are not monitored.
The total TP and TN load and yield from each 8-digit Hydrologic Unit Code (HUC8; Seaber et al., 1987) watershed throughout MRB3 and is provided in Table S5. The incremental nutrient load (load generated from within a HUC8 and delivered to its most downstream reach) was computed for each HUC8 watershed by summing the loads from each reach within the HUC8 watershed. The incremental yield was then computed by dividing the incremental load by the total area of the HUC8 watershed.
The lower and upper 95% confidence limits for the HUC8 loads and yields are also provided in Tables S5_TP and S5_TN (attached to the end of the Supporting information). The actual delivered fluxes computed in SPARROW are assumed to depend on a multiplicative error term that represents other sources and processes not included in the SPARROW analysis.
Because of this residual term, and because the determination of the predicted flux depends on coefficients that are estimated via statistical methods, the delivered yields across HUC8 watersheds are subject to uncertainty. Because of the nonlinear manner in which the estimated coefficients enter the model, it was necessary to use bootstrap methods (Schwarz et al., 2006) to assess the uncertainty. Bootstrap analyses were also used to correct for potential bias caused by log retransformations in the yield predictions. A brief summary of the bootstrap method is presented here; a full description of the bootstrap methodology is described by Robertson et al. (2009).
The bootstrapping method was implemented by performing 200 repeated calibrations of the SPARROW model using randomly selected integer weights (which sum to the total number of monitored reaches in both the models) applied to each of the squared residuals at monitored reaches, resulting in 200 realizations of the estimated coefficients, yields, delivered yields, and residuals. The distribution of the estimated delivered yields for each reach from the 200 iterations was used to estimate the standard errors in the yields from SPARROW. The estimated confidence interval for delivered yields required explicit consideration of the distribution of residuals in the model, rather than just the summary statistical properties of the residuals. The bootstrap method for incorporating the distribution of the model residuals is based on the empirical distribution of the combined bootstrap-iteration estimate of the modeled component of delivered flux and a randomly selected weighted error from the original monitored values obtained in the original calibration of the model. The 95% confidence interval for the delivered incremental yield from each HUC8 was then estimated using a ratio formulation of the hybrid bootstrap confidence limit (see Shao andTu, 1995, andSchwarz et al., 2006).
The total load and yield of TN and TP from each tributary (> 150 km 2 ) to each Great Lake is provided in Tables S6_TP and S6_TN (attached to the end of the Supporting information). The lower and upper 95% confidence limits for the loads are also in Tables S6_TP   and S6_TN. For each Great Lake, the tributaries are ranked on the basis of their respective loads and yields. It should be noted that SPARROW was only used to estimate the loads from the unmonitored part of each basin, and the confidence limits only represent the uncertainty in that portion of the load (model estimated load); no errors are assumed in the loads estimated with Fluxmaster. The percent contribution of each source of nutrients is also provided in Tables   S6_TP and S6_TN.

Comparison in yields from the Great Lake Basins with those from nearby large river basins
Previous large-scale loading studies (Smith et al., 1997)

D.
The TP yield from the Lake Erie Basin is slightly higher than those from the Upper Mississippi and Ohio River Basins primarily resulting from the Lake Erie Basin having more point-source inputs ( Figure S17). In general, point sources were more important to almost all of the Great Lakes than to the large nearby rivers (7-21% for the large rivers compared to 31-44% for all of the Great Lakes except Lake Superior). Although the land uses vary, TP yields from the Lake Superior and Huron Basins are similar to that of the more intensively farmed Red River Basin. The low yields from the Red River Basin may have resulted from the lower runoff in northern Minnesota and North Dakota than around Lakes Superior and Huron, rather than lower nutrient concentrations in the streams. TP concentrations in streams in the Red River Basin have been shown to be as high as those as many of the intense agricultural areas around the other lakes (Lorenz et al., 2009). The Red River Basin had the highest percentage of TP coming from fertilizer sources, but the lowest percentage coming from point sources.
TN yield from the Lake Erie Basin is almost twice those from the Upper Mississippi and Ohio River Basins, whereas TN yield from the Lake Ontario Basin is similar to those from the large river basins. The higher yield from the Lake Erie Basin was primarily due to it having more input from point sources and streams with shorter travel times and thus less time for losses, such as from denitrification. Point sources of N were more important to most of the Great Lakes than to the large nearby rivers (13-34% for the Great Lakes compared to 3-13% for these large rivers).
TN yield from the Red River Basin was less than those from all of the Great Lake Basins. Low TN yields from the Red River Basin, again, may have been caused by less runoff in northern Minnesota and North Dakota than in the forested areas around Lakes Superior and Huron (Lorenz et al., 2009). The Red River watershed had the highest percentage of N coming from agricultural sources, but the lowest percentage coming from point sources.   Table S6_TP. Total annual phosphorus and nitrogen loads and yields, with confidence intervals, for all tributaries to each lake with a drainage area greater than 150 square kilometers. Tributaries to each Great Lake are ranked based on their relative loads and yields. A value of 1 indicates it has the largest load or relatively largest yield. All loads and yields are adjusted to remove known spatial biases by only predicting areas that are not monitored. This results in very small standard errors and confidence intervals for basins for which a large percentage is monitored.  Table S6_TN. Total annual phosphorus and nitrogen loads and yields, with confidence intervals, for all tributaries to each lake with a drainage area greater than 150 square kilometers. Tributaries to each Great Lake are ranked based on their relative loads and yields. A value of 1 indicates it has the largest load or relatively largest yield.