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. Hoos1* ,‡ and Gerard McMahon2,‡

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

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

Abstract:

Understanding how nitrogen transport across the landscape varies with landscape characteristics is important for developing

sound nitrogen management policies. We used a spatially referenced regression analysis (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) among watersheds in the southeastern United States—higher in the western part

(Tennessee, Alabama, and Mississippi) than in the eastern part, and the average value for the region is lower compared to other

parts of the nation. When we model landscape delivery ratio as a continuous function of local-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 predicted load). The spatial pattern of modelled landscape delivery ratio and the spatial pattern

of residuals coincide spatially with Level III ecoregions and also with hydrologic landscape 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, and suggests that cross-scale processes affect nitrogen attenuation on the landscape.

The model-fitted coefficient values are logically consistent with the hypothesis that broad-scale classifications of hydrologic

response help to explain differential rates of nitrogen attenuation, controlling for local-scale 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; spatially referenced regression; watershed 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 programs 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 watershed, balanced by mass

of nitrogen exported by crop harvesting and also by the

interaction of environmental characteristics that influence

transport and attenuation of nitrogen during overland and

subsurface 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 analysis

is 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 these

processes 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 regional

patterns 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 than

100 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). The

map 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 studies

have discovered numerous examples of spatial variation in the set of controlling factors. For example,

Florinsky et al. (2004) demonstrated that topographic

properties, 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 streams 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 as a continuous function of

landscape characteristics because this approach cannot

capture spatial variation in the set of controlling factors

and their quantitative relation 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 from

other regions in terms of either a single type of process, such as a characteristic hydrologic response, or a

more 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 economic, climatic, and physiographic processes. These interactions tend to organize

the landscape into a distinct spatial mosaic of recognizable landscape features (Omernik, 1987, 2004; Bailey,

1988, 2004).

Wolock et al. (2004) developed hydrologic landscape

regions (HLRs) to differentiate areas with characteristic hydrologic response and primary flow paths. The US

Environmental Protection Agency (USEPA) developed

ecological regions that reflect a hypothesis that ecological processes, such as denitrification, may vary from one

region of the country to another due to broad geographic

variation 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 processes

that 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 regression 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 equation 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/hyp
SPATIAL 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 nitrogen 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 variables, 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 is

distinct 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 regression

equation 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 spatially referenced datasets

representing stream-channel networks, nitrogen sources,

physical watershed characteristics, and observations of

instream nitrogen load. The input 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 annual 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 ecoregions within the southeastern United States, 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/hyp
A. 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:

Linstreami D Lcatchmenti C Lupstreami !1"

where Linstreami D instream load at downstream node

of reach i; Lcatchmenti D load originating within the

catchment for reach i and delivered to the downstream

node of reach i; and Lupstreami 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

(Lcatchmenti) is determined by:

L catchmenti D

!NS

nD1

Sn,i ̨nDn

"

ZD

i ; #D

# A"

ZS

i , ZRi ; #S, #R

#

!2"

where n, Ns D source index where Ns is the total number

of individual sources;

Sn 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 (LDRavg), 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 LDRavg for the source.

Dn!Ð" 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 characteristics that influence nitrogen

attenuation. The DVF for catchment i is multiplied by the

LDRavg for source n (that is, by  ̨n) to calculate LDRi,n.

ZD 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.

ZS and ZR D vectors of measured stream and reservoir

variables, respectively (examples 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 (Lcatchmenti"1,

Lcatchmenti"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 Lcatchmenti 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 segments 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 : 500

000-scale Reach File 1 (RF1) with enhancements to support national and regional-scale 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 estimation 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/hyp
SPATIAL 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 on

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 km2) 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 agricultural 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 were

defined by multivariate-analysis to describe the combined effect of landscape and climate characteristics (soil

texture, aquifer permeability, land surface elevation and

slope, precipitation, and potential evaporation) on hydrologic response. The HLR is a geographically independent

framework 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 by

multivariate 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 relatively 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/hyp
A. B. HOOS AND G. MCMAHON

derivation of both frameworks (Figure 2). Differences

may 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 in

four 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 three

models are listed in Table I, along with estimated coefficients and calibration error statistics. To test the effect

of regional landscape variables 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 specification

for Model A: atmospheric deposition, commercial fertilizer applied to agricultural land, animal manure from

livestock production, impervious surface area, and pointsource discharge of wastewater. The model-estimated

coefficients  ̨ (Table I) represent, for the land-applied

sources, the LDRavg 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 kilogram

of 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 LDRavg) for fertilizer and animal manure compared to atmospheric deposition correspond to modelestimated, smaller fractions of the measured input from

these sources entering land-phase transport. The modelfitted values are scale-specific 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 km2 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-based

SPARROW 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 in

quantifying contributions from ammonia emissions associated with other sources in the model (animal manure

and commercial fertilizer). In contrast, we use total inorganic nitrogen to characterize atmospheric nitrogen deposition in the SAGT model because estimates of wet deposition of total inorganic nitrogen, animal manure, and

fertilizer are not correlated significantly (coefficients of

determination r2 < 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 LDRavg

for the southeastern United States compared with the

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009)

DOI: 10.1002/hyp
SPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Table I. Calibration results for three nitrogen SPARROW models (Models A, B, and C) for the southeastern U.S

Predictor variable Model-fitted coefficient

Model A Model B Model C

Description and unit Unit Value Standard error p Value Standard error p Value Standard error pS (Eq. 2), source input variables a (Eq. 2)a

Nitrogen mass in permitted wastewater discharge,2002; kg/yr

kg/kg 0⋅80 0⋅10

<0⋅005 0⋅79 0⋅09

<0⋅005 0⋅79 0⋅10

<0⋅005

Wet deposition of inorganic nitrogen (ammonia andnitrate), detrended to 2002; kg/yr

kg/kg 0⋅50 0⋅05

<0⋅005 0⋅50 0⋅05

<0⋅005 0⋅53 0⋅05

<0⋅005

Area of impervious surfaces, 2001; kg/yr kg/km2 1990 666 0⋅01 2470 649

<0⋅005 2230 399

<0⋅005

Nitrogen mass in commercial fertilizer applied toagricultural land, 2002; kg/yr

kg/kg 0⋅13 0⋅03

<0⋅005 0⋅11 0⋅02

<0⋅005 0⋅11 0⋅02

<0⋅005

Nitrogen mass in manure from livestock production,2002; kg/yr

kg/kg 0⋅05 0⋅02 0⋅02 0⋅05 0⋅02 0⋅01 0⋅04 0⋅01

<0⋅005

ZD (Eq. 2), physical landscape variables qD (Eq. 2 and 3)a

Ln 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⋅005

Ln of depth to bedrock; ln of cm −0⋅50 0⋅16

<0⋅005 −0⋅284 0⋅169 0⋅035 −0⋅094 (notretained) 0⋅710

Ln of mean annual precipitation; ln of mm 1⋅4 0⋅32

<0⋅005 1⋅2 0⋅32

<0⋅005 1⋅0 0⋅27

<0⋅005

Fraction of catchment in HLR 2; dimensionless −0⋅29 0⋅15 0⋅02

Fraction in HLR 7 −0⋅31 0⋅14 0⋅01

Fraction HLR 16b −0⋅14 0⋅14 0⋅11

Fraction in HLR 4 0⋅26 0⋅12 0⋅02

Fraction in HLRs 6, 9, or 11 0⋅26 0⋅11 0⋅02

Fraction in Ecoregion 45 −0⋅38 0⋅09

<0⋅005

Fraction in Ecoregion 66b −0⋅06 0⋅11 0⋅56

Fraction in Ecoregion 67b 0⋅17 0⋅12 0⋅13

Fraction in Ecoregion 68 0⋅50 0⋅15

<0⋅005

Fraction in Ecoregion 71b 0⋅02 0⋅16 0⋅89

Fraction in Ecoregion 75 0⋅52 0⋅11

<0⋅005

ZS (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; dayc 0⋅14e 0⋅05

<0⋅005 0⋅12e 0⋅04

<0⋅005

Time of travel in reach segments where meanQ>2.8 and

<28 m3/s; dayc

per day 0⋅13 0⋅05 0⋅005

Published in 2009 by John Wiley & Sons, Ltd. Hydrol. Process. (2009)

DOI: 10.1002/hyp
A. B. HOOS AND G. MCMAHON

Table I. (Continued)

Predictor variable Model-fitted coefficient

Model A Model B Model C

Description and unit Unit Value Standard error p Value Standard error p Value Standard error p

Time of travel in reach segments where meanQ>28 m3/s; dayc,d

per day 0⋅00 0⋅02 0⋅43 0⋅014 0⋅02 0⋅26 0⋅02 0⋅03 0⋅52

ZR (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⋅005

Model A Model B Model C

Standard Error of the Estimate (SEE), computed as

root mean square error (RMSE)

0⋅34 0⋅32 0⋅32

Coefficient of determination (R2) of load estimate 0⋅96 0⋅96 0⋅96

Coefficient of determination (R2) of yield estimate 0⋅68 0⋅72 0⋅73

Moran’s I of log residuals (z-score and significance of

departure from the null hypothesis

2⋅8 0 ⋅002 0⋅67 0⋅25 0⋅68 0⋅25

of spatially independent residuals)

Variables selected for these model specifications were significant at alpha

D 0.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. a Model 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 between

the spatial variation of landscape delivery ratio (spatial variation is modeled as the delivery variation factor, DVF) and the landscape variables. b The variables representing fraction of catchment in HLR 16 and in Ecoregions 66, 67, and 71 are retained despite lack of significance at alpha

D 0.05 in order to represent the complete regional classification

framework. HLR 1, and Ecoregions 63 and 65, have been excluded from Models B and C (respectively) to avoid colinearity. c Evaluated only for reaches in the indicated size class, otherwise set to 0.

d The variable time of travel in reach segments with mean discharge >28 m3/s is retained despite lack of significance at alpha

D 0.05, in order to complete the description of nitrogen removal across the full range

of stream sizes.

e When 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 alpha

D 0.05. The

two 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/hyp
SPATIAL 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 leading 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 LDRavg

for fertilizer, 0Ð10 kg/kg, estimated from nutrient budget

studies in small agricultural catchments in the southeastern coastal plain (Lowrance et al., 1985), although

the budget-based estimates are not directly comparable because of scale differences (average catchment size

about 1 km2

).

The estimate of  ̨ for impervious surface area,

1990 kg/km2

, associates the mass of nitrogen delivered

to a stream channel with unit area of impervious surface. In this application, impervious surface area serves

as a surrogate for many diffuse sources of nitrogen

in urban areas, such as vehicle emissions, lawn fertilizer, and onsite sewage-disposal systems. The estimate

of 1990 kg/km2 for the SAGT SPARROW model cannot be compared with values estimated for other SPARROW models because estimates of impervious surface

area 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/km2

, respectively, are about one-third the

magnitude of the SAGT SPARROW estimate of  ̨ for

impervious surface area (1990 kg/km2). 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 variable 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 attenuation cannot fully account

for point-source attenuation that may include higher rates

in 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 !"MD

mD1

ωnmZD

m i!Dm#

, "3#

where

MD D the number of landscape variables;

ZD

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, fertilizer applied to agricultural

land, animal manure, impervious surface area) with

all MD landscape variables and is set to 0 for the

interaction of point-source discharge of wastewater with

all MD 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 landscape

variables remain constant at their mean values). The relation of the individual landscape variables to the DVF

and 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 the

variable 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/hyp
A. B. HOOS AND G. MCMAHON

Figure 3. Variation of delivery variation factor and landscape delivery ratio for manure from livestock operations with mean annual precipitation and

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 Level

III 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 denitrification and sedimentation for

which 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!!ZS "S#,

where ZS 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 m3/s) streams, 0Ð13 per day of travel for

intermediate (mean annual flow 2Ð8–28 m3/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 m3/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/hyp
SPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Table II. Distribution of Model A estimated delivery variation factor and residuals for calibration sites, with respect to level III ecoregions and hydrologic landscape regions

Level III ecoregion Hydrologic landscape region with similar extent to corresponding ecoregionEcoregion numberand name

Number (density) of

sites with upstream

area pre-dominantly

in ecoregion

Mean value of

Model A

estimated delivery

variation factor

Tukey

class

Mean value

of Model A

residual

HLR name and

description

Primary

hydrologic

flow path

Number (density)of

sites with upstream

area pre-dominantly

in HLR

Mean value of

Model A

estimated delivery

variation factor

Tukey

class

Mean value

of Model A

residual

66-Blue Ridge 22 (0.6) 1.2 B

!0.10 16 S 33 (0.6) 1.2 A

!0.10

(Mt PS IB)

45-Piedmont 125 (0.8) 0.93 C

!0.06 7 S 56 (0.8) 0.97 AB

!0.10

(Pn PS IB)

2 S,D 36 (0.3) 0.82 B

!0.12

(Pn PS PB)

65-Southeastern Plains 72 (0.2) 0.96 C 0.02 4

S,D 29 (0.2) 1.2 A 0.16 (Pn PS PB)

a 6 O 4 (0.1) 0.94 AB 0.12 (Pn IS IB)63-Middle Atlantic Coastal Plain 4 (0.1) 0.82 C 0.11

1 S,D 46 (0.1) 0.79 B 0.08 (Pn PS PB)

75-Southern Coastal Plain 36 (0.1) 0.68 C 0.11

71-Interior Plateau 6 (0.2) 1.3 B 0.17 11 O 9 (0.4) 0.92 AB 0.06 (Pt IS IB)68-Southwestern Appalachians 10 (0.3) 1.5 A 0.15 b67-Ridge and Valley 20 (0.4) 1.1 B 0.17 9 O,D 23 (0.3) 1.2 A 0.13 (Pt IS IB)Multiple 26 0.96

C

!0.07 Multiple Multiple 83 0.96 AB 0.01

ANOVA r2 0.45 0.06 0.43 0.08[Model residual, natural log of observed load minus natural log of model-predicted load, negative values indicate overprediction; Level III ecoregions described in Omernik, 1987; density of sites expressed as ratio

of 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 Wolock et al, 2004: Mt

D mountains, Pn

Dplains, Pt

D plateau, PS

D permeable soils, IS

D impermeable soils, PB

D permeable bedrock, IB

D impermeable bedrock; primary hydrologic pathway: S

D shallow ground water, O

D overland flow, D

D deep

ground water; delivery factor estimated based on soil permeability, depth to bedrock, 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 the

regional framework].

a Ecoregions 45 and 65 contain small areas of HLR 6.b Ecoregion 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/hyp
A. 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) but

larger than observations from experimental studies (typically less than 0Ð15 per day, Howarth et al., 1996). We

suggest 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 size

classes 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 m3/s)

class, 117 load observations (about half) represent sites

with 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 ZR "R#,

where

ZR 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 model

closely 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 (r2) 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 (expressed in log units) that is

accounted for by the model; thus Model A explains 68%

of the variance in log-transformed yield observed for the

set of monitoring sites. The unexplained variance can be

attributed to errors in the input data sets, errors in model

specification (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/hyp
SPATIAL 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 spatial autocorrelation: a zscore of 2Ð8 and associated p-value of 0Ð002 for

Model A indicates spatial autocorrelation of residuals

and, potentially, shortcomings in the specification of

Model A (Table I). We suggest that the spatial bias

is due to incorrect simulation of the DVF; specifically, that the calculated DVF is underestimated in

areas 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 regional-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 overpredictions) 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/hyp
A. 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 attenuation were used in ecoregion

framework development, it is difficult to construct physically based hypotheses about the influences of individual

Level III ecoregions on the transport and attenuation of

nitrogen, 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 ecoregions 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 r2 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 autocorrelation 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 and

16 (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 independent 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/hyp
SPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

stream travel time, the model-estimated rate of nitrogen

removal from small streams (<2Ð8 m3/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 originally 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 for

local-scale variables used to delineate the regional boundaries and missing from the SPARROW model specification, explicit incorporation of 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 characteristics 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 residuals remained. Adding the

HLR regional variable to this model with the expanded

local-scale landscape variables reduced the spatial structure in the residuals (changed 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 mask

the spurious influence of measurement error in the finerresolution, landscape characteristics used to delineate the

regional boundaries, replacing the local-scale landscape

variables with rank-transformed 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 statistics and spatial structure

of residuals from this model indicate no improvement

compared to Model A. Testing was not possible for

ecoregions because of their lack of explicitly physical

definitions.

The fact that the HLR regional-scale variables do

not appear to be proxies for erroneous local-scale landscape information, or for variables included in the

regional classification but missing from the SPARROW

model specification, does not constitute proof that crossscale processes affect nitrogen 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-scale

classifications of hydrologic response help explain differential rates of nitrogen attenuation, controlling for other

local-scale landscape characteristics, 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 Tennessee, 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/hyp
A. B. HOOS AND G. MCMAHON

Figure 6. Delivery variation factor, representing variability in landscape delivery ratio, 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 actions 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  ̨ (LDRavg) 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 model 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 autocorrelation 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/hyp
SPATIAL ANALYSIS OF NITROGEN LOADS TO STREAMS IN SOUTHEASTERN USA

Figure 7. Instream nitrogen yield for each reach in the SAGT SPARROW model area (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 nitrogen 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 components of the nitrogen delivery

budget—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 for 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 remain 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/hyp
A. B. HOOS AND G. MCMAHON

reaches with mean streamflow less than 280 m3/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, and 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 model

error statistics. In the absence of such empirical evidence, we suggest that the greater interpretability of

the model specified with the HLR framework (Model

B) 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 scaled 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 characteristics

alone (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 the

Tombigbee River Basin, 2002, estimated from three model specifications

Tombigbee Rivera

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 m3/s

4.9 5.5 6.1

Average for reaches

28–280 m3/s

6.2 6.5 6.9

Average for reaches

>280 m3/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 representing Level III ecoregion (Omernik,

1987); m3/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/hyp
SPATIAL 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 used

to 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/hyp
A. 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/print/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/metadata/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: Theory, 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