Sources and Delivery of Nutrients to the Northwestern Gulf of Mexico from Streams in the South-Central United States1

Abstract SPAtially Referenced Regressions On Watershed attributes (SPARROW) models were developed to estimate nutrient inputs [total nitrogen (TN) and total phosphorus (TP)] to the northwestern part of the Gulf of Mexico from streams in the South-Central United States (U.S.). This area included drainages of the Lower Mississippi, Arkansas-White-Red, and Texas-Gulf hydrologic regions. The models were standardized to reflect nutrient sources and stream conditions during 2002. Model predictions of nutrient loads (mass per time) and yields (mass per area per time) generally were greatest in streams in the eastern part of the region and along reaches near the Texas and Louisiana shoreline. The Mississippi River and Atchafalaya River watersheds, which drain nearly two-thirds of the conterminous U.S., delivered the largest nutrient loads to the Gulf of Mexico, as expected. However, the three largest delivered TN yields were from the Trinity River/Galveston Bay, Calcasieu River, and Aransas River watersheds, while the three largest delivered TP yields were from the Calcasieu River, Mermentau River, and Trinity River/Galveston Bay watersheds. Model output indicated that the three largest sources of nitrogen from the region were atmospheric deposition (42%), commercial fertilizer (20%), and livestock manure (unconfined, 17%). The three largest sources of phosphorus were commercial fertilizer (28%), urban runoff (23%), and livestock manure (confined and unconfined, 23%).


SPARROW Model Details
SPARROW is a spatially explicit watershed model that uses a hybrid statistical/mechanistic 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, F i * , is given by The first summation term in the above equation represents the total load delivered to reach i from upstream reaches, where F j ' is the measured load if the upstream reach is monitored or the model-estimated load if it is not. The A( . ) term represents any stream delivery factors that cause load to be lost as it travels along the reach. Within this term, the Z S and Z R vectors represent losses in measured stream and reservoirs, respectively (with corresponding coefficient vectors, Ө s and Ө r ). The second summation term represents the amount of the within-reach load introduced to stream reach i. This term is composed of load originating in individual modeled sources, each source being indexed by n=1,…N S . Each source has a source variable, S n , and a source coefficient, α n , that measures the intensity of source contribution. The function D n ( . ) represents land-towater delivery factors, and, coupled with the coefficient, represents the rate at which the source variable is converted to nutrient mass delivered to streams. The last term in the equation, the function A'( . ), represents the fraction of the load originating in and delivered to reach i that is transported to the reach's downstream node. If reach i is a stream reach (as opposed to a reservoir reach), the N load or P load introduced to reach i from the incremental drainage for reach i is attenuated to receive the square root of the reach's full in-stream delivery. For reservoir reaches, the assumption is made that the nutrient mass receives the full attenuation, which is tabulated as a reach attribute. The multiplicative error term in equation (1), ε i , is applicable in cases where reach i is a monitored reach; the error is assumed to be independent and identically distributed across independent catchments in the intervening drainage between stream monitoring sites.
Statistical terms computed during each SPARROW model run that describe model performance and fit include the sum of squared errors (SSE), the mean square error (MSE), the root mean square error (RMSE), and three coefficients of determination, R 2 , adjusted R 2 , and yield R 2 (all based on log units, Schwarz et al., 2006). The SSE is the squared value of the estimated residual (actual minus predicted), times its weight, and summed over all monitored reaches. The MSE is the SSE divided by the number of degrees of freedom for the error, or DF Error. The DF error statistic is the difference between the number of observations (sampled locations) and the number of degrees of freedom in the model (number of terms used). The RMSE is the square root of the MSE.
The R 2 is the standard coefficient of determination related to the model residual, the measured flux for a particular observation, and the average flux over all observations (described by Judge et al., 1985). R 2 is the portion of the variance in the load data accounted by the independent variables in each model. The adjusted R 2 applies a degree of freedom adjustment to the R 2 statistic. The R 2 and adjusted R 2 terms tend to be large due to the fact that much of the variation in the dependent variable is associated with the size of the drainage area upstream of the monitored reach and that drainage area is highly correlated with source variables; therefore, a high R 2 does not necessarily mean good model fit for smaller basins. Goodness of fit is better described by yield R 2 , which is the logarithm of contaminant yield. Yield R 2 is based on applying a drainage area adjustment to the R 2 to account for scaling effects due to drainage area; therefore, yield R 2 is typically lower than R 2 (Schwarz et al., 2006).

Reach Network Details
There were 8,375 stream reaches in the reach network used for the TN and TP SPARROW models of the LMTG region. The reach network was based on the 1:500,000-scale enhanced River Reach File 2.0 (eRF1_2) reach network for the conterminous United States (U.S.) (Brakebill et al., this issue Sites were deleted if located in bays, estuaries, reservoirs, lakes, or inter-coastal canals, or if they could not be located on the eRF1_2 reach network. There were many cases where different agencies collected data at the same location or in close proximity to each other. Sites were considered collocated if they were within 1,000 meters of each other. There were 43 sites that were also considered collocated that were more than 1,000 meters apart but had no major tributaries entering between them (drainage area ratios were between 0.9 and 1.1). and then a Geographic Information System computer program was used to match the screened sites with a flow gage. Specifically, the sampling station must be on the same stream, and the drainage area ratio of the sampling station to the flow station must be within 0.5 to 2. Using these match criteria, 690 screened sites were matched with a flow gage.
All flow data were from USGS flow gages with the exception of two sites from the U.S. Army Corps of Engineers: flow data for the Mississippi River at Tarbert Landing was used for load calculations for USGS site 07373420 Mississippi River near St.
Francisville, LA, and flow data for the Atchafalaya River near Simmesport, LA, was used for load calculations for USGS site 07381496 Atchafalaya River at Melville, LA. Figure   S1 shows the ranges of drainage area and estimated loads for sites used in the TN and TP models. The drainage areas for the sites in both models ranged from about 23 to nearly 3 million square kilometers (Mississippi River), and the median drainage area for sites in both models was about 1,800 square kilometers indicating that the model calibration sites were representative of small headwater streams to very large rivers.
Loads of TN, TP, or both were estimated for the 690 matched sites using Fluxmaster software (Schwarz et al., 2006). Further manual screening criteria included, but were not limited to, deleting sites if the following occurred: • Large temporal gaps (greater than 5 years) in the concentration data set; • Concentration data were more than 70% censored; and • Flow data were not inclusive of the target year 2002.
Once these final data screenings were completed, there were a total of 468 unique sites, of which there were 344 for TN and 442 for TP load estimation (

Computations Used for Accumulated Delivered Load and Yield Estimates for LMTG Watersheds
A custom SAS program was written to provide accumulated delivered loads and yield estimates for each watershed, and standard errors and confidence limits for these load and yield estimates (Greg Schwarz, U.S. Geological Survey, written communication, July 2009). As part of this process, the analysis first delivers catchment loads to a

MAPS OF FINAL SOURCE AND LAND-TO-WATER DELIVER VARIABLES
Maps of final source and land-to-water delivery terms used in the LMTG TN and TP SPARROW models are based on aggregated datasets (Wieczorek and Lamotte, 2011; unless otherwise noted all spatial data in this article are from this source).. Figure S2 is a plot of land use in the LMTG region. The largest land use grouping is pasture and grassland at 32% followed by forested (19%), cropped (19%), barren and scrub (16%), wetlands (7%), developed residential (5%) and open water (2%). Developed residential land use was used as a surrogate for urban runoff in both models; forested, barren and scrub, and wetlands land use categories were combined and used as a surrogate for background P in the TP model.    Figure S5 -Nitrogen from livestock manure from confined animal feeding operations for the Lower Mississippi Texas-Gulf region. Figure S6 -Nitrogen from livestock manure from unconfined animal feeding operations for the Lower Mississippi Texas-Gulf region.

A B
animal feeding operations used in the TN model. Figure S7 is a plot of the combined livestock manure data used in the TP model. Figure S8 is a plot of total inorganic N from wet deposition detrended to 2002 used as the atmospheric deposition dataset for the TN model. Figure S9 is a plot of channel length for reaches that have streamflow greater than 1.4 cubic meters per second in the LMTG region. Channel length was used as a surrogate for P attached to sediment caused by in-channel erosion in the TP model.
The 30-year average rainfall for the period 1971-2000, which was a highly significant land-to-water delivery term in both models is plotted in Figure S10. The western part of the region is fairly arid (annual rainfall less than about 70 cm total per year), and the eastern part has a humid, subtropical climate with annual rainfall amounts greater than about 100 cm per year. Figure S11 is a plot of overland flow in excess of infiltration, which was highly significant in both models. Overland flow in excess of infiltration was considered a surrogate for runoff potential in both models. Figure S12 is a plot of soil erodibility (or K-factor from the Universal Soil-Loss Equation) for the LMTG region, which was statistically significant as a land-to-water delivery term in the TP model.

MODEL RESULTS
Graphical evidence related to goodness-of-fit for each model is shown in Figure   S13, where the natural log of the observed load is plotted with the natural log of the predicted load. These plots indicate that residuals for both models were normally distributed and homoscedastic (residuals are of constant variance and uniform scatter).
The patterns in both plots for each model also indicate that the model has better accuracy     These data were statistically significant as a land-to-water delivery term in the total nitrogen and total phosphorus Lower Mississippi Texas-Gulf SPARROW models.
Figure S12 -Soil erodibility (K-factor from the Universal Soil-Loss Equation), which was statistically significant as a land-to-water delivery term in the Lower Mississippi Texas-Gulf total phosphorus SPARROW model. A B Figure S14N. -a) Delivered incremental total nitrogen yield and b) delivered incremental total phosphorus yield for the Lower Laguna Madre watershed.