U.S. Department of the Interior

U.S. Geological Survey

Open-File Report 2009–1122

Update of Watershed Regressions for Pesticides (WARP) 

for Predicting Atrazine Concentration in Streams

Update of Watershed Regressions for 

Pesticides (WARP) for Predicting Atrazine 

Concentration in Streams

By Wesley W. Stone and Robert J. Gilliom

Open-File Report 2009–1122

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

KEN SALAZAR, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia 2009

For more information on the USGS—the Federal source for science about the Earth, its natural and living resources, 

natural hazards, and the environment, visit http://www.usgs.gov or call 1-888-ASK-USGS

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Stone, W.W.,and Gilliom, R.J., 2009, Update of watershed regressions for pesticides (WARP) for predicting atrazine 

concentration in streams: U.S. Geological Survey Open-File Report 2009–1122, 22 p.
iii

Foreword

The U.S. Geological Survey (USGS) is committed to providing the Nation with reliable scientific information that helps to enhance and protect the overall quality of life and that facilitates effective management of water, biological, energy, and mineral resources (http://www.usgs.gov/). Information on the 

Nation’s water resources is critical to ensuring long-term availability of water that is safe for drinking and 

recreation and is suitable for industry, irrigation, and fish and wildlife. Population growth and increasing 

demands for water make the availability of that water, measured in terms of quantity and quality, even 

more essential to the long-term sustainability of our communities and ecosystems.

The USGS implemented the National Water-Quality Assessment (NAWQA) Program in 1991 to support 

national, regional, State, and local information needs and decisions related to water-quality management and policy (http://water.usgs.gov/nawqa). The NAWQA Program is designed to answer: What is the 

quality of our Nation’s streams and groundwater? How are conditions changing over time? How do natural 

features and human activities affect the quality of streams and groundwater, and where are those effects 

most pronounced? By combining information on water chemistry, physical characteristics, stream habitat, 

and aquatic life, the NAWQA Program aims to provide science-based insights for current and emerging water issues and priorities. During 1991 to 2001, the NAWQA Program completed interdisciplinary 

assessments and established a baseline understanding of water-quality conditions in 51 of the Nation’s 

river basins and aquifers, referred to as Study Units (http://water.usgs.gov/nawqa/studyu.html). 

National and regional assessments are ongoing in the second decade (2001–2012) of the NAWQA Program as 42 of the 51 Study Units are selectively reassessed. These assessments extend the findings in the 

Study Units by determining water-quality status and trends at sites that have been consistently monitored 

for more than a decade, and filling critical gaps in characterizing the quality of surface water and groundwater. For example, increased emphasis has been placed on assessing the quality of source water and 

finished water associated with many of the Nation’s largest community water systems. During the second 

decade, NAWQA is addressing five national priority topics that build an understanding of how natural features and human activities affect water quality, and establish links between sources of contaminants, the 

transport of those contaminants through the hydrologic system, and the potential effects of contaminants 

on humans and aquatic ecosystems. Included are studies on the fate of agricultural chemicals, effects of 

urbanization on stream ecosystems, bioaccumulation of mercury in stream ecosystems, effects of nutrient 

enrichment on aquatic ecosystems, and transport of contaminants to public-supply wells. In addition, 

national syntheses of information on pesticides, volatile organic compounds (VOCs), nutrients, trace elements, and aquatic ecology are continuing. 

The USGS aims to disseminate credible, timely, and relevant science information to address practical and 

effective water-resource management and strategies that protect and restore water quality. We hope this 

NAWQA publication will provide you with insights and information to meet your needs, and will foster 

increased citizen awareness and involvement in the protection and restoration of our Nation’s waters. 

The USGS recognizes that a national assessment by a single program cannot address all water-resource 

issues of interest. External coordination at all levels is critical for cost-effective management, regulation, 

and conservation of our Nation’s water resources. The NAWQA Program, therefore, depends on advice 

and information from other agencies—Federal, State, regional, interstate, Tribal, and local—as well as 

nongovernmental organizations, industry, academia, and other stakeholder groups. Your assistance and 

suggestions are greatly appreciated.

Matthew C. Larsen

Associate Director for Water

v

Contents

Abstract...........................................................................................................................................................1

Introduction.....................................................................................................................................................1

Purpose and Scope ..............................................................................................................................2

Methods...........................................................................................................................................................2

Atrazine-Concentration Data Used for Model Development and Validation .............................2

Watershed Characteristics Used as Explanatory Variables..........................................................2

Estimation of Atrazine Use ........................................................................................................2

Other Watershed Characteristics ............................................................................................3

Statistical Analysis ...............................................................................................................................3

Regression Methods ...................................................................................................................3

Transformations of Response and Explanatory Variables ....................................................6

Selection of Explanatory Variables...........................................................................................6

 Analysis of Model Fit..................................................................................................................8

Estimation of Prediction Intervals.............................................................................................8

Atrazine Models.............................................................................................................................................8

Analysis of Significant Explanatory Variables .................................................................................8

Model Performance............................................................................................................................10

Model-Development Sites........................................................................................................10

Model-Validation Sites..............................................................................................................12

Uncertainty in Model Predictions ...........................................................................................12

Comparison to Previous WARP Models..........................................................................................18

Model Limitations................................................................................................................................19

Summary and Conclusions.........................................................................................................................21

References Cited..........................................................................................................................................21

Figures

1. Map showing location of sampling sites used for WARP model development 

(112 sites) and validation (25 sites). ...........................................................................................3

2. Boxplot explanation......................................................................................................................8

3–5. Diagrams showing:

 3. Atrazine-concentration statistics from concentrations observed at modeldevelopment sites in relation to values of the same statistics predicted by the 

WARP models..............................................................................................................................11

 4. Residual errors for atrazine-concentration statistics predicted by WARP 

models in relation to watershed area for the 112 model-development sites ...................13

 5. Residual errors for atrazine-concentration statistics predicted by the WARP 

models in relation to geographic region for the 112 model-development sites...............14

6. Map showing U.S. Department of Agriculture Farm Resource Regions 

(U.S. Department of Agriculture, 2000) and regions used for evaluation of the 

WARP models..............................................................................................................................15
vi

7–10. Diagrams showing:

7. Atrazine-concentration statistics observed at model-validation sites in 

relation to values of the same statistics predicted by the WARP models ...............16

8. Residual errors for atrazine-concentration statistics predicted by the 

updated WARP models and the previous WARP models using data from 25 

model-validation sites.......................................................................................................17

9. Potential prediction errors for concentration statistics predicted with the 

atrazine WARP models .....................................................................................................18

10. Residual errors for atrazine-concentration statistics predicted by the 

updated WARP models and the previous WARP models ...........................................20

Tables

1. Watershed characteristics considered as explanatory variables for WARP models ......4

 2. Example results from stepwise model development for selected WARP models.............7

3. Summary of statistics and coefficients for variables in the atrazine WARP models........9

4. Comparison of updated WARP models to previous WARP models. ..................................19

Conversion Factors

Multiply By To obtain

Length

centimeter (cm) 0.3937 inch (in.)

millimeter (mm) 0.03937 inch (in.)

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

Area

square kilometer (km2) 247.1 acre

square kilometer (km2) 0.3861 square mile (mi2)

Volume

liter (L) 0.2642 gallon (gal)

Mass

gram (g) 0.03527 ounce, avoirdupois (oz)

kilogram (kg) 2.205 pound avoirdupois (lb)
vii

Abbreviations and Acronyms

μg/L microgram per liter

L liter

ARP Acetochlor Registration Partnership

CI confidence interval

CIESIN Center for International Earth Science Information Network

DEM Digital Elevation Model

DMRK DMRKYNETIC, Inc.

K soil erodibility factor (K-factor)

pR2 pseudo R-squared (R-squared value for tobit regression)

NASQAN National Stream Quality Accounting Network

NASS National Agricultural Statistics Service

NAWQA National Water-Quality Assessment

NLCD National Land Cover data set

NRI National Resource Inventory

NOAA National Oceanic and Atmospheric Administration

PERDUN percent of streamflow due to Dunne overland flow

PI prediction interval

PMJN total precipitation during May and June of the year of sampling

R rainfall erosivity factor (R-factor)

R2 coefficient of multiple determination

STATSGO State Soil Geographic database

UI use intensity, annual agricultural atrazine use in the watershed

USGS U.S. Geological Survey

USLE Universal Soil Loss Equation

WA watershed area (square kilometers)

WARP Watershed Regressions for Pesticides

WQL Water Quality Laboratory
viii
Update of Watershed Regressions for Pesticides (WARP) 

for Predicting Atrazine Concentration in Streams 

By Wesley W. Stone and Robert J. Gilliom 

Abstract

Regression models for predicting atrazine concentrations in streams were updated by incorporating refined annual 

atrazine-use estimates and by adding an explanatory variable 

representing annual precipitation characteristics. The updated 

Watershed Regressions for Pesticides (WARP) models enable 

improved predictions of specific pesticide-concentration statistics for unmonitored streams. 

Separate WARP regression models were derived for 

selected percentiles (5th, 10th, 15th, 25th, 50th, 75th, 85th, 

90th and 95th), annual mean, annual maximum, and annual 

maximum moving-average (21-, 60-, and 90-day durations) 

concentration statistics. Development of the regression models 

involved the same model-development data, model-validation 

data, and regression methods as those used in the original 

development of WARP. The original WARP models were 

based on atrazine-use estimates from either 1992 or 1997. This 

update of the WARP models incorporates annual atrazine-use 

estimates. In addition, annual precipitation data were evaluated as potential explanatory variables. 

The updated WARP models include the same five explanatory variables and transformations that were used in the original WARP models, including the new annual atrazine-use data. 

The models also include a sixth explanatory variable, total 

precipitation during May and June of the year of sampling. 

The updated WARP models account for as much as 82 percent 

of the variability in the concentration statistics among the 112 

sites used for model development, whereas previous WARP 

models accounted for no more than 77 percent. Concentration 

statistics predicted by the 95th percentile, annual mean, annual 

maximum and annual maximum moving-average concentration models were within a factor of 10 of the observed 

concentration statistics for most of the model development and 

validation sites. 

Overall, performance of the models for the development 

and validation sites supports the application of the WARP 

models for predicting atrazine-concentration statistics in 

streams and provides a framework to interpret the predictions 

in terms of uncertainty. For streams where direct measurements of atrazine are lacking, the updated WARP model 

predictions can be used to characterize the probable values 

of atrazine-concentration statistics for comparison to specific 

water-quality benchmarks. 

Introduction

Pesticide concentrations in streams vary widely across 

the United States. Each pesticide has a unique pattern of 

occurrence because factors such as pesticide use, application 

practices and timing, climate, and watershed characteristics 

vary geographically (Gilliom and others, 2006). Because of 

the geographic and temporal complexity of pesticide occurrence and concentrations, adequate monitoring of pesticide 

concentrations in the numerous streams in the United States 

that are potentially affected by pesticides is prohibitively 

expensive—particularly at sampling frequencies high enough 

to reliably estimate the necessary concentration statistics for 

assessing potential adverse effects on aquatic ecosystems or 

humans. The lack of adequate data on pesticide concentrations 

in most streams creates a need for tools that can be used to 

predict pesticide concentrations for streams that are not sufficiently monitored. Although such estimates cannot replace 

direct measurements when reliability requirements are high, 

they are useful for the initial screening-level steps of risk 

assessment and for efficiently guiding intensive monitoring. 

Larson and others (2004), using watershed characteristics as explanatory variables, developed and applied regression models for predicting annual frequency distributions of 

atrazine concentrations in streams. These Watershed Regressions for Pesticides (WARP), which originated from the work 

of Larson and Gilliom (2001), have been expanded to predict 

annual maximum and annual maximum moving-average 

concentrations of atrazine for 21-, 60-, and 90-day durations 

(Stone and others, 2008). For streams with inadequate direct 

measurements of atrazine, the WARP model predictions can be 

used to characterize the probable values of atrazine-concentration statistics for comparison to specific water-quality benchmarks. For example, during the initial screening-level steps 

of risk assessment, water-quality benchmarks for evaluating 

potential concerns for human health or aquatic life are compared to estimated environmental concentrations, which for 

atrazine are represented by the various concentration statistics. 
2 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

WARP is based on empirical relations between pesticide concentrations observed at monitoring sites and selected 

watershed characteristics available nationally—such as 

pesticide use—and soil and hydrologic characteristics. The 

development of WARP by Larson and others (2004) was 

based on atrazine-use estimates that were derived for Census 

of Agriculture reporting years 1992 and 1997. Subsequently, 

annual pesticide-use estimates became available for the entire 

model-development period (1992–2001). Estimated pesticide 

use in a watershed is the most influential variable in the WARP 

models for predicting pesticide concentrations in a stream 

(Larson and others, 2004). Thus, incorporation of the annual 

pesticide-use estimates into the WARP models was expected 

to enhance model performance. 

Purpose and Scope

This report describes updated WARP models for predicting atrazine concentration in streams. Improvements over 

previous work (Larson and others, 2004; Stone and others, 2008) include the incorporation of annual atrazine-use 

estimates on crops and the addition of annual precipitation 

data as potential explanatory variables. The WARP models 

are based on the most complete year of observations for each 

of the 112 sampling sites. The updated WARP models predict 

atrazine concentrations for selected percentiles of the annual 

frequency distribution (5th, 10th, 15th, 25th, 50th, 75th, 85th, 

90th, and 95th), as well as the annual mean, annual maximum, 

and annual maximum moving averages (21-, 60-, and 90-day 

durations).

Methods

Methods used in this study included obtaining atrazineconcentration data, estimating annual agricultural atrazine 

use, estimating watershed characteristics, and developing and 

evaluating regression models to predict atrazine-concentration 

statistics.

Atrazine-Concentration Data Used for Model 

Development and Validation 

Atrazine-concentration data from the 112 sampling 

sites used by Larson and others (2004) and Stone and others 

(2008) were used in this study. The sampling-site selection 

is discussed in detail in Larson and others (2004), and site 

locations are shown in figure 1. Data from these 112 sites were 

collected for the U.S. Geological Survey (USGS) National 

Water-Quality Assessment (NAWQA) and National Stream 

Quality Accounting Network (NASQAN) programs from 1992 

through 2001. Sample collection for the NAWQA program 

followed procedures described by Shelton (1994); sample collection for the NASQAN Program followed protocols of the 

U.S. Geological Survey (1997 to present).

Model-validation data collected from 25 sampling sites 

by the Water Quality Laboratory (WQL) of Heidelberg College in Tiffin, Ohio, the Acetochlor Registration Partnership 

(ARP), and the USGS, previously used by Larson and others 

(2004) and Stone and others (2008), were also used in the 

present study. Information about the WQL and ARP monitoring programs is provided in Richards and Baker (1993), Richards and others (1996), and Acetochlor Registration Partnership (2009). Locations of model-validation sampling sites for 

WQL and USGS are shown in figure 1. ARP sampling sites are 

not shown in figure 1 because these sites are associated with 

drinking-water sources, which are confidential. Larson and 

others (2004) used 26 model-validation sites, but this study 

used 25 because 1 site was sampled during a year for which 

annual atrazine-use data were not available. The concentration 

statistics (annual percentiles, mean, maximum, and maximum 

moving averages) used in model development for the present 

study are the same as those used by Larson and others (2004) 

and Stone and others (2008). 

Watershed Characteristics Used as Explanatory 

Variables

WARP models are an empirical regression approach to 

modeling that evaluates a large number of potential explanatory variables that could affect or indicate an influence on pesticide transport and runoff. Statistical procedures are used in 

the WARP model development to select explanatory variables 

on the basis of their ability to reproduce observed pesticide 

concentrations.

Estimation of Atrazine Use 

Atrazine use, as expressed by the fourth-root of use intensity (amount of atrazine used in a watershed divided by the 

watershed area), was the most important explanatory variable 

for predicting atrazine concentrations in the original WARP 

models (Larson and others, 2004). Atrazine-use data used in 

these original models were temporally static. Specifically, the 

atrazine use for a watershed was represented by a single value 

for either the 1991–94 period or the 1994–98 period (Larson 

and others, 2004) because actual annual atrazine-use estimates 

were not available. This update of the WARP models uses 

annual atrazine-use estimates that were made possible through 

obtaining proprietary annual use data from DMRKYNETEC, 

Inc. (DMRK). 

The process of computing annual atrazine use in a 

watershed, based on county-level use estimates, followed 

the same methods as described in Larson and others (2004). 

Areal weights of specific types of agricultural land for each 

county intersected by the watershed were computed. The areal 

weights were multiplied by county-level atrazine use on each 

individual agricultural land type to determine the countylevel atrazine use within the watershed for the individual land 

types. Total atrazine use in each watershed was estimated by 
Methods 3

summing the apportioned county-level atrazine-use estimates 

in the watershed. 

County-level estimates of atrazine use were obtained by 

integrating annual crop-reporting-district atrazine-use data on 

individual crops obtained from DMRK with annual countylevel crop acreage information from public sources (Gail Thelin, U.S. Geological Survey, written commun., 2009). Annual 

county-level crop acreage data were obtained from the Census 

of Agriculture for 1992, 1997, and 2002 and from National 

Agricultural Statistics Service (NASS) annual reports for years 

not reported by the Census of Agriculture.

The estimates of atrazine use derived by the methods 

described above represent agricultural use. There are no 

nationally available county-level estimates for nonagricultural 

use of atrazine or other pesticides. The primary use of atrazine 

is application to agricultural crops (Kiely and others, 2004); 

therefore, the lack of information on nonagricultural use likely 

has a minimal effect on the regression analysis.

Other Watershed Characteristics 

The many watershed characteristics evaluated as potential 

explanatory variables in the WARP model development are 

listed in table 1. Potential explanatory variables representing 

land use and population, agricultural management practices, 

soil properties, physical watershed characteristics, weather 

and climate characteristics, and hydrological properties were 

considered. These watershed characteristics are the same used 

by Larson and others (2004), with the exception of the annual 

precipitation characteristics and the improved annual atrazineuse intensity, which were added for this update of the WARP 

models. All variables considered during model development 

are available for the conterminous United States.

Statistical Analysis

A regression approach was used in the development of 

the WARP models. Statistical methods used to support the 

regression approach include transformations of response and 

explanatory variables, selection of explanatory variables, 

analysis of model fit, and estimation of prediction intervals. 

Regression Methods

Some of the atrazine-concentration statistics computed 

for sites included in the regression analysis were censored 

at the long-term reporting level. Conventional least-squares 

methods for estimating parameters of the explanatory model, 

using either the entire sample or the subsample of uncensored 

observations, yield biased and inconsistent estimates when 

data are censored (Judge and others, 1985). However, a tobit 

regression model can be used with censored observations 

(Judge and others, 1985; Tobin, 1958). For this study, the 

parameters of the tobit regression model were estimated by 

using maximum likelihood methods, as described in Larson 

and others (2004). Maximum likelihood methods implemented 

Development site

Validation site

Figure 1. Location of sampling sites used for WARP model development (112 sites) and validation (25 sites).
4 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

Table 1. Watershed characteristics considered as explanatory variables for WARP models

[WARP, Watershed Regression for Pesticides; CIESIN, Center for International Earth Science Information Network; DEM, Digital Elevation Model; DMRK, 

DMRKYNETIC, Inc.; NLCD, National Land Cover data set; NRI, National Resources Inventory; NOAA, National Oceanic and Atmospheric Administration; 

STATSGO, State Soil Geographic data base; USLE, Universal Soil Loss Equation; cm, centimeter; d, day; km, kilometer; km2, square kilometer; m, meter; 

mm, millimeter; yr, year]

Variable Description

Pesticide Use

UI Atrazine agricultural-use intensity [(sum of use data derived from DMRK for applications to row crops, orchards and 

vineyards, and pasture and hay crops extrapolated to the basin scale)/basin area], kg/km2) 

Land use and population

AG Percent of basin with agricultural land use [(sum of NLCD categories row crops, small grains, pasture/hay, and orchards/

vineyards/other at 30-m cell resolution) × 100/basin area]

FOREST Percent of basin with forest land use [(sum of NLCD categories deciduous forest, evergreen forest, and mixed forest at 

30-m cell resolution) × 100/basin area]

URBAN Percent of basin with urban land use [(sum of NLCD categories low intensity residential, high intensity residential, 

commercial/industrial/transportation, and urban/recreational grasses at 30-m cell resolution) × 100/basin area]

POPDEN Mean 1990 population density in basin (people/km2) (CIESIN 1-km grid)

Agricultural management practices

ARTDRN Percent of the basin that is artificially drained [(aggregated representation of conservation practice categories 606, 607, and 

608 from the 1997 NRI polygons converted to 1-km cells) × 100/basin area] 

CONTILL Percent of the basin under conservation tillage [(1997 NRI conservation practice 329 polygons converted to 1-km cells) × 

100/basin area]

IRRI Percent of the basin that is irrigated [(aggregated representation of irrigation type categories 01, 02, and 03 from the 1997 

NRI polygons converted to 1-km cells) × 100/basin area]

TILE Percent of the basin that is drained by tiles [(1997 NRI conservation practice 606 polygons converted to 1-km cells) × 100/

basin area]

Soil Properties

AWC Mean available water capacity (fraction) in basin (STATSGO polygons converted to 1-km cells)

CLAY Mean percent clay in basin soils (STATSGO polygons converted to 1-km cells)

HGAB Mean percent of basin soils classified as hydrologic groups A and B (sum of hydrologic groups A and B from STATSGO 

polygons converted to 1-km cells)

HGCD Mean percent of basin soils classified as hydrologic groups C, D, and C/D (sum of hydrologic groups C, D, and C/D from 

STATSGO polygons converted to 1-km cells)

K Mean soil erodibility of uppermost soil horizon in basin (K-factor for USLE) (STASGO polygons converted to 1-km cells)

ORGM Mean percent organic matter in basin soils (STATSGO polygons converted to 1-km cells)

PERM Mean soil permeability, in cm/hr (STATSGO polygons converted to 1-km cells)

SAND Mean percent sand in basin soils (STATSGO polygons converted to 1-km cells)

SILT Mean percent silt in basin soils (STATSGO polygons converted to 1-km cells)

Physical watershed characteristics

WA Basin drainage area, in km2

ELEV Mean basin elevation, in m (DEM at 1-km cell resolution)

LATC Latitude of basin centroid, in decimal degrees

LONC Longitude of basin centroid, in decimal degrees

SLOPE Mean percent slope in basin (DEM at 1-km cell resolution)
Methods 5

Table 1. Watershed characteristics considered as explanatory variables for WARP models.—Continued

[WARP, Watershed Regression for Pesticides; CIESIN, Center for International Earth Science Information Network; DEM, Digital Elevation Model; DMRK, 

DMRKYNETIC, Inc.; NLCD, National Land Cover data set; NRI, National Resources Inventory; NOAA, National Oceanic and Atmospheric Administration; 

STATSGO, State Soil Geographic data base; USLE, Universal Soil Loss Equation; cm, centimeter; d, day; km, kilometer; km2, square kilometer; m, meter; 

mm, millimeter; yr, year]

Variable Description

Weather/climate characteristics

ADRY Mean annual number of consecutive dry days (NOAA data interpolated and converted to 1-km cells)

APPT Mean annual 1961–1990 precipitation, in cm/yr (NOAA data interpolated and converted to 1-km cells)

APPTI Mean annual precipitation intensity, in mm/d (NOAA data interpolated and converted to 1-km cells)

ATEMP Mean annual 1961–1990 temperature, in degrees Celsius (NOAA data interpolated and converted to 1-km cells)

WET Mean annual number of consecutive wet days (NOAA data interpolated and converted to 1-km cells)

R Mean annual 1971–2000 rainfall erosivity (R-factor for USLE) (NOAA station data analyzed and interpolated to 1-km cells)

PYEAR Total precipitation during the year of sampling (NOAA climate division precipitation data interpolated and converted to 

1-km cells)

PAPRJUN Total precipitation during April, May, and June of the year of sampling (NOAA climate division precipitation data 

interpolated and converted to 1-km cells)

PAPRSEP Total precipitation during April through September of the year of sampling (NOAA climate division precipitation data 

interpolated and converted to 1-km cells)

PMAY Total precipitation during May of the year of sampling (NOAA climate division precipitation data interpolated and converted 

to 1-km cells)

PMJN Total precipitation during May and June of the year of sampling (NOAA climate division precipitation data interpolated and 

converted to 1-km cells)

PMJL Total precipitation during May, June, and July of the year of sampling (NOAA climate division precipitation data interpolated and converted to 1-km cells)

PDYEAR Departure from the 1961–1990 average precipitation during the year of sampling (NOAA climate division precipitation data 

interpolated and converted to 1-km cells)

PDAJN Departure from the 1961–1990 average April-June precipitation during the year of sampling (NOAA climate division 

precipitation data interpolated and converted to 1-km cells)

PDMAY Departure from the 1961–1990 average May precipitation during the year of sampling (NOAA climate division precipitation 

data interpolated and converted to 1-km cells)

Hydrologic properties

CONTACT Mean subsurface contact time, in days (estimated by means of TOPMODEL hydrologic model [Wolock, 1993]

PERDUN Percent of basin streamflow contributed by Dunne overland flow (estimated by means of TOPMODEL hydrologic model)

PERHOR Percent of basin streamflow contributed by Horton overland flow (estimated by means of TOPMODEL hydrologic model)

PET Mean potential evapotranspiration, in cm (estimated using temperature data derived from Parameter-Elevation Regressions 

on Independent Slopes model (PRISM; Daly and others, 1997) and the Hamon PET equation (Hamon, 1961)

ROFF Mean annual 1951-1980 runoff , in cm/yr (USGS data interpolated and converted to 1-km cells)
6 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

in the survreg procedure (Therneau, 1999) in the statistical 

analysis program S-PLUS (Insightful Corporation, 1999) were 

used to estimate the parameters of the regression models. 

Measures of goodness of fit, such as the standard deviation of the residual error (commonly referred to as the root 

mean square error) or the coefficient of multiple determination 

(R2), used for conventional least-squares regression analysis, cannot be computed for the tobit regression model. The 

standard deviation of residual error is alternatively referred to 

as the “scale parameter” in maximum likelihood estimation. 

Estimates of the scale parameter from the maximum likelihood 

procedure provide only asymptotically unbiased estimates of 

the standard deviation of the residual error when estimated 

from sample data (Aitkin, 1981). These estimates, on average, underestimate the true standard deviation. The bias is 

a function of the sample size and degree of censoring, and 

it is expected to be minimal for models with lower percentages of censored observations. In the remainder of this report, 

biased estimates of the standard deviation of residual error 

are referred to as “scale” in figures and tables. Several pseudo 

R2 (pR2) measures suitable for use with the tobit regression 

model have been proposed in the literature as alternatives to 

R2. For this study, pR2 was calculated by using the method of 

Laitila (1993). As with conventional R2, the pR2 ranges from 

0 to 1 and is an estimate of the proportion of the variation in 

the response variable explained by the regression model (0 

indicates no variation explained; 1 indicates all variation is 

explained). 

Transformations of Response and Explanatory 

Variables

The maximum likelihood methods used for estimating the 

parameters of the regression models require that the relation 

between the variables be linear in the parameters and that the 

residual error be identically and normally distributed. Departures from these requirements can result in flawed estimates 

of model coefficients. One means of addressing departures 

from model assumptions is through transformation of either 

the response or the explanatory variables or both (Neter and 

others, 1985).

Various transformations were considered to minimize 

departures from the requirements of the maximum likelihood 

methods used. The logarithm of concentration was used as 

the response variable throughout model development. For the 

explanatory variables, logarithmic, square-root, and fourthroot transformations, as well as the untransformed value, were 

considered during development of the regression models.

Because the response variable is a logarithmic transformation, concentrations predicted by the model (after 

retransformation) are the median concentrations expected for 

sites that have a given set of explanatory values, rather than 

the mean concentrations. Predicted concentration statistics 

were not adjusted for transformation bias because estimates 

of median values of the statistics were considered more appropriate for the purposes of this study.

Selection of Explanatory Variables

A stepwise approach was used for selection of explanatory variables to include in the regression models. The StepAIC procedure of Venables and Ripley (1999), implemented 

in S-PLUS, was used for the initial selection of explanatory 

variables to include in the regression models. The StepAIC 

procedure, based on Akaike’s Information Criterion (Akaike, 

1974), balances model goodness of fit with the number of 

parameters needed to achieve that fit. The procedure attempts 

to quantify the concept of model parsimony by choosing 

simpler models over complex ones unless a complex model 

substantially improves the fit of the model.

Variable selection began with evaluation of transformations for each potential explanatory variable. The StepAIC 

procedure was used to evaluate each variable and its transformations together as a group. Untransformed and transformed 

explanatory variables were also evaluated by plotting them 

against the response variable to see whether the relation 

was approximately linear and whether residual errors were 

approximately identically and normally distributed. The StepAIC procedure and the plot evaluation determined whether 

the untransformed or transformed (logarithmic, square-root, 

or square) variable would be used in the explanatory variable 

selection process. 

Variable selection continued with 43 potential explanatory variables, too many to include in the StepAIC procedure 

at one time. So, groups of 12 variables were randomly selected 

and then evaluated with the StepAIC procedure. The random 

selection and evaluation of variables was repeated 1,000 times. 

Variables selected 100 percent of the time for all models by the 

StepAIC procedure were retained, and the random selection 

(groups of 12, including the retained variables) and StepAIC 

evaluation process was repeated again. Variables selected 

greater than 50 percent of the time by the StepAIC procedure 

for all repetitions were retained and typically represented less 

than half of the original number of explanatory variables. A 

subsampling, cross-validation approach was used to evaluate 

the remaining potential explanatory variables. The data were 

randomly divided into two subsamples of equal size, and the 

StepAIC procedure was used to evaluate the potential explanatory variables for both subsamples. The random subsampling 

and evaluation with the StepAIC procedure was repeated 

multiple times. Variables selected by the StepAIC procedure 

for a majority of the repetitions were retained and subjectively 

evaluated for reasonableness and their overall contribution to 

explaining the variation in atrazine-concentration statistics. 

In addition, potential explanatory variables not selected were 

individually evaluated to determine whether adding any of 

these variables would improve the ability of the model to 

explain variation in atrazine-concentration statistics without 

making the model overly complex. (See example in table 2.) 
Methods 7

Table 2. Example results from stepwise model development for selected WARP models.

[WARP, Watershed Regression for Pesticides; 21-day, annual maximum 21-day moving-average concentration model; pR2, pseudo R-square 

(R-squared value for tobit regression); Scale, tobit regression analogue of the root mean squared error obtained from ordinary least squares 

regression; model variables are defined in table 1]

Model

Annual mean 95th Percentile 21-day 

Scale pR2 Scale pR2 Scale pR2

(UI)1/4 0.54 0.68 0.62 0.62 0.60 0.67

(UI)1/4+PMJN 0.46 0.77 0.54 0.72 0.53 0.75

(UI)1/4+PMJN+K 0.44 0.79 0.51 0.74 0.51 0.76

(UI)1/4+PMJN+K+(WA)1/2 0.42 0.80 0.49 0.76 0.50 0.77

(UI)1/4+PMJN+K+(WA)1/2+log10R 0.41 0.81 0.48 0.77 0.50 0.78

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN 0.40 0.82 0.47 0.78 0.48 0.79

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN+WET 0.39 0.83 0.47 0.79 0.48 0.79

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN+ELEV 0.40 0.82 0.47 0.79 0.47 0.80

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN+HGAB 0.40 0.82 0.46 0.79 0.48 0.79

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN+PERM 0.40 0.83 0.46 0.79 0.48 0.79

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN+SLOPE 0.40 0.83 0.47 0.79 0.48 0.79

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN+HGCD 0.40 0.82 0.46 0.79 0.48 0.79

(UI)1/4+PMJN+K+(WA)1/2+log10R+PERDUN+ATEMP 0.40 0.82 0.47 0.79 0.48 0.79

The same explanatory variables were used in all the final 

models, and this approach simplifies the data requirements for 

application of the models. Extra consideration was given to 

variables selected for the 95th percentile, annual maximum, 

and the annual maximum moving-average models because the 

percentage of censored data was less for these concentration 

statistics than for lower-percentile concentration statistics. 

Coefficients of the models for the various concentration statistics were determined independently, and no constraint prevents 

the prediction of a concentration for a low percentile from 

exceeding the concentration prediction for a higher percentile. 

Use of the same explanatory variables in all the final models reduced how often this occurred. Statistical methods for 

estimating correlated models are not available for regressions 

with censored observations. 

Interactions between variables were not considered as 

potential explanatory variables because of the already large 

number of potential explanatory variables evaluated for the 

models. To have included interactions would have resulted in 

thousands of additional potential explanatory variables.
8 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

 Analysis of Model Fit

Diagnostics for censored regression (Escobar and 

Meeker, 1992) available in the survreg procedure in the 

S-PLUS program (Insightful Corporation, 1999) were used 

to assess influential observations. Variance inflation factors 

were used for detecting the presence of multicollinearity 

among explanatory variables (Neter and others, 1985). Box 

and whisker plots (Tukey, 1977) were used to qualitatively 

assess model performance (fig. 2). These boxplots summarize 

a group of data by showing a measure of central tendency (the 

median), the variation (interquartile range), the range (shown 

by the whiskers, which extend to 1.5 times the interquartile 

range when extreme values are present), and extreme values 

(lower/upper outside values shown by individual points). 

Boxplots were used for displaying the distribution of model 

residuals and comparing residuals among groups of data (for 

example, different regions of the country). 

Comparisons between predicted concentration statistics 

and concentration statistics computed from observations are 

made frequently in the discussion of model performance. 

Terms used for these comparisons are defined here for clarity. 

Concentration statistics generated by the WARP models are 

referred to as “predicted concentration statistics,” and concentration statistics computed from observations are referred to 

as “observed concentration statistics.” Predicted concentration 

statistics within a factor of 10, or “order of magnitude,” of the 

observed concentration statistics are between one-tenth and 10 

times the observed concentration statistic. For example, for an 

observed concentration statistic of 3 μg/L, predicted concentration statistics from 0.3 to 30 μg/L are within a factor of 10 

of the observed concentration statistic; predicted concentration 

statistics from 0.6 to 15 μg/L are within a factor of 5.



Outside value

(>1.5 x IQR)

Largest value

less than

1.5 times IQR

75th percentile

25th percentile

Median

IQR, interquartile range;

distance between 75th and

25th percentiles

Smallest value

less than

1.5 times IQR

Figure 2. Boxplot explanation.

Estimation of Prediction Intervals

Prediction intervals were approximated by using normal 

theory and the t-distribution; that is, methods for ordinary 

least square regression were used, but these methods yield 

only approximate results when applied to censored data. The 

standard errors were estimated from the maximum likelihood 

scale parameter by use of the adjustment suggested by Aitkin 

(1981).

Atrazine Models

The models for the selected percentiles, the annual mean, 

the annual maximum, and the annual maximum moving-average atrazine-concentration statistics have the form:

 (1)

where

UI is the atrazine-use intensity, annual 

agricultural use in the watershed (kilograms) 

divided by watershed area (square 

 kilometers);

PMJN is the total precipitation during May and 

June of the year of sampling (millimeters);

K is the soil erodibility factor (K-factor) from 

Universal Soil Loss Equation (USLE);

WA is the watershed area (square kilometers);

R is the rainfall erosivity factor (R-factor) 

from USLE; and

PERDUN is the percentage of streamflow due to 

Dunne overland flow.

Statistics for all models are given in table 3. Regression 

coefficients and statistics were based on model fit using tobit 

regression. Values of pR2 ranged from 0.63 for the 5th percentile model to 0.82 for the annual mean model, meaning that 

the models accounted for 63 to 82 percent of the variability in 

the concentration statistics among the 112 sites used for model 

development. Concentration statistics predicted by these models represent the expected median value of the concentration 

statistic for all sites with the same values for the explanatory 

variables. 

Analysis of Significant Explanatory Variables

Larson and others (2004) found that fourth-root of 

atrazine-use intensity in the watershed, as the fourth-root 

transformation, was the most important explanatory variable 

in the original WARP models. The same is true in the present update of the WARP models. Models using just (UI)1/4 as 
Table 3. Summary of statistics and coefficients for variables in the atrazine WARP models.[WARP, Watershed Regression for Pesticides; Pseudo R-square, R-squared value for tobit regression; Scale, tobit regression analogue of the root mean squared error obtained froregression; Max, maximum; MA, moving average; model variables are defined in table 1; <, less than]

Atrazine-concentration statistic

Regression coefficients (p-value) Pseudo

R-square 

Intercept UI1/4 PMJN K (WA)1/2 log10R PERDUN

5th -3.79

(<0.05)

0.71

(<0.05)

0.0010

(0.19)

1.21

(0. 14)

0.00035

(0.10)

0.239

(0.27)

-0.124

(<0.05)

0.63 

10th -3.91

(<0.05)

0.67

(<0.05)

0.0010

(0.15)

1.44

(0.06)

0.00038

(0.06)

0.336

(0.10)

-0.135

(<0.05)

0.66 

15th -3.91

(<0.05)

0.63

(<0.05)

0.0010

(0.12)

1.38

(<0.05)

0.00042

(<0.05)

0.429

(<0.05)

-0.141

(<0.05)

0.69 

25th -3.79

(<0.05)

0.61

(<0.05)

0.0008

(0.17)

1.09

(0.09)

0.00043

(<0.05)

0.469

(<0.05)

-0.122

(<0.05)

0.70 

50th -3.73

(<0.05)

0.58

(<0.05)

0.0012

(<0.05)

1.25

(<0.05)

0.00047

(<0.05)

0.473

(<0.05)

-0.100

(<0.05)

0.74 

75th -3.77

(<0.05)

0.63

(<0.05)

0.0020

(<0.05)

1.62

(<0.05)

0.00063

(<0.05)

0.489

(<0.05)

-0.083

(<0.05)

0.79 

85th -3.91

(<0.05)

0.71

(<0.05)

0.0024

(<0.05)

2.22

(<0.05)

0.00066

(<0.05)

0.506

(<0.05)

-0.080

(<0.05)

0.79 

90th -3.91

(<0.05)

0.69

(<0.05)

0.0027

(<0.05)

2.39

(<0.05)

0.00070

(<0.05)

0.554

(<0.05)

-0.086

(<0.05)

0.80 

95th -4.01

(<0.05)

0.71

(<0.05)

0.0024

(<0.05)

2.89

(<0.05)

0.00058

(<0.05)

0.665

(<0.05)

-0.091

(<0.05)

0.78 

Annual mean -4.23

(<0.05)

0.75

(<0.05)

0.0022

(<0.05)

2.40

(<0.05)

0.00046

(<0.05)

0.571

(<0.05)

-0.096

(<0.05)

0.82 

Annual maximum -3.36

(<0.05)

0.91

(<0.05)

0.0025

(<0.05)

1.77

(<0.05)

0.00016

(0.48)

0.591

(<0.05)

-0.122

(<0.05)

0.76 

Annual maximum 21-day MA -3.60

(<0.05)

0.86

(<0.05)

0.0026

(<0.05)

2.07

(<0.05)

0.00039

(0.06)

0.582

(<0.05)

-0.112

(<0.05)

0.79 

Annual maximum 60-day MA -3.86

(<0.05)

0.81

(<0.05)

0.0026

(<0.05)

2.26

(<0.05)

0.00049

(<0.05)

0.610

(<0.05)

-0.110

(<0.05)

0.81 

Annual maximum 90-day MA -3.91

(<0.05)

0.80

(<0.05)

0.0025

(<0.05)

2.23

(<0.05)

0.00051

(<0.05)

0.597

(<0.05)

-0.105

(<0.05)

0.81 
10 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

an explanatory variable account for 56 to 68 percent of the 

variability in concentration statistics among the 112 sites. The 

coefficient of the transformed use-intensity variable is positive 

in all models.

Total precipitation during May and June of the year of 

sampling (PMJN) is significant (p<0.05) in the higher percentile, annual mean, annual maximum, and annual maximum 

moving-average models (table 3). For the majority of sites, 

the months of April through June include the highest atrazine 

application period. High precipitation during May and June, 

after atrazine application, increases atrazine runoff to streams. 

This increase is consistent with the positive coefficient 

observed for PMJN in all models.

Soil erodibility factor (K-factor) is significant (p<0.05) in 

all the models except for the 5th-, 10th-, and 25th-percentile 

models (table 3). The coefficient for K-factor is positive in all 

models. The K-factor is an index that indicates the susceptibility of soil to erosion (Brooks and others, 1991). The positive 

coefficient for K-factor in the regression models implies that 

more runoff of atrazine occurs in areas susceptible to soil erosion. As discussed in Larson and others (2004), this does not 

mean that atrazine is being transported on eroded soil particles, because atrazine is primarily transported in the dissolved 

phase in runoff water. The positive relation between K-factor 

and atrazine concentrations in the stream may indicate that 

areas with high K-factor soils have a relatively high potential 

for the formation of runoff during precipitation, resulting in an 

increased potential for transport of atrazine from application 

areas to streams.

Watershed area, as the square-root transformation, was 

significant (p<0.05) in all models except for the 5th- and 

10th- percentile, annual maximum, and annual maximum 

21-day moving-average models (table 3). The coefficient for 

watershed area was positive in all models. Larson and Gilliom 

(2001) listed three factors that may contribute to the positive 

relation between watershed area and atrazine concentration in 

streams. First, the contribution of water from multiple small 

streams to a large river can result in elevated concentrations 

that are sustained for a longer time than in the individual 

streams because the timing of pesticide application and subsequent runoff can vary among the individual streams. Second, 

in some small streams, the highest concentrations may not 

have been sampled because concentrations of pesticides may 

remain elevated for relatively short periods. The computed 

concentration statistics for the higher percentiles and annual 

maximum for these small streams may be biased low, further 

strengthening the positive relation between watershed area 

and concentration reflected in the coefficients for the watershed-area variable in the regression models. Third, for some 

streams with large watersheds, pesticide use is concentrated 

in a relatively small farmed area of the watershed, wheras the 

remainder of the watershed is largely untreated. The use-intensity value for these watersheds may be relatively low, despite 

high use in the farmed areas. It is likely that more water is 

contributed to the stream through runoff from the farmed areas 

than from nonagricultural areas, especially in arid regions 

where crops are irrigated. For these cases, the watershed-area 

variable may act as a correction factor, accounting for the 

high concentrations in these streams despite relatively low 

use-intensity values for the watersheds as a whole (Larson and 

others, 2004). 

Rainfall erosivity (R-factor), as the logarithm transformation, was significant (p<0.05) in all models except for the 

5th- and 10th- percentile models (table 3). The coefficient for 

R-factor was positive in all models. The R-factor is an index 

that characterizes the energy of storms in a specific area, averaged over a number of years (Brooks and others, 1991). The 

positive coefficient implies that more runoff of atrazine occurs 

in areas of high-energy rainstorms. 

Percent Dunne overland flow (PERDUN) was significant 

(p<0.05) in all the models (table 3). The coefficient for PERDUN was negative in all models. PERDUN is the percentage 

of streamflow generated by overland flow from precipitation 

falling on saturated soils (Dunne overland flow). The negative 

coefficient for PERDUN implies that atrazine concentrations 

were lower in areas where Dunne overland flow contributions to streamflow were high. Generally, near-stream lowland 

and wetland areas that quickly become saturated contribute 

more Dunne overland flow to streams than areas that are 

well drained. Lowland and wetland areas are less likely to 

be cropped and less likely to have been treated with atrazine. 

Dunne overland flow may be a source of relatively atrazinefree water and may act to dilute concentrations of atrazine in 

streams that receive a higher percentage of Dunne overland 

flow (Larson and Gilliom, 2001). 

Model Performance

The discussion, tables, and figures describing model 

performance in this report are for the 95th-percentile, annual 

mean, annual maximum, and annual maximum 21-, 60-, and 

90-day moving-average WARP models. Model performance 

was evaluated by assessment of (1) goodness of fit and residual errors for the model-development sites, (2) residual errors 

for the model-validation sites, and (3) uncertainty in model 

predictions. Each is described below.

Model-Development Sites

Regression-model results for the model-development 

sites are shown in figure 3. For the models shown, most of the 

predicted concentration statistics are within a factor of 10 of 

the observed concentration statistics at the development sites. 

Residual errors were calculated by subtracting the 

predicted concentration statistics from the observed concentration statistics; therefore, residual errors less than zero indicate 

overprediction of the concentration statistic, and residual 

errors greater than zero indicate underprediction. A residual 

error of zero indicates exact agreement, and residuals of -1 

and +1 indicate that the prediction for a given site was 10 

times and one-tenth of the observed concentration statistic, 
Atrazine Models 11

–4

–2

0

2

pR–square = 0.82

Scale = 0.40

Annual mean

–4

–2

0

2

pR–square = 0.78

Scale = 0.47

95th percentile

–4

–2

0

2

pR–square = 0.76

Scale = 0.54

Annual maximum

LO

G

A

RIT

H

M

OF O

B

SER

VED

A

T

R

A

ZINE-C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC, IN

MIC

R

O

G

R

A

M

S PER LITER

–4

–2

0

2

pR–square = 0.79

Scale = 0.48

Annual maximum

21–day moving average

–4

–2

0

2

pR–square = 0.81

Scale = 0.45

Annual maximum

60–day moving average

–4 –3 –2 –1 0 1 2

LOGARITHM OF PREDICTED ATRAZINE-CONCENTRATION STATISTIC,

IN MICROGRAMS PER LITER

–4

–2

0

2

pR–square = 0.81

Scale = 0.43

Annual maximum

90–day moving average

–4 –3 –2 –1 0 1 2

Figure 3. Atrazine-concentration statistics from concentrations observed at model-development sites in relation to values of the 

same statistics predicted by the WARP models. The black line is a 1:1 line, indicating exact agreement of the observed and predicted 

concentration statistics. The red lines are plus and minus 1 log unit from the 1:1 line; the area between a red and black line represents 

predicted concentration statistics within a factor of 10 of the observed concentration statistics. The green lines represent the 

95-percent prediction interval. Filled symbols indicate censored observed concentration statistics. pR-square is the pseudo R-squared 

(R-square value used for tobit regression).
12 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

respectively. Boxplots of the residual errors from the models 

show the performance of the models with respect to watershed 

area (fig. 4) and geographic region (fig. 5) of the model-development sites. 

Predictions of atrazine-concentration statistics are not 

substantially biased with respect to watershed area, as indicated by figure 4. For this plot, watershed areas of the 112 

model-development sites were divided into quintiles (that is, 

the first quintile contains the smallest 20 percent of watershed 

areas, the second quintile contains the second smallest 20 percent, and so on). Residual errors are similar for all five groups, 

which together cover more than five orders of magnitude in 

watershed area (17 to 2,965,000 km2).

Predictions of atrazine-concentration statistics have no 

substantial regional bias, as indicated by figure 5. The regional 

groupings used for this plot were based on the U.S. Department of Agriculture Farm Resource Regions (U.S. Department 

of Agriculture, 2000), which categorized agricultural regions 

of the conterminous United States on the basis of climate, 

topography, soil types, and dominant agricultural activities. 

The nine Farm Resource Regions were consolidated into five 

regions for the present study (fig. 6) so that each region would 

have sufficient sites for the computation of statistics summarizing model fit. 

Model-Validation Sites

The selected models were applied to concentration data 

from the 25 validation sites, and plots of the predictions for 

these sites are shown in figure 7. Predicted concentration statistics from the selected models are within a factor of 10 of the 

observed concentration statistics for most of the sites. 

The updated WARP model predictions for the validation 

sites are not substantially biased (fig. 8), even though predictions from the previous WARP models were biased low for 

most of the validation sites (Larson and others, 2004; Stone 

and others, 2008). The incorporation of the annual atrazine-use 

estimates and the PMJN explanatory variable into the WARP 

models appears to have reduced the low prediction bias seen 

with application of the previous WARP models to the validation-site data. 

Uncertainty in Model Predictions

Uncertainty in the prediction of a concentration statistic 

can be expressed in terms of a prediction interval (PI) for a 

specified confidence level—the confidence level used in this 

study is 95 percent. Conceptually, each predicted concentration statistic is the median estimate of the particular concentration statistic (percentiles, annual mean, annual maximum, and 

annual maximum moving averages) for all the stream sites that 

have the same combination of values for the explanatory variables. The PI is the range of values for a concentration statistic 

within which 95 percent of the actual concentration-statistic 

values are expected to occur for all stream sites with the same 

values of explanatory variables. In addition, the PI can be 

interpreted as the range within which the actual concentration 

statistic for an individual site and year is expected to fall 95 

percent of the time. 

The PIs for the selected models are shown as green lines 

in figure 3. Concentration statistics are expressed as logarithms, resulting in symmetrical intervals for the PI. (The high 

and low bounds of the intervals are the same distance from 

the predicted value.) However, expressing the concentration 

statistics as logarithms obscures the fact that the intervals are 

skewed—the upper part of the PI interval covers a wider range 

of values than the lower part. Comparison of the PIs (green 

lines) to the predicted concentration statistics within a factor of 10 of the observed concentration statistics (red lines) 

shows that the PIs are largest for the annual maximum model 

and smallest (PIs less than a factor of 10) for the annual mean 

model (fig. 3). The PI is a function of the fit of the model and 

the amount of variability explained by the model. The annual 

mean model has a better model fit and explains more variability (scale = 0.40; pR2

 = 0.82) than the annual maximum model 

does (scale = 0.54; pR2 = 0.76). 

The levels of uncertainty for the 95th-percentile, annual 

mean, annual maximum, and annual maximum movingaverage (21-, 60-, and 90-day) models are compared in figure 

9. PI size is represented as the ratio of the upper boundary of 

the interval to the predicted concentration statistic. (This is the 

same as the ratio of the predicted concentration statistic to the 

lower boundary of the interval.) The extreme values (shown as 

asterisks) represent the same sites in each model, which have 

one or more explanatory variables that are relatively extreme 

in value when compared to values for the rest of the sites. For 

example, the highest ratio for each of the models, shown by 

highest asterisks in each boxplot, is for the Mississippi River 

at St. Francisville, La., which has a watershed area far greater 

than the rest of the sites. PIs for the annual maximum model 

extend to a factor of about 13 to 14 above and below the predicted annual maximum concentration statistic for most sites, 

whereas PIs for all other models extend to a factor of only 

7 to 11.
Atrazine Models 13

–3

–1

1

3

Annual mean

–3

–1

1

3

95th percentile

–3

–1

1

3

Annual maximum

–3

–1

1

3

Annual maximum 21–day moving average

–3

–1

1

3

Annual maximum 60–day moving average

–3

–1

1

3

Annual maximum 90–day moving average

All data

N

=

112

Quintile

1

N

=

22

Quintile

2

N

=

23

Quintile

3

N

=

22

Quintile

4

N

=

23

Quintile

5

N

=

22

LO

G

A

RIT

H

M

OF T

HE O

B

SER

VED

C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC

MIN

U

S

LO

G

A

RIT

H

M

OF PREDIC

TED

C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC

Figure 4. Residual errors for atrazine-concentration statistics predicted by WARP models in relation to watershed area for the 112 

model-development sites. First boxplot in each row shows the residual errors for all 112 sites. Remaining boxplots show the residual 

errors for sites grouped into five classes based on quintiles of watershed area. Residual error is [log10(observed concentration statistic) 

– log10(predicted concentration statistic)]. See figure 2 for boxplot explanation. 
14 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

Figure 5. Residual errors for atrazine-concentration statistics predicted by the WARP models in relation to geographic region for 

the 112 model-development sites. First boxplot in each row shows the residual errors for all 112 sites. Remaining boxplots show the 

residual errors for sites grouped by region. Regions are shown in figure 6. Residual error is [log10(observed concentration statistic) – 

log10(predicted concentration statistic)]. See figure 2 for boxplot explanation. 

–3

–1

1

3

Annual mean

–3

–1

1

3

95th percentile

–3

–1

1

3

Annual maximum

–3

–1

1

3

Annual maximum 21–day moving average

–3

–1

1

3

Annual maximum 60–day moving average

–3

–1

1

3

Annual maximum 90–day moving average

All data

N

=

112

Region

1

N

=

30

Region

2

N

=

31

Region

3

N

=

17

Region

4

N

=

11

Region

5

N

=

23

LO

G

A

RIT

H

M

OF T

HE OBSER

VED

C

O

N

CEN

T

R

A

TIO

N ST

A

TISTIC

MIN

US

LO

G

A

RIT

H

M

OF PREDIC

TED

C

O

N

CEN

T

R

A

TIO

N ST

A

TISTIC
Atrazine Models 15

Northern Crescent

Eastern Uplands

Heartland

Southern Seaboard

Mississippi Portal

NorthernGreat Plains

Prairie Gateway

Basin and Range

Fruitful Rim

Region for WARP analysis

EXPLANATION

4

Farm Resource Region

Figure 6. U.S. Department of Agriculture Farm Resource Regions (U.S. Department of Agriculture, 2000) and regions used for 

evaluation of the WARP models.
16 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

–4

–2

0

2

Annual mean

–4

–2

0

2

95th percentile

–4

–2

0

2

Annual maximum

LO

G

ARITHM

OF O

BSERVED

A

TRA

ZINE-CO

NCEN

TRA

TIO

N ST

A

TISTIC, IN

MICRO

GRA

MS PER LITER

–4

–2

0

2

Annual maximum

21–day moving average

–4

–2

0

2

Annual maximum

60–day moving average

–4 –3 –2 –1 0 1 2

LOGARITHM OF PREDICTED ATRAZINE-CONCENTRATION STATISTIC,

IN MICROGRAMS PER LITER

EXPLANATION

Model development sites

Censored observations for model development sites

Model validation sites

Censored observations for model validation sites

–4

–2

0

2

Annual maximum

90–day moving average

–4 –3 –2 –1 0 1 2

Figure 7. Atrazine-concentration statistics observed at model-validation sites in relation to values of the same statistics predicted 

by the WARP models. Values for the model-development sites are shown for comparison. The black line is a 1:1 line, indicating exact 

agreement of the observed concentration statistics and predicted concentration statistics; the area between a red and black line 

represents predicted concentration statistics within a factor of 10 of the observed concentration statistics. 
Atrazine Models 17

–2

0

2

New models and validation sites

–2

0

2

Previous models and validation sites

Annual

mean

95th

percentile

Annual

maximum

Annual

maximum

21–day

moving

average

Annual

maximum

60–day

moving

average

Annual

maximum

90–day

moving

average

LO

G

A

RIT

H

M

OF T

HE O

B

SER

VED

C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC

MIN

U

S

LO

G

A

RIT

H

M

OF PREDIC

TED

C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC

Figure 8. Residual errors for atrazine-concentration statistics predicted by the updated WARP models 

and the previous WARP models using data from 25 model-validation sites. Residual error is [log10(observed 

concentration statistic) – log10(predicted concentration statistic)]. See figure 2 for boxplot explanation.
18 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

Comparison to Previous WARP Models

The updated WARP models include the same five 

explanatory variables and transformations that were used in 

the original WARP models: use intensity, K-factor, watershed 

area, R-factor, and percentage of streamflow from Dunne 

overland flow. However, the atrazine-use-intensity variable in 

the updated WARP model is now based on annual atrazine-use 

estimates rather than estimates for 1992 or 1997. The updated 

WARP models also include a sixth explanatory variable, total 

precipitation during May and June of the year of sampling. 

Addition of this sixth explanatory variable explains more variability without making the models overly complex. 

Figure 9. Potential prediction errors for concentration statistics predicted with the atrazine WARP models. 

Potential prediction error is represented by the ratio of the upper boundary of the 95-percent prediction 

interval to the predicted atrazine-concentration statistic. Each boxplot shows the ratios among the 112 modeldevelopment sites. See figure 2 for boxplot explanation. 

Annual

mean

95th

percentile

Annual

maximum

Annual

maximum

21–day

moving

average

Annual

maximum

60–day

moving

average

Annual

maximum

90–day

moving

average

5

10

15

20

25

R

A

TIO

OF UPPER

B

O

U

N

D

A

R

Y

OF PREDIC

TIO

N IN

TER

V

AL T

O

PREDIC

TED

A

T

R

A

ZINE C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC

Model performance values (pR2 and scale parameter) 

for the 5th-percentile updated WARP model and the previous 

WARP model are comparable. However, the updated models 

for the other atrazine-concentration statistics show improvement over the previous WARP models in terms of pR2

 and the 

scale parameter (table 4). 

Differences in model performance between the updated 

and previous WARP models can be seen in terms of their 

residual errors, as shown in figure 10. The updated WARP 

models generally have a smaller variation (interquartile range 

shown by the box) and range (shown by the whiskers) of 

residual errors than the previous WARP models. 
Atrazine Models 19

Model Limitations

Use of the regression models for predicting atrazine-concentration statistics, and the WARP methodology in general, 

are subject to several limitations that are important to consider 

when applying the models.

1. The sampling frequencies of the model-development sites 

were not sufficient to reliably characterize the highest 

concentrations during a year. Thus, application of the 

models to predict the annual maximum concentration is 

expected to underpredict the actual annual maximums.

2. The regression models are designed for prediction of atrazine-concentration statistics for streams of the conterminous United States. Although the 112 sites used for model 

development represent a wide variety of environmental 

settings and a large range of watershed areas, it is likely 

that other watersheds within the United States have one or 

more characteristics outside the ranges of the watershed 

parameters used to develop the models. Application of the 

models to streams draining such watersheds would result 

in increased uncertainty in predicted concentrations.

3. The models were developed by using concentration 

data from streams. Application of the models to lakes or 

reservoirs will result in biased predictions because of the 

influence of water storage on the temporal distribution of 

concentrations.

4. The atrazine-use data used in the models are estimates for 

applications to agricultural land only. Substantial nonagricultural use of atrazine in a watershed could result in 

underprediction of atrazine concentrations in a stream, if 

such use cannot be estimated. 

Table 4. Comparison of updated WARP models to previous WARP models.

[WARP, Watershed Regression for Pesticides; MA, moving average ; Pseudo R-square, R-squared value for tobit 

regression; Scale, tobit regression analogue of the root mean squared error obtained from ordinary least squares 

regression. Previous WARP pseudo R-square and scale values were obtained from Larson and others (2004) and 

Stone and others (2008).]

Atrazine-concentration statistic

Updated WARP Previous WARP

Pseudo

R-square Scale PseudoR-square Scale

5th 0.63 0.50 0.62 0.50

10th 0.66 0.47 0.63 0.49

15th 0.69 0.43 0.65 0.45

25th 0.70 0.40 0.64 0.44

50th 0.74 0.37 0.68 0.41

75th 0.79 0.37 0.70 0.44

85th 0.79 0.43 0.71 0.50

90th 0.80 0.43 0.73 0.50

95th 0.78 0.47 0.73 0.52

Annual mean 0.82 0.40 0.77 0.46

Annual maximum 0.76 0.54 0.72 0.58

Annual maximum 21-day MA 0.79 0.48 0.74 0.54

Annual maximum 60-day MA 0.81 0.45 0.75 0.51

Annual maximum 90-day MA 0.81 0.43 0.75 0.50
20 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

Figure 10. Residual errors for atrazine-concentration statistics predicted by the updated WARP models 

and the previous WARP models. Residual error is [log10(observed concentration statistic) – log10(predicted 

concentration statistic)]. See figure 2 for boxplot explanation. 

–2

0

2

New models and development sites

–2

0

2

Previous models and development sites

Annual

mean

95th

percentile

Annual

maximum

Annual

maximum

21–day

moving

average

Annual

maximum

60–day

moving

average

Annual

maximum

90–day

moving

average

LO

G

A

RIT

H

M

OF T

HE O

B

SER

VED

C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC

MIN

U

S

LO

G

A

RIT

H

M

OF PREDIC

TED

C

O

N

CEN

T

R

A

TIO

N

S

T

A

TIS

TIC
References Cited 21

Summary and Conclusions

Regression models were developed for predicting 

selected percentiles (5th, 10th, 15th, 25th,50th, 75th, 85th, 

90th, and 95th), annual mean, annual maximum, and annual 

maximum moving-average (21-, 60-, and 90-day) atrazineconcentration statistics for streams, using nationally available 

data on watershed characteristics and atrazine use. Concentration statistic predictions generated by these models can be 

used to characterize the levels of atrazine for comparison to 

specific water-quality benchmarks for evaluation of potential 

concerns regarding human health and aquatic life.

Separate models were developed for the selected percentiles, annual mean, annual maximum, and the annual maximum moving-average concentration statistics. The updated 

WARP models include the same five explanatory variables and 

transformations that were used in the original WARP models. 

The updated WARP models also include a sixth explanatory 

variable, total precipitation during May and June of the year 

of sampling. Comparison of pR2, scale parameter, and residual 

errors between the updated WARP models and previous 

WARP models show better model performance for the updated 

models.

The models accounted for 63 to 82 percent of the variability in concentration statistics among the 112 sites used 

for model development. Uncertainty in predicted concentration statistics was expressed in terms of prediction intervals. 

For the 95th percentile, annual mean, annual maximum, and 

annual maximum moving-average (21-, 60-, and 90-day) 

models, 95-percent prediction intervals extend to a factor of 7 

to 14 above and below the predicted concentration statistic in 

most cases. 

Results for the 25 model validation sites with the 95th 

percentile, annual mean, annual maximum, and annual maximum moving-average (21-, 60-, and 90-day) models show that 

concentration statistics were predicted within a factor of 10 of 

the observed concentration statistics in nearly all cases. The 

updated WARP models show no substantial bias for the validation sites, an improvement over the previous WARP models 

that showed a low bias for the same validation sites. 

Adequate monitoring of pesticide concentrations in 

the numerous streams in the United States is prohibitively 

expensive, particularly for sampling frequencies that are high 

enough to reliably estimate the concentrations statistics needed 

for risk assessments. The WARP models are tools for predicting pesticide concentrations in unmonitored or inadequately 

monitored streams. Overall, performance of the models for 

the development and validation sites supports the application 

of the WARP models for predicting atrazine concentrations in 

streams and provides a framework to interpret the predictions 

in terms of uncertainty. For streams where direct measurements of atrazine are lacking, the WARP model predictions for 

atrazine concentrations can be used to characterize the probable levels of atrazine for comparison to specific water-quality 

benchmarks. 

References Cited

Acetochlor Registration Partnership, 2009, Acetochlor Registration Partnership (ARP) Web site, accessed April 13, 2009 

at http://www.arpinfo.com

Aitken, Murray, 1981, A note on the regression analysis of 

censored data: Technometrics, v. 23, no. 2, p. 161–163.

Akaike, Hirotugu, 1974, A new look a the statistical model 

identification: IEEE Transactions on Automatic Control, 

v. 19, no. 6, p. 716–723.

Brooks, K.N., Ffolliott, P.F., Gregersen, H.M., and Thames, 

J.L., 1991, Hydrology and the management of watersheds: 

Ames, Iowa, Iowa State University Press, 392 p.

Daly, C., Taylor, G., and Gibson, W., 1997, The PRISM 

approach to mapping precipitation and temperature: Proceedings of the Conference on Applied Climatology, Reno, 

Nev., American Meteorology Society, p. 208–209.

Escobar, L.A., and Meeker, W.Q., Jr., 1992, Assessing influence in regression analysis with censored data: Biometrics, 

v. 48, no. 2, p. 507–528.

Gilliom and others, 2006, The quality of our Nation’s waters—

Pesticides in the Nation’s streams and ground water, 

1992–2001: U.S. Geological Survey Circular 1291, 172 p.

Hamon, W.R., 1961, Estimating potential evapotranspiration: 

Journal of the Hydraulics Division, Proceedings of the 

American Society of Civil Engineers, v. 87, no. HY3, 

p. 107–120.

Insightful Corporation, 1999, S-PLUS 2000 user’s guide: 

Seattle, Wash., Insightful Corporation, Data Analysis 

Products Division, 558 p. 

Judge, C.G., Griffiths, W.E., Carter-Hill, R., Lutkepohl, H., 

and Lee, T.C., 1985, The theory and practice of econometrics (2d ed.): New York, John Wiley and Sons, 1,019 p.

Kiely, Timothy; Donaldson, David; and Grube, Arthur, 2004, 

Pesticide industry sales and usage 2000 and 2001 market 

estimates: U.S. Environmental Protection Agency Office of 

Pesticide Programs, EPA–733–R–04–001, 33 p.

Laitila, Thomas, 1993, A pseudo-R2 measure for limited and 

qualitative dependent variable models: Journal of Econometrics, v. 56, no. 3, p. 341–355.

Larson, S.J., Crawford, C.G., and Gilliom, R.J., 2004, Development and application of Watershed Regressions for 

Pesticides (WARP) for estimating atrazine concentration 

distributions in streams: U.S. Geological Survey WaterResources Investigations Report 03–4047, 68 p.
22 Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams

Larson, S.J., and Gilliom, R.J., 2001, Regression models for 

estimating herbicide concentrations in U.S. streams from 

watershed characteristics: Journal of the American Water 

Resources Association, v. 37, no. 5, p. 1349-1367.

Neter, John; Wasserman, William; and Kutner, M.H., 1985, 

Applied linear statistical models (2d ed.): Homewood, Ill., 

Irwin, 1,127 p.

Richards, R.P., and Baker, D.B., 1993, Pesticide concentration 

patterns in agriculture drainage networks in the Lake Erie 

Basin: Environmental Toxicology and Chemistry, v. 12, 

no. 1, p. 13–26.

Richards, R.P., Baker, D.B., Kramer, J.W., and Ewing, D.E., 

1996, Annual loads of herbicides in Lake Erie tributaries in 

Ohio and Michigan: Journal of Great Lakes Research, v. 22, 

no. 2, p. 414–428.

Shelton, L.R., 1994, Field guide for collecting and processing stream-water samples for the National Water-Quality 

Assessment Program: U.S. Geological Survey Open-File 

Report 94–455, 42 p.

Stone, W.W., Gilliom, R.J., and Crawford, C.G., 2008, Watershed Regressions for Pesticides (WARP) for predicting 

annual maximum and maximum moving-average concentrations of atrazine in streams: U.S. Geological Survey OpenFile Report 2008–1186, 19 p.

Therneau, T.M., 1999, A package of survival routines for S: 

Rochester, Minn., Mayo Clinic Department of Health 

Sciences Research, 83 p.

Tobin, J., 1958, Estimation of relationships for limited 

dependent variables: Econometrics, v. 26, p. 24–36.

Tukey, J.W., 1977, Exploratory data analysis: Reading, Mass., 

Addison-Wesley, 688 p.

U.S. Department of Agriculture, 2000, Farm resource regions: 

U.S. Department of Agriculture, Economic Research 

Service, Agricultural Information Bulletin Number 760 

[also available at http://www.ers.usda.gov/publications/

aib760/aib-760.pdf]

U.S. Geological Survey, 1997 to present, National field manual for the collection of water-quality data: U.S. Geological 

Survey Techniques of Water-Resources Investigations, 

book 9, chaps. A1–A9, 2 v., [variously paged].

Venables, W.N., and Ripley, B.D., 1999, Modern applied 

statistics with SPLUS (3d ed.): New York, Springer-Verlag, 

501 p.

Wolock, D.W., 1993, Simulating the variable-source-area 

concept of streamflow generation with the watershed model 

TOPMODEL: U.S. Geological Survey Water-Resources 

Investigations Report 93–4124, 33 p.

Stone and Gilliom—Update of Watershed Regressions for Pesticides (WARP) for Predicting Atrazine Concentration in Streams—Open-File Report 2009–1122