HYDROLOGICAL PROCESSES Hydrol. Process. (2009) Published online in Wiley InterScience (www.interscience.wiley.com) DOI: 10.1002/hyp.7323



# Spatial analysis of instream nitrogen loads and factors controlling nitrogen delivery to streams in the southeastern United States using spatially referenced regression on watershed attributes (SPARROW) and regional
classification frameworks †

Anne B. Hoos 1 

### *
, ‡ 

### and Gerard McMahon 2

, ‡

1 U.S. Geological Survey, Nashville, TN, USA

2 U.S. Geological Survey, Raleigh, NC, USA



### Abstract:

Understanding how nitrogen transport acros s the landscape varies with landscape char acteristics is important for developing sound nitrogen management policies. We us ed a spatially referenced regression an alysis (SPARROW) to examine landscape characteristics influencing delivery of nitrogen from sources in a watershed to stream channels. Modelled landscape delivery ratio varies widely (by a factor of 4) am ong watersheds in the southeastern Unite d States—higher in the western part (Tennessee, Alabama, and Mississippi) than in the eastern part, a nd the average value for the region is lower compared to other parts of the nation. When we model landscap e delivery ratio as a continuous function of l ocal-scale landscape characteristics, we estimate a spatial pattern that varies as a function of soil and climate characteristics but exhibits spatial structure in residuals (observed load minus pr edicted load). The spatial pattern of modelled landscape delivery ratio and the spatial pattern of residuals coincide spatially with Leve l III ecoregions and also with hydrologic lands cape regions. Subsequent incorporation into the model of these frameworks as regional scale variables improves estimation of landscape delivery ratio, evidenced by reduced spatial bias in residuals, a nd suggests that cross-scale processes a ffect nitrogen attenuation on the landscape.
The model-fitted coefficient values are l ogically consistent with the hypothesis that broad-scale classifications of hydrologic response help to explain differential rates of nitrogen attenuation, controlling for local-s cale landscape characteristics. Negative model coefficients for hydrologic landscape regions where the primary flow path is shallow ground water suggest that a lower fraction of nitrogen mass will be delivered to streams; this relation is reversed for regions where the primary flow path is overland flow. Published in 2009 by John Wiley & Sons, Ltd.

Additional Supporting information may be found in the online version of this article.

KEY WORDS nutrients; total nitrogen; spatia lly referenced regression; waters hed models; landscape attenuation

Received 24 November 2008; Accepted 10 March 2009

INTRODUCTION

Riverine and coastal eutrophication arising from natural and human-derived nitrogen sources is an important water-quality issue at scales ranging from a stream reach (Duff et al ., 2008) to large waterbodies such as the Gulf of Mexico (e.g. Bricker et al ., 1999; Alexander et al ., 2000; Mitsch et al ., 2001; Rabelais et al

., 2002; Alexander and Smith, 2006; Scavia and Donnelly, 2007). Nitrogen inputs to the environment are expected to increase with population and economic growth and in response to specific policies and economic forces, such as the promotion of ethanol-based fuels (Booth and Campbell, 2007; Cox, 2007; Jackson, 2007). Improved understanding of the location and amounts of nitrogen

* Correspondence to: Anne B. Hoos, U.S. Geological Survey, 640 Grassmere Park, Suite 100, Nashville, TN 37211, USA. E-mail: abhoos@usgs.gov
† This article is a US Government work and is in the public domain in
the USA. ‡ Anne B. Hoos and Gerard McMahon have authored this work as part of their official duties as employees of the US Government. entering the environment and of environmental characteristics that influence the amount and timing of nitrogen delivery to waterbodies is needed to assess the effects of nitrogen load reduction pr ograms on receiving waters (Destouni et al ., 2006).

The mass of nitrogen delivered from sources in a watershed to a stream channel is determined by the mass of nitrogen inputs to the wate rshed, balanced by mass of nitrogen exported by crop harvesting and also by the

interaction of environmental c haracteristics that influence transport and attenuation of n itrogen during overland andsubsurface transport. We use the term ‘landscape delivery ratio’ (LDR) in this paper to refer to the capacity of a watershed to deliver nitrogen to a stream channel as a result of the rates of controlling processes in that watershed. LDR at the watershed scale of this analysisis expressed as the fraction of nitrogen input that completes the overland and subsurface phase of transport to

the stream channel.

Estimation of the spatial characteristics of LDR requires a means to represent and link the sources of

Published in 2009 by John Wiley & Sons, Ltd.A. B. HOOS AND G. MCMAHON

nitrogen, the processes controlling attenuation (such as sources of energy to fuel denitrification, extent of the

anoxic terrestrial conditions that favour denitrification, availability of denitrifiers, and residence time in soils and ground water), and the amount of nitrogen that reaches a stream (Boyer et al ., 2006). Direct measurement of theseprocesses is difficult at any scale; instead, landscape characteristics that are related to attenuation processes, such as soil, landform, and climate characteristics, have been used in numerous catchment- and basin-scale nitrogen attenuation models to delineate and map the conditions for varying rates of attenuation (Smith et al ., 1997; De Wit, 2001; Arnold and Fohrer, 2005).

Three challenges exist in the use of these landscape characteristics to simulate nitrogen attenuation:

1. Nitrogen attenuation is affected by environmental factors interacting at multiple scales. For example, microbial denitrification in soils is influenced by the interaction of factors such as soil moisture and pore size, which vary at the scale of a hillslope, with regionalpatterns in temperature and precipitation. Understanding the interplay between fine- and broad-scale patterns and processes is key to understanding ecosystem dynamics (Peters et al ., 2008; Lamon and Qian, 2008).2. Spatial patterns in landscape characteristics can be recognized and mapped at multiple scales. Soil characteristics, for example, can be mapped in cells less than100 m on a side, a resolution that can capture relatively fine-grained spatial variation reflected in countyscale soil databases (e.g. the Soil Survey Geographic Database, SSURGO). These same soil characteristics can be represented at a much coarser resolution (e.g. the State Soil Geographic Database, STATSGO). Themap scale and associated spatial resolution of landscape data used for modelling nitrogen attenuation on the landscape should correspond with the scale at which the attenuation processes associated with those factors are expected to vary (McMahon et al ., 2004;

Wolock et al ., 2004). 3. Different factors limit nitrogen attenuation processes in

different areas (Boyer et al ., 2006). Field-scale studieshave discovered numerous examples of spatial variation in the set of controlling factors. For example, Florinsky et al . (2004) demonstrated that topographicproperties, such as land-surface slope and relative elevation, influenced denitrification rates in relatively wet soils, but had little or no influence on denitrification in dry soils. At the watershed scale, soil depth influences nitrogen delivery to strea ms through denitrification along groundwater flow paths (Clement et al ., 2002) and therefore a lower LDR is expected for deeper soils.In regions where mass transport is predominantly overland, however, variation in soil depth has less effect on LDR. These examples suggest that LDR cannot be adequately modelled a s a continuous function of landscape characteristics because this approach cannot capture spatial variation in the set of controlling factors and their quantitative rela tion to nitrogen attenuation. Regional frameworks, such as physiographic, geologic, or ecological regions, may be useful in modelling the effects of relatively broad-scale spatial processes that affect nitrogen attenuation. The area within any region

in such a framework has a characteristic mosaic of landscape features and processes that is presumably distinctive from adjoining regions; in effect, the framework represents a hypothesis that each region is distinct fromother regions in terms of either a single type of process, such as a characteristic hydrologic response, or amore integrated set of processes, such as those that produce characteristic biota in undeveloped areas (McMahon

et al ., 2004). The landscape mosaic is defined by interactions that occur at multiple scales among human and natural influences, such as e conomic, climatic, and physiographic processes. These interactions tend to organizethe landscape into a distinct spatial mosaic of recognizable landscape features (Omernik, 1987, 2004; Bailey, 1988, 2004).

Wolock et al . (2004) developed hydrologic landscaperegions (HLRs) to differentiate areas with characteristic hydrologic response and primary flow paths. The US

Environmental Protection A gency (USEPA) developedecological regions that reflect a hypothesis that ecological processes, such as denitrification, may vary from one region of the country to anothe r due to broad geographicvariation in a matrix of factors, such as topography, geology, soils, land cover, land use, and climate (Omernik,1987). In either case, if a regional framework differentiates areas with characteristic biotic and abiotic processesthat influence the transformation of nitrogen in the environment, such a framework may provide information useful in fitting and applying nutrient attenuation models (McMahon et al ., 2004).

In this paper we examine an approach to using regional- and local-scale landscape variables to predict nitrogen LDR as part of a broader regression analysis [spatially referenced regre ssion on watershed attributes (SPARROW)] of nitrogen input and attenuation in 321

basins in the southeastern United States. SPARROW uses nonlinear regression to quantify the relation among nitrogen inputs or sources, attenuating environmental characteristics, and measured instream nitrogen load. The term instream nitrogen is used in this paper to refer to the mass of total nitrogen transported in stream; that is, the total mass of dissolved and suspended fractions of inorganic and organic nitrogen. The regression analysis also provides a prediction e quation for mean annual loads and concentrations of nitrogen for each stream reach in the model area. The SPARROW model has been applied

to assess the effect of sources and attenuation factors on stream nutrient loading at the national scale (Smith et al

., 1997; Alexander et al ., 2000) and for individual regions and river basins, such as the Mississippi River Basin (Alexander et al ., 2008), the Chesapeake Bay watershed (Preston and Brakebill, 1999), New England river basins (Moore et al ., 2004), eastern North Carolina river basins (McMahon et al ., 2003), and Tennessee, Kentucky, and Alabama river basins (Hoos, 2005). Most of these model

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

applications, with the exception of Chesapeake Bay,

estimate LDR as a continuous function of local-scale landscape characteristics.

We incorporate regional-scale landscape variables in the SPARROW model to address two questions:

1. How does incorporation of regional landscape variables affect model error?

2. How does inclusion of such variables affect model

predictions of instream nitr ogen load, nitrogen delivery budgets, and the relative importance of nitrogen

sources?

We hypothesize that nitrogen delivery from the landscape to the stream is affected by processes operating at multiple spatial scales and that SPARROW models that include regional-scale varia bles, such as hydrologic landscape regions or ecoregions, will have a better fit between predicted and observed values of instream nitrogen load,and less spatial structure in model residuals, than models without regional variables. Regions, characterized in the model using nominal variables, are hypothesized to influence landscape nitrogen processing in a way that isdistinct from measured, local-scale landscape characteristics summarized at a stream-reach scale, such as soil permeability, depth to bedrock, and annual precipitation.

The area included in this investigation includes the river basins draining to the south Atlantic Coast, the eastern Gulf Coast, and the Tennessee River (Figure 1), referred to collectively as the SAGT river basins. This area is one of the eight large geographical regions across the nation (referred to as ‘major river basins’) identified by the National Water-Quality Assessment (NAWQA) program of the US Geological Survey as the basis for assessments of status and trends. The NAWQA program has integrated the SPARROW modelling approach in the interpretation of nutrient transport in six of these major

river basins.

METHODS

The SPARROW model uses a nonlinear regressionequation to describe the relation between spatially referenced watershed and channel characteristics (predictors) and instream load (response) (Schwarz et al ., 2006, p. 2). Model input consists of spa tially referenced datasets representing stream-channel networks, nitrogen sources, physical watershed character istics, and observations of instream nitrogen load. The i nput datasets developed for model application in the SAGT area are presented in the work by Hoos et al . (2008).

The model

For each reach in a hydrologic network, SPARROW predicts long-term mean a nnual instream nitrogen load as a function of nitrogen sources, nitrogen attenuation on the landscape, and nitrogen losses that occur within the

Figure 1. Location of the SPARROW model study area and Level III ecoregi ons within the southeastern United S tates, and instream load estimates for 2002 used to calibrate the nitrogen model

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

stream. Conceptually, the instream nitrogen load or flux at the downstream node of a reach can be expressed as the sum of two components:

L instream i D L catchment i C L upstream i ! 1 "

where L instream i D instream load at downstream node

of reach i ; L catchment i D load originating within the

catchment for reach i and delivered to the downstream

node of reach i ; and L upstream i D load generated within catchments for upstream reaches and transported to the

downstream node of reach i via the stream network.

The load originating within the catchment for reach i ( L catchment i ) is determined by:

L catchment i D N S !

n D 1 S n,i  ̨ n D
n " Z D i ; # D # A " Z S i , Z R i ; # S , # R #

! 2 "

where n , N
s D source index where N s is the total number of individual sources;

S
n D vector of source variables (for example, a measurement of mass placed in the watershed, or the area of a particular land cover); and

 ̨ n D vector of coefficients, estimated by the model

, in units that convert source variable units to flux units. For land-applied sources,  ̨ n is the model estimate of the average LDR (LDR avg ), across all catchments in the model area. For land-applied sources represented by characteristics other than mass input (for example, by impervious surface area),  ̨ expresses the conversion of source units to mass applied to the watershed as well as

the LDR avg for the source.

D
n !

Ð " D the delivery variation factor (DVF), defining the variation among catchments in nitrogen landscape attenuation processes and therefore in LDR. The DVF is modelled as a series of exponential functions of physical landscape characteris tics that influence nitrogen

attenuation. The DVF for catchment i is multiplied by the

LDR avg for source n (that is, by  ̨ n ) to calculate LDR i,n

.

Z D

D vector of physical landscape variables (for example measured landform or soil characteristics); and

# D D vector of coefficients, estimated by the model , for the physical landscape variables.

A !

Ð " D the stream delivery function, representing the result of attenuation processes acting on flux as it travels along the stream channel. Modelled as first-order decay, the stream delivery function defines the fraction of flux originating in and delivered to reach i that is transported

to the reach’s downstream node.

Z S and Z R

D vectors of measured stream and reservoir variables, respectively (ex amples include stream-water depth or velocity and reservoir areal hydraulic loading);

and # S and # R D vectors of coefficients, estimated by the model , for the stream and reservoir variables, respectively.

The DVF allows the model to simulate variation in LDR among catchments. The median value of DVF for all catchments in the model area is approximately 1 when the DVF is modelled as exponential functions of the departure of the landscape variables from their respective means. DVFs greater than 1 indicate a larger fraction of nitrogen reaching streams than the median for the model area; DVFs less than 1 indicate a smaller fraction of nitrogen reaching streams than the median for the model

area.

The second component in Equation (1), the flux entering reach i from upstream reaches, is the sum of the flux from any upstream catchment ( L catchment i

" 1 , L catchment i

" 2 , etc.) adjusted for losses due to stream and reservoir attenuation processes acting on flux along the reach pathway to and including reach i . For headwater reaches, Equation (1) is simplified to include only

the Lcatchment i term. More information about the model form and assumptions is available in the work by Schwarz

et al . (2006).

The stream network

The data framework for SPARROW is the network of stream- or reservoir-reach segm ents and associated catchments. The hydrologic network used for the SPARROW

model of the SAGT river basins is the Enhanced River

Reach File 2

Ð 0 (ERF1 2), based on USEPA’s 1 : 500000-scale Reach File 1 (RF1) with enhancements to support national and regional-sc ale water-quality modelling (Nolan et al ., 2002). The digital datasets describing the

ERF1 2 network for the SAGT area are presented by

Hoos et al . (2008) and include stream discharge, time of travel, and, for reaches associated with a reservoir,reservoir areal hydraulic loading (ratio of reservoir outflow to surface area). The NHDPlus digital network (US Environmental Protection Agency and US Geological Survey, 2006) represents the stream network at a finer spatial scale; its implementation as infrastructure for a SPARROW model will require additional attributes that are currently (2008) in development (Terziotti and Hoos, 2008).

Annual instream nitrogen load

Measurements of nutrient water quality at stream monitoring sites collected by Federal, State, and local agencies during 1975–2004 were used to develop observations of mean annual nitrogen load as the response variable in the SPARROW regression equation. Mean annual load was estimated as the product of daily streamflow and estimated daily concentration, which was modelled from nutrient water-quality data and streamflow data. Siteselection criteria and estim ation methods used to develop the set of estimates of mean annual nitrogen load for 2002 are described in the study by Hoos et al . (2008). Of the 637 stations with estimates of 2002 nitrogen load in the SAGT SPARROW model area, 331 stations (Figure 1) were placed on the SAGT ERF1 2 digital segmented network and used to calibrate a nitrogen SPARROW model;

the other 306 sites were located on tributaries too small to be represented in the relatively coarse 1 : 500 000 ERF1 2 network or lacked independent information for calibration due to proximity (for example, within 1 km) to another site with a nitrogen load estimate.

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Figure 2. Hydrologic landscape regions in the southeastern United States, in relation to Level III ecoregion boundaries (ecoregions are identified o n figure 1)

The number of stations available for model calibration affects the statistical power of the regression: models with more stations generally have greater power to detect the effect of an explanatory variable on instream load. The SAGT calibration set (331 stations) is large compared with previously developed regional SPARROW nitrogen models (with fewer than 80 stations) and corresponds to a calibration site density for the SAGT model area (780 000 km 2 ) of about 2500 square kilometers per site. About one-half (163) of the stations are situated downstream from one or more stations; that is, about

one-half of the monitored basins are nested within other monitored basin areas. Analyses by Schwarz et al . (2006, p. 69) have shown that some degree of nesting may raise the precision of the estimated coefficients; precision is bounded above, however, so that at some point greater nesting does not improve the accuracy of the estimate.

Nitrogen source and landscape variables

Nitrogen source variables tested in the model include 2002 estimates of atmospheric deposition, commercial fertilizer applied to agricu ltural land, animal manure from livestock production, point-source discharge of wastewater, impervious surface area, population density, and urban, agricultural, and forested land cover (Hoos

et al ., 2008). Impervious surface area, population density, and land cover are surrogate indicators of nitrogen mass, which can be considered to be proportional to the actual mass loadings generated by human activities.

Local-scale landscape variables tested in the model include hydrologic soil group (infiltration rate), soil permeability, available water holding capacity of the soil, average percentage of saturation overland flow in total streamflow, soil sand content, soil organic carbon content, soil erodibility, average depth to bedrock, wetland area, land-surface slope, percentage of flatland, and mean annual precipitation and temperature. These data vary over relatively small spatial scales.

One regional-scale variable tested in the model is the set of HLRs developed by Wolock et al . (2004) to provide a framework for distinguishing primary hydrologic flow path and response characteristics in 20 distinct regions across the United States, 14 of which are within the SAGT area (Figure 2). The boundaries of the HLRs weredefined by multivariate-analysis to describe the combined effect of landscape and climate characteristics (soil texture, aquifer permeability, land surface elevation andslope, precipitation, and potential evaporation) on hydrologic response. The HLR is a geographically independentframework in that a region is defined by a combination of characteristics and may be composed of several non-contiguous areas.

Another set of regional-scale variables tested are the USEPA Level III ecoregions, developed to assist resource managers and researchers in structuring programs related to water quality and biological criteria (Omernik, 1987; 1995; 2004; McMahon et al ., 2001). The boundaries of Level III ecoregions (Figure 1) are also defined bymultivariate analysis of a set of biotic and abiotic factors. Regional boundary decisions, however, are arrived at heuristically through extensive discussions among a diverse team of experts, followed by field verification. Ecoregions distinguish rela tively homogeneous spatial mosaics of biotic and abiotic resources, ecosystems, and

human influences.

Level III ecoregions and HLRs have a broad spatial coincidence, reflecting overlap in the set of the factors—topography, climate, soils—common to the

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

derivation of both frameworks (Figure 2). Differencesmay be attributed to differences of scale in the derivation of the frameworks—HLRs have smaller mapping units—as well as differences in objectives of the two mapping efforts—Level III ecoregions reflect patterns in many biotic and abiotic factors hypothesized to influence ecosystem functions, while HLRs distinguish patterns infour factors that influence hydrology. Some of the finerscaled detail in the HLRs is captured in the Level IV ecoregions (US Environmental Protection Agency, 2007).

Model specification and testing

We tested a variety of model specifications to evaluate which sources and landscape characteristics among those that can be reasonably represented and described within the construct of SPARROW are important in controlling nitrogen transport. Variables identified as significant in explaining nitrogen transport (using the model-computed

p -value as the test for significance) were retained and

combined with additional attributes in a series of model runs, until a model specification was achieved that was optimal in terms of model fit, model-estimated coefficients, and residual plots.

Certain sites within the SAGT river basins were

excluded from the calibration set due to concern that actual hydrologic boundaries upstream from a site did not correspond with the apparent watershed determined from surface topography; thus, the SPARROW approach of explaining instream loads based on watershed attributes would be inappropriate. River basins identified with this concern included those in south Florida (where surface-water flow paths have been extensively altered) and the Oklawaha, Crystal, Lower Sante Fe, Lower Suwanee, St. Marks, and Chipola river basins in central and northern Florida (where flow exchange with the underlying regional aquifer may contribute substantial nitrogen influx to and outflux from the surface-water basins; Rumenik, 1988; Miller, 1990). Ten of the 331 sites with estimates of 2002 nitrogen load were thus excluded, reducing the set calibration sites to 321. River

basins in south Florida are excluded from both model prediction and calibration; model predictions for the river basins in central and north Florida are presented but may

be less reliable due to this unmodelled component of flux.

RESULTS AND DISCUSSION

To test the effect of regional landscape variables on model error, we developed three different model specifications, presented here as Models A, B, and C. In Model A, the

DVF was simulated as a continuous function of localscale landscape variables; in Models B and C, the modelled relation between DVF and landscape variables is modified for each HLR and ecoregion, respectively. The attributes selected as explanatory variables in the threemodels are listed in Table I, along with estimated coefficients and calibration error s tatistics. To test the effect of regional landscape variab les on model predictions, we compare nutrient budgets estimated by Models A, B, and C for an individual river basin (Tombigbee River) in the

southeastern United States.

Source variables and coefficients, Model A

Five source variables were selected in the specificationfor Model A: atmospheric deposition, commercial fertilizer applied to agricultural land, animal manure fromlivestock production, impervious surface area, and pointsource discharge of wastewater. The model-estimated

coefficients  ̨ (Table I) represent, for the land-applied sources, the LDR avg for each source; LDR for a specific catchment ranges upward or downward depending on the model-estimated DVF for the catchment. For example,

the estimate of  ̨ for commercial fertilizer applied to agricultural land, 0

Ð 13 kg/kg, means that, for each kilogramof nitrogen applied to the land surface in fertilizer, estimated from fertilizer sales data, the model predicts that

0

Ð 13 ( C /

" 0

Ð 03) kg is delivered to the adjacent stream

channel from a catchment with estimated DVF of 1. The

balance of 0

Ð 87 ( C /

" 0

Ð 03) kg is removed either at the point of application (for example by crop harvest) or during land-phase transport. Smaller estimates of  ̨ (and

thus of LDR avg ) for fertilizer and animal manure compared to atmospheric deposition correspond to modelestimated, smaller fractions of the measured input fromthese sources entering land-phase transport. The modelfitted values are scale-specifi c and describe the fractions delivered to catchments similar in size to the segmented

reach-catchment network used to calibrate the model. The average catchment size is 87 km 2 in the 1 : 500 000-

scale-based model network.

The estimate of  ̨ for wet deposition of total inorganic nitrogen, 0

Ð 50 kg/kg, is not directly comparable

to the estimates of  ̨ from other 1 : 500 000-scale-basedSPARROW models that use measures of wet deposition of nitrate-nitrogen (1

Ð 0 kg/kg from the Chesapeake Bay watershed model, Preston and Brakebill, 1999;

0

Ð 69 kg/kg in the national model, Alexander et al ., 2008).Wet deposition of nitrate-nitrogen alone was used to characterize atmospheric nitrogen deposition in the Chesapeake Bay and national models to avoid uncertainties inquantifying contributions from ammonia emissions associated with other sources in the model (animal manureand commercial fertilizer). In contrast, we use total inorganic nitrogen to characteriz e atmospheric nitrogen deposition in the SAGT model because estimates of wet deposition of total inorganic nitrogen, animal manure, and

fertilizer are not correlated s ignificantly (coefficients of

determination r 2 < 0

Ð 1) in the SAGT area, and because this approach permits quantification of the contribution of total atmospheric inputs for the watershed to the stream.

Estimates of  ̨ for fertilizer and animal manure, 0

Ð 13

and 0

Ð 05 kg/kg, respectively, are smaller than estimates from SPARROW models that used comparable input measures (the Chesapeake Bay watershed and national models). Smaller values of SPARROW-estimated LDR avg for the southeastern United States compared with the

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USATable I. Calibration results for three nitrogen SPARROW models (Models A, B, and C) for the southeastern U.SPredictor variableModel-fitted coefficientModel A Model B Model CDescription and unit Unit Value Standard error p Value Standard error p Value Standard error pS (Eq. 2), source input variablesa(Eq. 2)aNitrogen mass in permitted wastewater discharge,2002; kg/yrkg/kg 0⋅80 0⋅10<0⋅005 0⋅79 0⋅09<0⋅005 0⋅79 0⋅10<0⋅005Wet deposition of inorganic nitrogen (ammonia andnitrate), detrended to 2002; kg/yrkg/kg 0⋅50 0⋅05<0⋅005 0⋅50 0⋅05<0⋅005 0⋅53 0⋅05<0⋅005Area of impervious surfaces, 2001; kg/yr kg/km21990 666 0⋅01 2470 649<0⋅005 2230 399<0⋅005Nitrogen mass in commercial fertilizer applied toagricultural land, 2002; kg/yrkg/kg 0⋅13 0⋅03<0⋅005 0⋅11 0⋅02<0⋅005 0⋅11 0⋅02<0⋅005Nitrogen mass in manure from livestock production,2002; kg/yrkg/kg 0⋅05 0⋅02 0⋅02 0⋅05 0⋅02 0⋅01 0⋅04 0⋅01<0⋅005ZD(Eq. 2), physical landscape variablesqD(Eq. 2 and 3)aLn of soil permeability, low value; ln of cm/day−0⋅12 0⋅065 0⋅04−0⋅09 (notretained) 0⋅17−0⋅40 0⋅069<0⋅005Ln of depth to bedrock; ln of cm−0⋅50 0⋅16<0⋅005−0⋅284 0⋅169 0⋅035−0⋅094 (notretained) 0⋅710Ln of mean annual precipitation; ln of mm 1⋅4 0⋅32<0⋅005 1⋅2 0⋅32<0⋅005 1⋅0 0⋅27<0⋅005Fraction of catchment in HLR 2; dimensionless−0⋅29 0⋅15 0⋅02Fraction in HLR 7−0⋅31 0⋅14 0⋅01Fraction HLR 16b−0⋅14 0⋅14 0⋅11Fraction in HLR 40⋅26 0⋅12 0⋅02Fraction in HLRs 6, 9, or 110⋅26 0⋅11 0⋅02Fraction in Ecoregion 45−0⋅38 0⋅09<0⋅005Fraction in Ecoregion 66b−0⋅06 0⋅11 0⋅56Fraction in Ecoregion 67b0⋅17 0⋅12 0⋅13Fraction in Ecoregion 680⋅50 0⋅15<0⋅005Fraction in Ecoregion 71b0⋅02 0⋅16 0⋅89Fraction in Ecoregion 750⋅52 0⋅11<0⋅005ZS(Eq. 2), stream variables:qS(Eq. 2)Time of travel in reach segments where meanQ per day 0⋅23 0⋅11 0⋅02<2.8 m3/s; dayc0⋅14e0⋅05<0⋅005 0⋅12e0⋅04<0⋅005Time of travel in reach segments where meanQ>2.8 and<28 m3/s; daycper day 0⋅13 0⋅05 0⋅005

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHONTable I. (Continued)Predictor variableModel-fitted coefficientModel A Model B Model CDescription and unit Unit Value Standard error p Value Standard error p Value Standard error pTime of travel in reach segments where meanQ>28 m3/s; dayc,dper day 0⋅00 0⋅02 0⋅43 0⋅014 0⋅02 0⋅26 0⋅02 0⋅03 0⋅52ZR(Eq. 2), reservoir variable:qR(Eq. 2)Inverse of areal hydraulic loading; yr/m m/yr 13⋅1 2⋅8<0⋅005 10⋅7 2⋅2<0⋅005 10⋅3 3⋅3<0⋅005Model A Model B Model CStandard Error of the Estimate (SEE), computed asroot mean square error (RMSE)0⋅34 0⋅32 0⋅32Coefficient of determination (R2) of load estimate 0⋅96 0⋅96 0⋅96Coefficient of determination (R2) of yield estimate 0⋅68 0⋅72 0⋅73Moran’sIof log residuals (z-score and significance ofdeparture from the null hypothesis2⋅8 0⋅002 0⋅67 0⋅25 0⋅68 0⋅25of spatially independent residuals)Variables selected for these model specifications were significant at alphaD0.05, except where coefficient value is shown with shading; HLR, hydrologic landscape region, described in Wolock and others (2004);Level III ecoregions described in Omernik (1987); kg,kilogram; km2, square kilometer; Ln, natural logarithm transformation; cm, centimeter; mm, millimeter; m, meter; s, second; yr, year; not retained, variable notretained in model specification; Moran’s I weights were calculated as the inverse of the distance between each pair of sites.aModel variables and coefficients are defined and explained in equations 2 and 3. ̨(Eq. 2) is interpreted for land-applied sources as average landscape delivery ratio.!D(Eq. 2 and 3) defines the relation betweenthe spatial variation of landscape delivery ratio (spatial variation is modeled as the delivery variation factor, DVF) and the landscape variables.bThe variables representing fraction of catchment in HLR 16 and in Ecoregions 66, 67, and 71 are retained despite lack of significance at alphaD0.05 in order to represent the complete regional classificationframework. HLR 1, and Ecoregions 63 and 65, have been excluded from Models B and C (respectively) to avoid colinearity.cEvaluated only for reaches in the indicated size class, otherwise set to 0.dThe variable time of travel in reach segments with mean discharge>28 m3/s is retained despite lack of significance at alphaD0.05, in order to complete the description of nitrogen removal across the full rangeof stream sizes.eWhen stream classes<2.8 and 2.8–28 m3/s are specified separately in Model B, fitted coefficients are 0.15 and 0.14, respectively, and the coefficient for the<2.8 m3/s class is not significant at alphaD0.05. Thetwo classes were therefore combined to a single variable (<28 m3/s) (Similarly for Model C.).

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Chesapeake Bay watershed or the national average may be explained by higher rates of plant uptake or higher microbial activity leadin g to greater denitrification or immobilization on the landscape. Alternatively, the smaller coefficient values may reflect the SAGT model

treatment of ammonia losses related to volatilization from animal manure as atmospheric deposition, rather than

as direct losses to the stream. The estimate of  ̨ for fertilizer, 0

Ð 13 kg/kg, is similar to the reported LDR avg for fertilizer, 0

Ð 10 kg/kg, estimated from nutrient budgetstudies in small agricultural catchments in the southeastern coastal plain (Lowrance et al ., 1985), althoughthe budget-based estimates are not directly comparable because of scale differences (average catchment size

about 1 km 2 ).

The estimate of  ̨ for impervious surface area, 1990 kg/km 2 , associates the mass of nitrogen deliveredto a stream channel with unit area of impervious surface. In this application, impervious surface area serves as a surrogate for many diffuse sources of nitrogenin urban areas, such as vehicle emissions, lawn fertilizer, and onsite sewage-disposal systems. The estimate of 1990 kg/km 2 for the SAGT SPARROW model cannot be compared with values estimated for other SPARROW models because estimates of impervious surfacearea have not been widely available as broad-scale continuous measurements for testing in previous SPARROW models. In previous models (Chesapeake Bay watershed and New England river basins), urban and suburban land classifications (US Geological Survey, 2001) have been used as a surrogate for diffuse nitrogen sources. In the SAGT SPARROW area, urban–suburban land area is about 3 times greater than impervious surface area. The

estimates of  ̨ for the Chesapeake Bay watershed (Preston and Brakebill, 1999) and New England river basins (Moore et al ., 2004) for urban–suburban land area, 784 and 895 kg/km 2 , respectively, are about one-third the magnitude of the SAGT SPARROW estimate of  ̨ for impervious surface area (1990 kg/km 2 ). The estimate of nitrogen mass transported to streams from urban areas, calculated as the product of source variable and  ̨ , therefore appears to be comparable among these three models. The impervious surface area variable is preferred over the

urban–suburban land area var iable to quantify nitrogen mass transported from urban areas because it represents a single type of physical surface and thus the estimate of

 ̨ is interpretable as an export coefficient from a specific

land cover.

An  ̨ value of 1 kg/kg is expected for point-source discharges of wastewater because point sources are discharged directly to streams. The model estimate (0

Ð 80 kg/kg) and the upper bound of the 90% confidence interval (0

Ð 97 kg/kg) are both less than 1 kg/kg, suggesting either that mass inputs are overestimated or

that modelled instream atte nuation cannot fully account for point-source attenuation that may include higher ratesin local mixing zones near the discharge point. This finding is consistent with the SPARROW model for North Carolina basins (McMahon et al ., 2003) and with models that incorporate different stream-loss functions for point sources versus diffuse sources (Destouni et al ., 2006).

Variables and coefficients describing spatial variation in landscape delivery ratio, Model A

The Model A estimates of the coefficients ! D (Table I) represent the best fit between the inferred gradient across calibration sites in the DVF and the gradient across calibration sites in the combination of tested landscape variables. Soil permeability, depth to bedrock, and mean annual precipitation were the most significant predictors

of DVF from the set of tested characteristics.

The DVF is modelled in this application as:

D "

Ð # D exp ! M D "

m D 1 ω nm Z D m i ! Dm #

, " 3 #

where

M D D the number of landscape variables;

Z D

m i D landscape variable m for catchment i (expressed as the departure from its mean value);

! Dm D the corresponding coefficient, estimated by the

model ;

ω nm D an indicator variable that is 1 if landscape

variable m affects source n and is 0 otherwise. ω nm is set to 1 for the interaction of land-applied sources (atmospheric deposition, fertili zer applied to agricultural land, animal manure, impervious surface area) with

all M D landscape variables and is set to 0 for the interaction of point-source discharge of wastewater with

all M D landscape variables (Schwarz et al ., 2006). The

coefficient vector ! D m defines the relation between landscape variables and LDR by quantifying the marginal change in DVF (and therefore in LDR) for a marginal change in each landscape variable. Because the variables are log transformed, the magnitude of the coefficient defines the percent change in DVF given a percent change in the landscape variable.

The Model A estimates of ! D (

" 0

Ð 12,

" 0

Ð 50, and 1

Ð 4 for log-transformed soil permeability, depth to bedrock, and mean annual precipitation, respectively, Table I) define the separate effect of each landscape variable on the DVF (Figure 3). For example, the estimate of ! D for mean annual precipitation, 1

Ð 4, means that a 1%difference between two catchments in mean annual precipitation (1010 vs 1000 mm) causes a 1

Ð 4% difference in DVF (0

Ð 667 vs 0

Ð 658 provided the other landscapevariables remain constant at their mean values). The relation of the individual landscape variables to the DVFand LDR can be inferred from the sign of their coefficients. Negative coefficients for soil permeability and depth to bedrock indicate that the DVF and LDR are higher for catchments with lower soil permeability and depth to bedrock; the positive coefficient for mean annual precipitation indicates that the DVF and LDR are higher for catchments with higher precipitation. Because thevariable depth to bedrock (derived from the US Department of Agriculture’s STATSGO dataset as described in the study by Hoos et al ., 2008) is truncated at 150 cm,

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

Figure 3. Variation of delivery variation factor and landscape delivery ra tio for manure from livestock operations with mean annual precipitation a nd hydrologic landscape region for (A) Model A and (B) Model B

the estimated coefficient for depth to bedrock defines the percent change in DVF for depths less than 150 cm only.

One-way analysis of variance (ANOVA) and Tukey multiple comparison tests were used to evaluate the degree of correspondence between the distribution of

Model A estimated DVF and the boundaries of LevelIII ecoregions and HLRs (Table II). The Level III ecoregion framework explains 40% of the total variation in estimated DVF, with distributions sufficiently different (at the 5% significance level) between some ecoregions to enable division into three distinct groupings. The HLR framework explains 42% of total variation in DVF, but

differences in distributions between HLRs were not as distinct as for ecoregions. These comparisons provide

evidence that the estimated distribution of DVF and LDR for individual catchments has a spatial structure that is moderately correlated with the regions within both

frameworks.

Stream and reservoir loss coefficients, Model A

Instream removal of nitrogen is influenced by processes such as denitrificatio n and sedimentation forwhich mean annual rates likely decrease with increasing stream size (Alexander et al ., 2000; Schwarz et al

., 2006). Stream delivery of nitrogen is modelled in this application as a first-order decay function:

exp !

! Z S " S #

,

where Z S

D stream travel time, summed along the reach pathway separately for each stream size class; and " S D first-order stream loss coefficients for each stream size class, in units of inverse time, estimated by the model

.

Although channel characteristics other than stream size, such as channelization, ditching/draining, condition of riparian vegetation, also may influence instream loss rates, these are not included in the model due to lack of regionally extensive and consistent datasets.

The estimated loss-rate coefficients for the SAGT

Model A are 0

Ð 23 per day of travel for small (mean

annual flow < 2

Ð 8 m 3 /s) streams, 0

Ð 13 per day of travel for intermediate (mean annual flow 2

Ð 8–28 m 3 /s) streams, and negligible ( < 0

Ð 005, p -value D 0

Ð 43 and thus not significant at the 5% significance level) for large (mean

annual flow > 28 m 3 /s) streams. The SPARROW-fitted inverse relation between nitrogen loss-rate coefficient and stream size is consistent with the concept that attenuation decreases with increasing stream size, and with the results of experimental studies (Howarth et al ., 1996;

Mulholland et al ., 2002; Gibson and Meyer, 2007) and of other SPARROW nitrogen models (Preston and Brakebill,

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USATable II. Distribution of Model A estimated delivery variation factor and residuals for calibration sites, with respect to level III ecoregions and hydrologic landscape regionsLevel III ecoregionHydrologic landscape region with similar extent to corresponding ecoregionEcoregion numberand nameNumber (density) ofsites with upstreamarea pre-dominantlyin ecoregionMean value ofModel Aestimated deliveryvariation factorTukeyclassMean valueof Model AresidualHLR name anddescriptionPrimaryhydrologicflow pathNumber (density)ofsites with upstreamarea pre-dominantlyin HLRMean value ofModel Aestimated deliveryvariation factorTukeyclassMean valueof Model Aresidual66-Blue Ridge 22 (0.6) 1.2 B!0.1016S 33 (0.6) 1.2 A!0.10(Mt PS IB)45-Piedmont 125 (0.8) 0.93 C!0.067S 56 (0.8) 0.97 AB!0.10(Pn PS IB)2S,D 36 (0.3) 0.82 B!0.12(Pn PS PB)65-Southeastern Plains 72 (0.2) 0.96 C 0.024S,D 29 (0.2) 1.2 A 0.16(Pn PS PB)a6O 4 (0.1) 0.94 AB 0.12(Pn IS IB)63-Middle Atlantic Coastal Plain 4 (0.1) 0.82 C 0.111S,D 46 (0.1) 0.79 B 0.08(Pn PS PB)75-Southern Coastal Plain 36 (0.1) 0.68 C 0.1171-Interior Plateau 6 (0.2) 1.3 B 0.1711O 9 (0.4) 0.92 AB 0.06(Pt IS IB)68-Southwestern Appalachians 10 (0.3) 1.5 A 0.15b67-Ridge and Valley 20 (0.4) 1.1 B 0.179O,D 23 (0.3) 1.2 A 0.13(Pt IS IB)Multiple 260.96C!0.07Multiple Multiple 830.96AB0.01ANOVA r20.45 0.06 0.43 0.08[Model residual, natural log of observed load minusnatural log of model-predicted load, negative values indicate overprediction; Level III ecoregions described in Omernik, 1987; density of sites expressed as ratioof number of sites to area of region, in square kilometers x 103; Multiple, no single region dominant in upstream area; Hydrologic landscape regions (HLRs) described in Wolocket al, 2004: MtDmountains, PnDplains, PtDplateau, PSDpermeable soils, ISDimpermeable soils, PBDpermeable bedrock, IBDimpermeable bedrock; primary hydrologic pathway: SDshallow ground water, ODoverland flow, DDdeepground water; delivery factor estimated based on soil permeability, depth tobedrock, and mean annual precipitation; negative value of residuals indicates model overprediction; Tukey class, means with the sameTukey class letter are not significantly different at 0.05; ANOVA r2, coefficient of determination from one-way analysis of variance, interpreted as the percentage of total variation of the variable explained by theregional framework].aEcoregions 45 and 65 contain small areas of HLR 6.bEcoregion 68 does not correspond closely with an HLR or set of HLRs.

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

1999; Alexander et al ., 2000; Alexander et al ., 2002;

McMahon et al ., 2003).

The estimated stream loss coefficients are similar to values estimated for these stream classes by the national SPARROW models (range of 0

Ð 05–0

Ð 32 per day, Alexander et al ., 2008; 0

Ð 38 per day, Smith et al ., 1997) butlarger than observations from experimental studies (typically less than 0

Ð 15 per day, Howarth et al ., 1996). Wesuggest that this discrepancy is partly due to comparing rates modelled for mean annual flux with measured instantaneous rates, and partly due to the coarse-scale (1 : 500 000) hydrography used to define flow-path length

and time of travel for these SPARROW models. Timeof-travel estimates for individual reaches in SPARROW models are calculated from mean water velocity and flow-path length estimated from digital line graph (DLG) hydrography. The 1 : 500 000-scale hydrography used for the SPARROW models cited above yields estimates of flow-path length and corresponding time of travel between fixed points that are biased low (due to less detail of flow-path curvature) compared to measured time of travel used in reach-scale experimental studies. This

low bias in estimated reach time of travel for the models results in a proportionally high bias in estimated loss rate.

The number and the boundaries of the stream sizeclasses vary among the SPARROW models. The boundaries of stream size classes for each model were selected

to ensure that a sufficient number of load observations in the calibration set is represented in each flow class. In the SAGT SPARROW model-calibration set, 34 load observations (about 10%) represent sites for which all contributing stream reaches are in the small ( < 2

Ð 8 m 3 /s) class, 117 load observations (a bout half) represent siteswith contributing stream reaches in the small and intermediate size classes, and 170 load observations represent sites with contributing reaches in all classes.

We simulated reservoir processing of nitrogen as the first-order mass transfer rate expression

1 / ! 1 C Z R " R #

,

where

Z R

D the inverse of the reservoir attribute areal hydraulic loading, in units of time per distance; and

" R D the reservoir loss coefficient, in units of distance per time, estimated by the model

.

The model-estimated reservoir loss coefficient of

13

Ð 1 m per year, describing the mean water column length from which nitrogen is removed annually, is highly significant. It is larger than most values reported in the literature for lakes where denitrification (rather than algal uptake) is known to be the predominant removal process (Alexander et al ., 2002), but is similar to the estimates

from the SPARROW models for eastern North Carolina river basins (16–18 m per year; McMahon et al ., 2003).

The spatial structure of residuals, Model A

The nitrogen load and yield predicted by the modelclosely match the observed values as indicated by a standard error of the estimate (SEE) of 0

Ð 34, expressed in Figure 4. Observed and predicted (Model A) total nitrogen flux, as (A) load and (B) yield, at 321 monitoring sites in the southeastern United
States

log units, and by coefficients of determination ( r 2 ) of

0

Ð 96 and 0

Ð 68 for load and yield, respectively (Table I and Figure 4). The SEE is a measure of the size of the typical error of the model-estimated load or yield for a reach compared to the observed load for a reach; the

value of 0

Ð 34 is roughly equivalent to a mean percentage error of 34%. The coefficient of determination for nitrogen yield measures the fraction of variance in the observed nitrogen yield (expr essed in log units) that is accounted for by the model; thus Model A explains 68% of the variance in log-transf ormed yield observed for the set of monitoring sites. The unexplained variance can be attributed to errors in the input data sets, errors in modelspecification (structure or selection of explanatory variables), errors in the monitored yield estimates, or some

combination of these errors.

The spatial structure in residuals from Model A(Figure 5) reveals the tendency for Model A to overpredict in some areas (North and South Carolina) and

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Figure 5. Model residuals (observed yield minus predicted yield ) at 321 stream monitoring sites, from two model specifications

underpredict in other areas (Mississippi and Alabama).

The statistical test Moran’s I (Moran, 1950) provides a quantitative estimate of s patial autocorrelation: a z

score of 2

Ð 8 and associated p -value of 0

Ð 002 for Model A indicates spatial aut ocorrelation of residuals and, potentially, shortcomings in the specification of Model A (Table I). We suggest that the spatial biasis due to incorrect simulation of the DVF; specifically, that the calculated DVF is underestimated inareas of positive residuals (underprediction) and overestimated in areas of negative residuals (overprediction). To examine whether this spatial bias could be reduced and overall model fit improved, alternative combinations of all the local-scale landscape variables listed in the Methods section were incorporated to improve simulation of DVF. No combination of these landscape characteristics resulted in raising DVF in areas of underprediction and lowering DVF in areas of overprediction.

An alternative explanation for the observed spatial pattern of residuals is that re gional-scale spatial differences may exist in the effects of local-scale landscape characteristics, such as soil permeability, on nitrogen attenuation. Incorporating regional frameworks in the model as step functions allows for scaling the relation between DVF and landscape characteristics differently for different regions. The fact that the spatial pattern of residuals from Model A corresponds to mapped units in both the Level III ecoregion and HLR frameworks suggests that each of these regional frameworks provides information that may improve model simulation

of DVF.

Incorporating hydrologic landscape framework as a regional-scale variable to simulate spatial variation in landscape delivery ratio, Model B

Because the HLR framework explicitly hypothesizes that regions differ in terms of characteristic hydrologic response and primary flow path, a priori hypotheses can be constructed with respect to expected differences among HLRs in the quantitative relation between LDR and landscape characteristics. For example, where the primary flow path is shallow groundwater (HLRs 1, 2, 4, 7, and 16), a continuous decrease in LDR (and therefore in model estimates of DVF) with increasing soil depth is expected, but where primary flow path is overland (HLRs 6, 9, and 11), LDR may be unresponsive to changes in soil depth. The negative bias of Model A residuals (model overpredic tions) for calibration sites in three regions where flow paths are predominantly shallow groundwater (HLRs 2, 7, and 16; Table II) supports the hypothesis that the differing relation between LDR and local-scale landscape variables depends on primary flow path and that incorporating the HLR framework will improve simulation of the pattern of DVF and LDR beyond the Model A estimates.

A model specification including a set of variables describing the areal boundaries of all the HLRs in the SAGT model area (that is, HLRs 2, 4, 6, 7, 9, 11, and

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

16; HLR 1 was excluded to avoid colinearity), in addition to the three local-scale landscape variables included in Model A, was used to test this hypothesis. Each HLR

variable was defined as the fraction of the area of the catchment within the HLR. HLR 6 was under-represented by monitoring sites (four sites; Table II) to support calibration of a separate coefficient. HLR 6 therefore was combined with HLR 11 based on the expected similarity in DVF (primary hydrologic flow path for both HLRs is overland flow). When this combined variable failed as a statistically significant predictor, HLR 6 and 11 were combined with HLR 9. The resulting model specification (Model B) is listed in Table I.

Incorporating ecoregion framework as a regional-scale variable to simulate spatial variation in landscape delivery ratio, Model C

Because many factors other than those that influence nitrogen transport and attenua tion were used in ecoregionframework development, it is difficult to construct physically based hypotheses about th e influences of individual Level III ecoregions on the transport and attenuation ofnitrogen, or to predict expected improvements to estimation of delivery resulting from incorporating the ecoregion framework in model specification. The strong bias

of Model A residuals towards underprediction for calibration sites in the Interior Plateau, Southwestern Appalachians, and Ridge and Valley (Table II), however, gives an empirical basis for expecting that estimation of DVF will be improved by incorporating the ecoregion framework, despite the fact that the mechanisms associated with the improvement in model fit are not apparent.

A model specification that included a set of variables describing the areal boundaries of all the ecoregions in the SAGT model area (that is, ecoregions 45, 66, 67, 68, 71, and 75; ecoregion 65 was excluded to avoid colinearity), in addition to the three local-scale landscape variables included in Model A, was used to test the hypothesis that ecore gions can provide additional explanatory information and improve modelling of LDR. Each ecoregion variable was defined as the fraction of the area of the catchment within the ecoregion. Ecoregion 63 was under-represented by monitoring sites (four sites, Table II) in the calibration set and therefore was combined with ecoregion 65; this combination was based on geography (adjoining ecoregions) rather than on physically based hypotheses about the expected effect of the individual ecoregions on DVF. The resulting model specification (Model C) is listed in Table I.

Comparison of overall fit and coefficients—Models B and C compared with Model A

In addressing the first study question, ‘How does incorporation of regional landscape variables affect model error?’, we find that the yield r 2 and root mean square error (RMSE) for both Models B and C are improved, by similar amounts, compared with Model A (Table I). The change in these fit statistics is relatively small (less than 10%), however, and cannot be taken, by itself, to indicate substantial change in model performance by inclusion of the regional-scale variables. Change in the statistic

Moran’s I quantifying spatial autoc orrelation of residuals, however, is substantial and significant for Models B and C. Whereas Model A residuals exhibit significant spatial autocorrelation ( z -score for Moran’s I D 2

Ð 8 with p -value

0

Ð 002), Model B and C residuals are spatially independent ( z -score for Moran’s I < 0

Ð 7 with p -value > 0

Ð 2, Table I). Overprediction bias is reduced in HLRs 7 and16 (the Blue Ridge and Piedmont ecoregions), and underprediction bias is reduced in HLRs 4, 6, and 9 (the Ridge and Valley and Southwestern Appalachians ecoregions) (Figure 5B).

The coefficients ( ! D ) for the regional-scale variables effectively adjust the model-predicted relation between DVF and the other landscape variables (soil permeability, soil depth, and precipitation) upward or downward for catchments in each region. Four of the five coefficients for the HLR variables (Model B, Table I) were found to be statistically significant at the 5% significance level. Coefficient signs match closely with the interpretation that regions where the primary flow path is shallow groundwater deliver a smaller fraction of nitrogen mass

to the stream channel. The coefficients for variables representing HLRs 2, 7, and 16 (primary flow path is shallow groundwater) are negative, and the coefficient for the variable representing HLRs 6, 9, and 11 (primary flow path is overland flow) is positive. Only the coefficient for HLR 4 does not match the expected sign. The

modelled interaction between the HLR variables and the other landscape variables is illustrated, for mean annual precipitation, in Figure 3B.

Coefficients for the ecoregion-based variables (Model C) are not as significant, as a group, as those for the HLR-based variables (Model B). Only three of the six regional variables are found to be statistically significant at the 5% significance level (Table I). The coefficient for ecoregion 45 (Piedmont) is negative and significant at the 5% significance level; coefficients for ecoregion 68 and 75 (Southwestern Appalachians and Southern Coastal Plain) are positive and significant at 5% while coefficients for ecoregions 66, 67, and 71 (Blue Ridge, Ridge and Valley, and Interior Plateau) are not significant at 5%.

As previously noted, we have no theory-based reason to assume a particular relation between ecoregions and landscape nitrogen attenuation. Empirical results from Model C—spatially independe nt model residuals—only suggest that the broad-scale mosaic of landscape features in ecoregions 45, 68, and 75 (Piedmont, Southwestern Appalachians, and Southern Coastal Plain) produce biotic and abiotic characteristics in each of these regions that affect nitrogen attenuation.

Based on comparing the coefficient, ! R , for reservoir inverse hydraulic loading, the model-estimated rate of nitrogen removal in reservoirs is smaller in Models B and C (10

Ð 7 and 10

Ð 3 m per year, respectively) than in Model A (13

Ð 1 m per year), and the difference between Models A and C (2

Ð 8) equals the standard error for ! R in Model A (Table I). Based on comparing the coefficient, ! S , for

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

stream travel time, the model-estimated rate of nitrogen removal from small streams ( < 2

Ð 8 m 3 /s) is also smaller, by about half, in Models B and C (0

Ð 12 and 0

Ð 14 per day, respectively) compared with Model A (0

Ð 23 per day) and equals the standard error for the coefficient in Model A (Table I). This strongly suggests that some fraction of nitrogen attenuation ori ginally assigned in Model A to processes in small streams and in reservoirs has been attributed instead to landscape processes in Models B and C. The number of sites (34) representing the small stream class may not be sufficient to constrain model estimation

of the loss coefficient and stabilize the estimation of partitioning of loss in headwater watersheds between landscape and instream attenuation. This uncertainty may be resolved by using a finer-scale stream network, for example the 1 : 100 000 NHDPlus (see Methods section), as infrastructure for the SPARROW model, thus including monitoring sites on smaller stream sites in the

calibration set.

Understanding the usefulness of regional-scale variables for modelling nitrogen transport

Incorporating regional-scale landscape variables into the SPARROW model specification reduces the spatial correlation of model residuals, compared with a model that includes only local-scale landscape variables. This is true for both regional frameworks used here even though they differ in several important regards: the regional units are different in size; the regions were delineated based on different sets of processes; and regional boundaries were defined using different classification approaches.

The spatial structure in the residuals of Model A suggests that the model is mis-specified. The addition of regional variables to the model specification removes the spatial structure of the model residuals; however, it is unclear whether the regional variables serve as a proxy for missing or erroneous local-scale variables or whether there are, in fact, processes operating at a broader scale that affect nitrogen and need to be incorporated into the model. If the regional variables were serving as proxy forlocal-scale variables used to delineate the regional boundaries and missing from the SPARROW model specification, explicit incorporation o f these local-scale landscape characteristics in the model would be expected to reduce the significance of the regional variables as predictors of DVF. We tested this for HLRs by specifying DVF as a function of all the landscape ch aracteristics used to delineate the HLRs (percentage of flatland, soil permeability,hydraulic conductivity of the surficial aquifer, and precipitation excess; Wolock et al ., 2004), as well as the characteristics identified in Model A that explain variation in delivery (depth to bedrock). All except conductivity were significant predictors; however, the model fit was not improved compared with Model A (RMSE D 0

Ð 34) and the spatial structure in re siduals remained. Adding the HLR regional variable to this model with the expanded local-scale landscape variab les reduced the spatial structure in the residuals (change d autocorrelation statistic from significant to not significant). Unfortunately, this test cannot be performed for ecoregions because of their lack of explicit physical definitions.

If the regional variables were serving primarily to maskthe spurious influence of measurement error in the finerresolution, landscape characteristics used to delineate the regional boundaries, replacing the local-scale landscape

variables with rank-transforme d surrogate variables that are less affected by measurement error should have the same effect as incorporation of regional variables (i.e. reduced spatial structure). We tested this for HLR by modelling DVF as a function of rank-transformed values of each of the landscape characteristics used to

delineate the HLRs. The fit sta tistics and spatial structure of residuals from this model indicate no improvement compared to Model A. Testing was not possible for ecoregions because of the ir lack of explicitly physical

definitions.

The fact that the HLR regional-scale variables donot appear to be proxies for erroneous local-scale landscape information, or for variables included in the regional classification but missing from the SPARROWmodel specification, does not constitute proof that crossscale processes affect nitroge n attenuation. The empirical results from incorporating regional-scale variables do,however, suggest this possibility. This empirical finding is supported by the logically consistent interpretation of ! D estimates for the HLR variables: negative

! D estimates for HLRs where the primary flow path is shallow groundwater suggest that a lower fraction of nitrogen mass will be delivered to streams, and the reverse relation for regions where the primary flow path is overland flow. The empirical finding that broad-scaleclassifications of hydrologic response help explain differential rates of nitrogen attenuation, controlling for other local-scale landscape charact eristics, is consistent with the hypothesis presented in the Introduction section.

Effect of the regional landscape variables when using SPARROW to estimate delivery of nitrogen to streams

The different specifications of DVF in Models A, B, and C result in different spatial distributions of predicted DVF and predicted LDR for each of the land-applied sources (Figure 6). For all three models, areas with higher values of DVF (greater than 1

Ð 3, or the 75th percentile) occur in eastern T ennessee, northern Alabama, and southern Mississippi, which are areas with low values (relative to the average for the model area) of soil permeability and depth and relatively high values of mean annual precipitation. For Models B and C, however, these areas extend more widely and coincide with the extent of HLRs 4, 6, 9, and 11 (Model B) and ecoregions 67 and 68 (Model C). Catchments in these areas, such as the south Mouse Creek Basin (Figure 6), are estimated to transport a greater proportion of nitrogen to the stream, given equal nitrogen inputs, than the rest of the model area.

Areas with lower values of DVF (less than 0

Ð 75, or the 25th percentile) coincide with areas with relatively high values of soil permeability and depth and relatively low values of mean annual precipitation, primarily in

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

Figure 6. Delivery variation factor, repre senting variability in landscape delivery r atio, from three model specifications

central North Carolina, South Carolina, and Georgia. For Models B and C, these areas also extend more widely and coincide with the extent of HLRs 2, 7, and 16 (where the primary hydrologic pathway is shallow groundwater, Model B), and with ecoregion 45 (Model C). The Little Pee Dee River in South Carolina is typical of catchments in this group (Figure 6).

Spatial variability in DVF and LDR has implications for water-resource managers and planners in prioritizing areas for management acti ons because such variability can lead to uneven results in a uniform implementation of load reduction. For example, given an implementation plan for reducing inputs by 2

Ð 5 kg/ha for atmospheric deposition, for which estimated  ̨ (LDR avg ) is 0

Ð 50 kg/kg, the estimated decrease in load reaching streams in south Mouse Creek Basin, Tennessee (DVF D 1

Ð 5) would be

2

Ð 5 ð 0

Ð 50 ð 1

Ð 5 kg/ha, or 1

Ð 9 kg/ha, whereas the estimated decrease in load reaching streams in the Little Pee Dee River Basin, South Carolina (DVF D 0

Ð 68) would

be about half that amount

# 2

Ð 5 ð 0

Ð 50 ð 0

Ð 68 kg/ha, or

0

Ð 85 kg/ha. Stated as the inverse, achievement of a specified reduction in loading of land-applied sources to the stream would require about twice the reduction in source

mass for catchments in the Little Pee Dee River Basin as compared with the South Mouse Creek Basin.

The different spatial distributions of DVF estimated by Models A, B, and C result in markedly different predictions of instream load for certain areas in the

SAGT model area. Differences at the reach level are illustrated in Figure 7; the reach-level predictions of instream load are also provided in data files in the Supporting Information. The differences when compared

with Model A for 6300 of the 8028 reaches are within 23% (Model B) and 21% (Model C), but for almost 100 reaches, the differences exceed 35% (Model B) and 50% (Model C). Differences are most pronounced in HLR 16 or Blue Ridge, where Model B and C estimates are smaller than Model A estimates, and in HLR 4 and HLR 9, or Ridge and Valley, where Model B and C estimates are larger. While these results address the second research question in this paper, ‘How does inclusion of regional landscape variables affect mode l predictions?’, they raise an additional question: ‘Which of the three sets of model predictions is most accurate?’

As discussed in the previous section, model error statistics (specifically spatial a utocorrelation of residuals) favour Models B and C over Model A for prediction accuracy but do not differentiate between Models B and C. The predicted stream conditions (such as instream load) from Models B and C are not identical, however. Additional monitoring data for streams in areas where predictions from Models B and C diverge (HLRs 4, 9, and 16) are needed to differentiate performance. In the absence of empirical evidence, we suggest that the

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Figure 7. Instream nitrogen yield for each reach in the SAGT SPARROW model a rea (southeastern United States), 2002, estimated from three model specifications

greater interpretability of the Model B specification and calibration results with respect to processes affecting landscape attenuation of n itrogen favours the use of Model B predictions for model applications.

To examine differences at the basin scale, we examined variation among the three models in model predictions of several component s of the nitrogen deliverybudget—mass of nitrogen delivered to the stream network, mass removed by processing in streams and reservoirs, and mass delivered to the basin outlet—as well as in predictions of source share contribution, for the Tombigbee River Basin (Table III and Figure 8). The reach-level predictions fo r the full set of 8028 reaches are also provided (for Model B only) in data files in the Supporting Information.

The different model specifications of DVF result in varying estimates of nitrogen entering the stream network. Basin-average estimates for the Tombigbee

River Basin from Models B and C exceed the Model A estimate by 0

Ð 6 and 1

Ð 1 kg/ha, or about 12 and 20%, respectively. The larger Model B and C estimates in the Tombigbee River Basin result from the scaling up of DVF in HLRs 4, 6, 9, and 11 (Model B) and in the Southwestern Appalachians (Model C). Model B and C estimates

of mass removed in streams by aquatic processing are also larger (by 0

Ð 3 and 0

Ð 7 kg/ha) compared with Model A; this result is expected based on larger Model B and C

estimates of mass delivered to the stream and on similar (except for smallest stream size) Model B and C estimates of stream loss coefficients. In contrast, the estimates of

mass removed in reservoirs re main nearly constant across the three models, due to the smaller Model B and C

estimates of the reservoir loss coefficient.

Because the incremental change in Model B and C estimates of mass removed in streams almost equals and offsets the incremental change in Model B and C estimates of mass entering the stream network, the estimates of mass delivered to the basin outlet are nearly equal for all models, varying by less than 10%. The

increased values of LDR estimated for catchments in the Tombigbee River Basin in the Model B and C simulations do not translate, therefore, to proportional increases in

instream load estimates for all reaches. Estimates increase proportionally compared with Model A estimates for

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

reaches with mean streamflow less than 280 m 3 /s, but are almost unchanged for reaches on larger rivers.

The source shares (percentage of total mass) associated

with estimates for mass delivered to the stream network are similar (varying by less than 10% in most cases) to
those for mass delivered to basin outlet and are similar among the different models. This finding suggests that within the Tombigbee River Basin, the spatial distribution of each source variable [illustrated in the work by

Hoos et al . (2008)] is relatively unbiased with respect to the regional frameworks that modify simulated DVF in Models B and C. An exception is the spatial coincidence of higher input rate of atmospheric wet deposition (Hoos

et al ., 2008) with regions of increased (relative to Model A) Model C-estimated DVF (Figure 6). This spatial coincidence explains the modest increase in the share of atmospheric deposition in Model C (59%) compared with Model A (54%).

CONCLUSIONS

This study improves understanding of the relation between LDR for nitrogen and local- and regional-scale landscape characteristics, an d of the spatial arrangement of landscape delivery rates in the southeastern United States. Incorporating regional-scale landscape variables into a SPARROW model specification reduces the spatial correlation of model residuals, compared with a model that includes only local-scale landscape variables. We conclude from these empirical results that regional-scale differences in landscape processes result in differential effects of local-scale landscape characteristics, such as soil permeability, on nitrogen attenuation.

The models specified with the two alternate regional frameworks were not distinguishable based on modelerror statistics. In the absence of such empirical evidence, we suggest that the greater interpretability of the model specified with the HLR framework (ModelB) with respect to processes affecting landscape attenuation of nitrogen favours use of Model B predictions for model applications. The mapping units of the two regional frameworks differ in many parts of the model area primarily because of the finer-scaled detail in the HLR framework, but the similarity in error statistics for separate models derived from these frameworks highlights the need for a denser calibration network in the areas where predictions from the alternate models diverge(HLRs 4, 9, and 16), and points to the general conclusion that density of calibration data is a major factor in distinguishing among regionalization frameworks.

Both regional frameworks sca led the predicted relation

between LDR and local-scale soil and climate characteristics downwards in the Blue Ridge and Piedmont regions, and upwards in the Ridge and Valley, Interior Plateau, and Southwestern Appalachian regions. Estimates of LDR in the latter three regions were among the highest in

the model area based on soil and climate characteristicsalone (Model A) and were further increased in the models calibrated using the regional-scale variables (Models Table III. Total nitrogen budget and source shares for theTombigbee River Basin, 2002, estimated from three model specifications

Tombigbee River a

Model A Model B Model C

Mass delivered from

catchment to adjacent
stream channel 5

. 1 5

. 7 6

. 2

Contribution from individual

sources, in percent Point-source discharge of

wastewater 987

Atmospheric deposition 54 57 59 Impervious surface (urban land) 677

Commercial fertilizer applied to agricultural land 18 15 16

Animal manure from livestock production 13 13 11

Mass removed by

processing in streams
and reservoirs 1

. 2 1

. 5 1

. 9

Fraction removed by

processing 24 26 30

Streams 0

. 7 1

. 0 1

. 4
Reservoirs 0

. 5 0

. 5 0

. 5

Mass delivered to the basin

outlet (Mobile Bay) b 3

. 9 4

. 2 4

. 3

Contribution from individual

sources, in percent Point-source discharge of

wastewater 10 9 8

Atmospheric deposition 55 57 59 Impervious surface (urban land) 687

Commercial fertilizer applied to agricultural land 18 15 16

Animal manure from livestock production 11 11 10

Instream load

Average for reaches

< 28 m 3 /s 4

. 9 5

. 5 6

. 1

Average for reaches
28–280 m 3 /s 6

. 2 6

. 5 6

. 9

Average for reaches

> 280 m 3 /s 5

. 5 5

. 5 5

. 5

[Boldface values are reported as kilograms per hectare, other values are reported as percentage of total mass; estimates of mass delivered or removed represent the sum of reach-level estimates for the entire basin, divided by total basin area; spatial variation in landscape delivery
ratio is modelled in Model A as a function of three local-scale soil and climate variables (soil permeability, depth to bedrock, and mean annual precipitation), in Model B as a function of the local-scale variables as well as nominal variables representing hydrologic landscape regions (Wolock
et al ., 2004), and in Model C as a function of the local-scale variables as well as nominal variables represe nting Level III ecoregion (Omernik, 1987); m 3 /s, cubic meter per second.

a The Tombigbee River Basin budget includes results for all of hydrologic subregion 0316; this includes the area draining to the mouth of the Tombigbee River (confluence with the Alabama River) as well as the area contributing directly to the Mobile River and Mobile Bay.

b Estimated without adjusting predicted values with observed values, to allow for mass comparisons; therefore, these estimates do not exactly match model estimates of loading to Mobile Bay.

B and C). The modelled high LDR in these regions, however, are subject to greater uncertainty because of

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypSPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Figure 8. Total nitrogen budget for the Tombigbee River Basin, 2002, estimated from three model specifications

their sparser representation in the calibration set. The addition of nutrient monitoring sites in these regions may improve understanding of regional controls on nitrogen
transport rates.

Loss coefficients for small streams estimated by all three models are similar to values reported by other

SPARROW model studies for small streams but are substantially larger than observations from experimental studies; we suggest the larger modelled values may be a result of the coarse-scale (1 : 500 000) hydrography usedto define flow-path length for all these models. Comparison of estimated stream-loss coefficients among the three models points to the uncertainty in partitioning nitrogen losses between landscape and instream attenuation. This uncertainty may be resolved by using a finer-resolution

stream network that allows inclusion in the calibration set of monitoring sites on smaller streams.

ACKNOWLEDGMENTS

We are grateful to David Wolock, Gregory Schwarz, Jerad Bales, and William Wolfe of the US Geological Survey and two anonymous reviewers for their comments and suggestions on the manuscript. This work was supported by the National Water-Quality Assessment (NAWQA) Program of the US Geological Survey.

REFERENCES

Alexander RB, Smith RA, Schwarz GE. 2000. Effect of stream channel

size on the delivery of nitrogen to the Gulf of Mexico. Nature 403 : 758– 761. Alexander RB, Elliott AH, Shankar U, McBride GB. 2002. Estimating

the sources and transport of nutrients in the Waikato River Basin, New
Zealand. Water Resources Research 38 : 1268– 1290. Alexander RB, Smith RA. 2006. Trends in the nutrient enrichment of

U.S. rivers during the late 20th century and their relation to changes in probable stream trophic conditions. Limnology and Oceanography
51 : 639– 654.Alexander RB, Smith RA, Schwarz GE, Boyer EW, Nolan JV, Brakebill JW. 2008. Differences in phosphorus and nitrogen delivery to the Gulf of Mexico from the Mississippi River Basin. Environmental Science and Technology 42 (3): 822–830. Arnold JG, Fohrer N. 2005. SWAT2000: Current capabilities and research

opportunities in applied watershed modeling. Hydrological Processes
19 (3): 563– 572. Bailey RG. 1988. Problems with using overlay mapping for planning and

their implications for geographic information systems. Environmental Management 12 : 11–17. Bailey RG. 2004. Ecoregions and Subregions of the United States, Puerto

Rico, and the U.S. Virgin Islands (Bailey) , accessed August 2008 at http://national.atlas.gov/mld/ecoregp.html. Booth MS, Campbell C. 2007. Spring nitrate flux in the Mississippi

River Basin: a landscape model with conservation applications. Environmental Science and Technology 41 (15): 5410– 5418.Boyer EW, Alexander RB, Parton WJ, Li C, Butterbach-Bahl, K, Donner SD, Skaggs RW, Del Grosso SJ. 2006. Modeling denitrification in terrestrial and aquatic ecosystems at regional scales. Ecological Applications 16 : 2123– 2142. Bricker SB, Clement CG, Pirhalla DE, Orlando SP, Farrow DRG. 1999.

National estuarine eutrophication assessment; effects of nutrient
enrichment in the Nation’s estuaries , National Oceanic and Atmospheric Administration, National Ocean Service, Special Projects Office and the National Centers for Coastal Ocean Science, Silver Spring, Maryland, 71. Clement JC, Pinay G, Marmonier P. 2002. Seasonal dynamics of

denitrification along topohydrosequences in three different riparian
wetlands. Journal of Environmental Quality 31 : 1025– 1037. Cox T. 2007. As ethanol demand grows, so does ‘dead zone’ in Gulf

of Mexico, International Herald Tribune, July 23, 2007, accessed January 2, 2008, at http://www.iht.com /articles/2007/07/22/bloomberg/ bxdead.php. Destouni G, Lindgren GA, and Gren IM. 2006. Effects of inland

nitrogen transport and attenuation modeling on coastal nitrogen load
abatement. Environmental Science and Technology (Policy Analysis)
40 : 6208– 6214. De Wit MJM. 2001. Nutrient fluxes at the river basin scale, I: the PolFlow

model. Hydrological Processes 15 : 743– 759. Duff JH, Tesoriero AJ, Richardson WB, Straus EA, Munn MD. 2008.

Whole-stream response to nitrate loading in three streams draining agricultural landscapes. Journal of Environmental Quality 37 : 1133– 1144. Florinsky IV, McMahon S, Burton DL. 2004. Topographic control of soil

microbial activity: a case study of denitrifiers. Geoderma 119 : 33– 53. Gibson CA, Meyer JL. 2007. Nutrient uptake in a large urban river.

Journal of the American Water Resources Association 43 (3): 576–587.

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hypA. B. HOOS AND G. MCMAHON

Hoos AB. 2005. Evaluation of the SPARROW model for estimating

transport of nitrogen and phosphorus in streams of the Interior Low Plateau Ecoregion, Tennessee, Kentucky, and Alabama, during 1992– 2002 [abs]. In Proceedings of the Tennessee Water Resources Symposium, April 2005 , Burns, Tennessee, 2A–47.Hoos AB, Terziotti SE, McMahon G, Savvas K, Tighe KC, AlkonsWolinsky R. 2008. Data to support statistical modeling of instream nutrient load based on watershed attributes, southeastern United States,
2002 , U.S. Geological Survey Open-File Report 2008-1163, 50. Howarth RW, Billen G, Swaney D, Townsend A, Jaworski N, Lajitha K,

Downing JA, Elmgran R, Caraco N, Jordan T, Berendse F, Freney J, Kudeyarov V, Murdoch P, Zhao-Liang Z. 1996. Regional nitrogen budgets and Riverine N and fluxes for drainages to the North Atlantic
Ocean—natural and human influences. Biogeochemistry 35 : 75–139.
Jackson HC. 2007. Ethanol craze worsens Gulf ‘dead zone’ , The News

and Observer, December 18, 2007, accessed January 2, 2008, at http://www.newsobserver.com/prin t/tuesday/front/story/835729.html. Lamon EC, Qian SS. 2008. Regional scale stressor-response models

in aquatic ecosystems. Journal of the American Water Resources

Association 44 (3): 771–781. Lowrance RR, Leonard RA, Asmussen LE, Todd RL. 1985. Nutrient

budgets for agricultural watersheds in the southeastern coastal plain. Ecology 66 (1): 287–296. McMahon G, Gregonis S, Waltman S, Omernik JM, Thorson T, Freeouf J,

Rorick A, Keyes J. 2001. Developing a spatial framework of common ecological regions for the conterminous United States. Environmental Management 28 : 293– 316. McMahon Gerard, Alexander RB, Qian S. 2003. Support of TMDL

programs using spatially referenced regression models. ASCE Journal of Water Resources Planning and Management 129 : 315–329.McMahon G, Wiken EB, Gauthier DA. 2004. Toward a scientifically rigorous basis for developing mapped ecological regions. Environmental Management 34 (Suppl 1): S111–S124.
Miller JA. 1990. Ground water atlas of the United States; Alabama,

Florida, Georgia, and South Carolina . U.S. Geological Survey Hydrologic Atlas 730–73G. Mitsch WJ, Day JW, Gilliam JW, Groffman PM, Hey DL, Randall GW,

Naiming W. 2001. Reducing nitrogen loading to the Gulf of Mexico from the Mississippi River Basin: Strategies to counter a persistent ecological problem. Bioscience 51 (5): 373– 388. Moore RB, Johnston CM, Robinson KW, Deacon JR. 2004. Estimation

of total nitrogen and phosphorus in New England streams using spatially referenced regression models , U.S. Geological Survey Scientific Investigations Report 2004– 5012, 42. Moran, PAP. 1950. Notes on continuous stochastic phenomena.

Biometrika 37 : 17–21. Mulholland PJ, Tank JL, Webster JR, Bowden WB, Dodds WK,

Gregory SV, Grimm NB, Hamilton SK, Johnson SL, Marti E, McDowell WH, Merriam JL, Meyer J, Peterson BJ, Valett HM, Wolheim WM. 2002. Can uptake length in streams be determined by nutrient addition experiments? Results from an interbiome comparison study. Journal of the North American Benthological Society 21 : 544– 560. Nolan JV, Brakebill JW, Alexander RB, Schwarz GE. 2002.

ERF1 2—Enhanced River Reach File 20, accessed August 2007 at http://water.usgs.gov/GIS/me tadata/usgswrd/XML/erf1 2.xml. Omernik JM. 1987. Ecoregions of the conterminous United States,

Map (scale 1 : 7500,000). Annals of the Association of American Geographers 77 (1): 118–125. Omernik JM. 1995. Ecoregions: a spatial framework for environmental

management. In Biological assessment and criteria: tools for water resource planning and decision-making , Davis, WS and Simon, TP (eds). Lewis Publishers: Boca Raton; 49–62. Omernik JM. 2004. Perspectives on the nature and definition of ecological

regions. Environmental Management 34 (Suppl 1): S27–S38. Peters DPC, Groffman PM, Nadelhoffer KJ, Grimm NB, Collins SL,

Michener WK, Huston MA. 2008. Living in an increasingly connected
world: a framework for continental-scale environmental science. Frontiers in Ecology and the Environment 6 (5): 229–237.
Preston SD Brakebill JW. 1999. Application of spatially referenced

regression modeling for the evaluation of total nitrogen loading in the Chesapeake Bay watershed , U.S. Geological Survey Water-Resources Investigations Report 99–4054, 12. Rabelais NN, Turner RE, Scavia, D. 2002. Beyond science into policy:

Gulf of Mexico hypoxia and the Mississippi River. BioScience 52 (2):
129–141. Rumenik, RP. 1988. Runoff to streams in Florida: Florida Geological

Survey Map Series 122 , Tallahassee, Florida, Florida Geological Survey, 1 sheet. Scavia D, and Donnelly KA. 2007. Reassessing Hypoxia Forecasts for

the Gulf of Mexico. Environmental Science and Technology 41 (23):
8111– 8117. Schwarz GE, Hoos AB, Alexander RB, Smith RA. 2006. The SPARROW

surface water-quality model: Theo ry, application, and user documentation . Book 6, Chapter B3, Techniques and Methods Report, U.S. Geological Survey. Smith RA, Schwarz GE, Alexander RB. 1997. Regional interpretation

of water-quality monitoring data. Water Resources Research 33 (12):
2781– 2798. Terziotti SE, Hoos AB. 2008. Comparison of two stream networks

as framework for the SPARROW model for the southeastern U.S. Proceedings of the NHDPlus National Workshop, February 26–28,
2008 , Denver, Colorado. U.S. Environmental Protection Agency. 2007. Level IV Ecoregions,

accessed January 23, 2008, at http://www.epa.gov/wed/pages/ ecoregions/level iv.htm. U.S. Environmental Protection Agency and U.S. Geological Survey.

2006. NHDPlus user’s Guide, accessed November 25, 2007, at http://www.horizon-systems.com/nhdplus/. U.S. Geological Survey. 2001. National Land Cover Database 2001

(NLCD 2001), U.S. Geological Survey, accessed August 2007 at http://www.mrlc.gov/mrlc2k nlcd map.asp/. Wolock DM, Winter TC, McMahon G. 2004. Delineation and evaluation

of hydrologic-landscape regions in the United States using geographic information system tools and multivariate statistical analysis. Environmental Management 34 (Suppl 1): S71–S88.

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009) DOI: 10.1002/hyp