Abstract:

A process for integrating field studies, information from remotely sensed data, watershed models, and GIS, was used to accurately map and monitor the dispersion of pyrite rich mine waste and tailings, resulting from historic small to medium scale mining in the Patagonia and southern Santa Rita Mountains District, located in southern Arizona. Areas selected for this study include the drainage basins for Upper Harshaw Creek, Lower Harshaw Creek, Cox Gulch_Three-R Canyon, Alum Gulch_ Flux Canyon, and Providencia Canyon. Numerous underground mines, surface prospects, and associated waste dumps and tailing piles are abandoned within these basins however, large parts of the region are mineralized and remain largely undisturbed. Altered but unmined exposures contain assemblages characterized by pyrophyllite-sericite, clay alteration zones dominated by alunite-kaolinite and limited quart-sericite-pyrite zones. Maps of soils, sediments, water, topography and vegetation were used as input into models to account for the dispersion of waste piles and erosion potential of mineralized rock and soils in mountainous watersheds.

Water is an important and valuable resource in this semiarid environment. The issue is exemplified by problems associated with the Comprehensive Environmental Response Compensation and Liability Act of 1980 (CERCLA) efforts to provide cleanup to degraded areas. Sites of degraded water (low pH, elevated metals) are coincident and influenced by manmade and natural processes. Distinct enhancements in contaminant levels of water and sediments are indicated down-gradient of sulfide-bearing waste piles situated in or near stream channels. Erosion of historic piles has been a primary factor in the mechanical breakup and distribution of detrital material deposited within the streambed. The anthropogenic imprint is manifested by the documented increase in dissolved metal cations in water and increase in acid generation potential along the length of the stream channel. This study will attempt to identify erosion rates and net sediment delivery of soil and waste/tailings and compare with measured pH ranges within several watershed regions. A watershed scale erosion prediction model (USLE) and sediment delivery model (SEDMOD) are described and used within a GIS to estimate sediment yield and deposition within 5 separate watershed systems. Spatial analysis tools allow for the calculations to be repeated across the watersheds at a 30-meter resolution. The model results will provide insight into the non-point sources and sinks of sediment yield to the drainage system. Several alternative ways in which the pattern could be used to develop mitigation strategies are discussed.

Introduction

Patagonia and the southern Santa Rita Mountains area, located in southern Arizona (see Figure 1) was mined intermittently from the 1600’s to the mid- 1960’s primarily for silver, lead, copper and zinc as well as gold. The studies being done require an understanding of the movement of water through these mined areas and how alterations may affect that flow. Several of the mines have been identified in a list of high priority environmental degradation sites by CERCLA within the Coronado National Forest. Research has found an average pH reading for the disturbed areas of 3.5. This predominates in most surface water sampled as well as the groundwater sampled at residents’ wells. A low pH means that the water is acidic. Mine waters are often "acid" because of the common association of the sulfide pyrite with most metal ores and many solid fuels. Pyrite, as well as a number of other ore materials, rapidly decomposes when broken and exposed to moisture and air, eventually producing sulfuric acid. This chemical reaction occurs spontaneously and the acid mine water then has the ability to leach other pollutants in the solution (USDA, U. S. Forest Service 1977).

This study focuses on locating areas contributing pollutants to streams, areas containing large deposits of this type of pollution, and an estimate of concentrations of these sediment-loading and deposition sites. Basic processes of hydrological modeling in the area were determined using the spatial analysis tools available in Geographic Information Systems (GIS). Once these basic characteristics were known, watershed models could be applied to create accurate estimates of erosion potential and sediment delivery. This data is then queried and analyzed to derive the specific Non Point Source (NPS) sources and sinks within the watershed.

Watershed Parameters

A watershed is defined as the total drainage basin or catchment area flowing into a given outlet (pour point). A digital elevation data set is the necessary to delineate watersheds. A set of six 30-meter resolutions Digital Elevation Model’s (DEM’s) were mosaiced together in ARC/INFO to form the surface model DEM that would encompass the study area. The DEM data came from the United States Geological Survey (USGS 1999) archive on CD-ROM. Digital elevation data sometimes come with inherent error due to the resolution of the data or spatial limitations due to rounding of elevations in creating the dataset (USGS). In ARC/INFO’s GIS terminology, points that actually may not exist or points that create inaccurate surface flow conditions on a DEM are referred to as ‘sinks’ or ‘peaks’ depending on over or under estimating values. These depressions in DEM data were filled to ensure proper representation of water flow and watershed systems. In this study area, using the ‘FILL’ command in ARC/INFO 7.0, a total of 650 sinks were filled and 1720 peaks were leveled.

GIS programming allows for user-specified outlets in which the computer calculates the total contributing surface area using a series of complex algorithms. A script was downloaded from the Esri’s ArcView Hydrological Modeling help pages to delineate watersheds from the elevation data. It requires the use of Esri’s ArcView Spatial Analyst and Geoprocessing extensions, based on a point that is specified with a cursor in the view. The shape of the surface the elevation change at any given point determines the direction of which water will flow. The hydrological analysis tools available in GRID, help portray this type of natural system. The pour points or outlets specified in this study were determined with multiple stipulations. A need to compare prior studies in the area and update these findings as well as the priority identified by CERCLA.

Watershed boundaries are a key requirement for developing statistical analysis in hydrologic modeling. The entire study area was named the Patagonia Experimental Watershed, which is comprised of these 5 sub-watersheds. Table 1 describes the size of these (sub-) watersheds:
 

The watersheds that were effectively "poured" into the DEM were labeled according to their corresponding canyons (see Figure 2). The different thematic layers used to describe the environmental variables within this derived watershed boundary can be integrated as input into a lumped parameter model for predicting sediment loss and carried even further into quantification of sediment amounts and location of deposits.

Stream networks have been developed using the FLOWDIRECTION and FLOWACCUMULATION commands in ARC/INFO’s GRID module. The tools within GRID allowed for the determination of down slope water flow. This digital representation of likely channels was compared to known stream networks observed from a Digital Line Graph (DLG) created by the Arizona Land Resource Information System (ALRIS) within the Arizona State Land Department. This data was downloaded from the AZGENREF library, a digital library at the University of Arizona Advanced Resource Technology (UA/ART) Group. This indicated the derived stream channels were relatively accurate, see Figure 3. Once the stream channels and direction of overland flow were determined, these, along with the 30-meter DEM, were used in describing watershed characteristics and are primary features in the applied hydrological models.

Universal Soil Loss Equation

One application of a GIS is the prediction of spatial variability in surface erosion is the Universal Soil Loss Equation (USLE) which can be used to assess watershed conditions (Brooks, et al. 1997). The USLE is an empirical formula used to predict average annual soil loss in tons per acre per year (Wischmeier and Smith 1978). In this study, the location of sites with high potential erosion allows for identification of critical areas. The amount of erosion calculated on a cell by cell basis also acts as input to derive sediment delivery calculations in a second model. The USLE formula to calculate soil loss is as follows:

A = R * K * L * S * C * P

Where:

This formula is appropriately suited to be applied in a GRID based environment where map algebra can be performed. The analyses were broken up into separate tasks by data theme. These were determined after careful inventory and evaluation of the current digital database. The models chosen for use in this study defined the necessary objectives to be met using a GIS. Acquisition of the necessary factors is described below.

The R factor was the easiest data to retrieve as the Patagonia Experimental Watershed is small enough (only 42,306.5 acres) that it fell within one designated vicinity as defined by the USDA, Soil Conservation Service (1976) for areas of strong relief, to be a constant value of 80.

The K factor requires acquisition of accurate geo-spatial soils data. This proves to be the most difficult parameter to quantify in this study area as with most hydrological models requiring soils data. The USDA, NRCS, offers a couple of digital soils databases. The State Soil Geographic (STATSGO) Data Base is the only digital publication of an area within Santa Cruz County. This data was downloaded as Arc Info coverages and unzipped and untarred from its compressed deliverable. The projection was altered from Albers, NAD27 to UTM, Zone 12, NAD83.

The study area was digitized with the use of a Digital Raster Graphic (DRG) downloaded from the AZGENREF library, of the Nogales, AZ quadrangle, which included the Patagonia Experimental Watershed and the town of Patagonia. This was then used to subset the soils map extents to the watershed boundary. The frequency command was then issued and eight known soil types emerged as the composition of this area (see Figure 4). Each of the soils types, known as MUID’s (Map Unit Identifiers), were composed of multiple soils series (compnames). For example, MUID: AZ032, is composed of 55% Comoro, 25% Riveroad, and 20% Arizo type soil components. Each of these components in turn has many different soil descriptors assembled in a complex relational database. Of these were two items of interest. The first was represented by the term ‘kfact’, described as the soil erodibility factor, which includes rock fragments. The second was under the guise, ‘kffact’, which was described as the soil erodibility factor that was fragment free for use in the USLE. Unfortunately, this number was listed predominantly as zero, which when multiplied into the USLE would predict a zero amount of erosion, which is highly unlikely.

Therefore, both factors were interpreted for all components of the MUID’s and weighted averages were calculated for kfact and kffact accordingly for purposes of comparison. These values are available in Table 2:
 

The STATSGO data soils maps were created by interpretation of other, more detailed soils maps in the area and reasonable estimates were calculated. These were then digitized using USGS 1:250,000 scale, 1 by 2-degree quadrangle series maps for reference. Because of the idiosyncrasy within the soils database and due to the poor resolution that they were developed with, alternative data was needed for this study.

Preexisting high-resolution soils data could not be found in digital form. The most accurate soil information for the study area was in hard copy. This was available in the ‘Soil Survey of Santa Cruz and Parts of Cochise and Pima Counties, Arizona, 1979’ also a product of the USDA, NRCS, formerly the Soil Conservation Service and the Forest Service in cooperation with the Arizona Agricultural Experiment Station. This data was created according to the site conditions in 1971, when soil scientists drew the boundaries of identifies soils types onto aerial photographs. The scale at which these were published is 1:20,000.

The task of automating these soil maps was extremely tedious. The aerial photos had not been orthoganalized, and contained distortion. A total of 15 maps composed the study area as laid out 5 by 3 on a dining room table. These maps were scanned using an 8-bit black and white drum scanner at 100 dpi into TIFF format. The images were imported into ERDAS IMAGINE and the white borders were removed through subset decollaring processes. Five CD-ROM’s containing Digital Orthophoto Quarter Quads (DOQQ’s) were used to register and rectify the scanned soils maps. Known points were identified on the aerial photo and matched to points on the DOQQ’s, these were referred to as Ground Control Points (GCP’s). This was the most time consuming portion of this project as the aerial photos were taken some 30 years prior to the DOQQ’s and buildings, trees, and waterways had changed considerably. The easiest and most accurate objects to identify were roads and intersections of roads with other features. These appeared to have the same shape throughout time, although some forest roads were now out of use, or had been paved or widened.

A 3^rd order polynomial transformation requires a minimum of 10 GCP’s to be identified. However, the level of accuracy increases as more points are entered and widely distributed. The GCP prediction tool within ERDAS IMAGINE uses the current transformation parameters to guess where the user will locate GCP’s from the work in progress to source data, this enables the user to identify when enough points have been entered to ensure that the transformation is accurate (ERDAS 1997). An average of 80 GCP’s were identified on each aerial photo and cross-referenced with the source data for this study (see Figure 5). The cubic convolution method of resampling was performed to effectively pierce the aerial photo with pinpoints to known real time coordinates and stretch or fold the picture to accurate proportions. This sampling method is suggested for aerial photos in which the cell size is dramatically changed (ERDAS 1997).

This transformed the abstract piece of paper into an accurate representation of real time and space with registered known coordinates. The cubic convolution method resamples using an algorithm which recognizes the data files of 16 pixels in a 4 by 4 window, and this creates the most accurate output when ortho-rectifying aerial photos (ERDAS 1997). Error still exists despite the high number of GCP’s used to control the transformation. Error existed in the DOQQ’s and new error was introduced in the resampling process. However, the photos edge-matched positively and roads, rivers, trees and soil polygons merged together seamlessly when mosaiced to create the big picture. The raster geometric correction was successful for use in this project. The final .IMG file was converted and compressed within ARC/INFO to TIFF format and laid out onscreen with known vector coverages of digitized roads and rivers overlaid to check for accuracy and error. The most useful was the road coverage downloaded the AZGENREF library, which identified error to be within 0-40 meters. A small portion of this analysis is available in Figure 6. This seemed to be acceptable error for the project.

The soils data that had been inscribed on the aerial photos was then automated through the process of on-screen digitizing in ARCEDIT. The distance command identified acceptable tolerances, node snap to closest 100 meters and weed and grain tolerances to 15 meters. The user-friendly graphical user interface (GUI) called ARCTOOLS was employed for the initial digitizing (see Figure 7). The coverage was then cleaned manually using command line editing and topology was built.

User defined items were added to the newly digitized soil coverage feature attribute table to define the map unit descriptions: soil series, slope angle and previous erosion. Labels were created and attribution of the new soils coverage was completed utilizing the ARCEDIT GUI called ‘forms’. This allowed for regular segmentation of space allowing for 443 polygons to be attributed against the labeled polygons of the final aerial TIFF, as a backdrop. The image was then subset to the size of the study area, leaving 305 polygons which was displayed in ARCPLOT (see Figure 8) and is also available as a soil series map to interested others (see Figure 9).

The k factor values derived from the USDA, SCS Technical Notes in Phoenix Arizona, Sept. 1, 1976, were added as an item to this coverage’s attribute table. Thirty-four different soil types are represented in this area; some have two k factors depending on multiple associations or complexes. These are averaged to calculate the total k factor per soil type. These are listed in Table 3:
 

The projection was defined according to its origin and the Patagonia Experimental Watershed boundary was used to clip the extents of the new soil map to the area of interest, this was then converted to GRID format.

The watershed contains 73 known mine sites according to a vector point coverage downloaded from the AZGENREF library, created by the Bureau of Mines Mineral Availability System (MAS) dataset. These mines are relatively small, they are not visible even on a 1- meter resolution DOQQ, yet are very visible in the field and pose potential to be highly toxic to the watershed at hand. Therefore the point coverage was converted to a grid cell size 30 meter with k factor of 0.55, at these point locations. This was laid over the Soil Conservation grid of k factors to accurately represent erosion potential of the soils due to mine, dump, and tailing pile presence. This was the final K factor grid used in the USLE equation. A comparison of the three derived K factor values for the study area is available in Figure 10. High resolution proved to be critical for this portion of the study.

The S factor is very closely associated with the L factor. The S is the slope gradient factor and the L is the length of that slope. The slope was calculated from the 30 meter DEM discussed earlier. This was calculated in percentrise in order to fit into the equation properly. This percentrise was then plugged into the formula:

The USLE was created to predict soil erosion delivered to the base of a 22 meter agricultural plot (Wischmeier, 1976). In this study, each 30 meter cell’s flow length was calculated as & = 96.24 feet and plugged into the following formula:

L = (&/ 72.6)^m, where m= 0.5

The S and L factors are then combined to form the LS factors using the following formula:

The C factor is the cropping or vegetation management factor. This is a user-defined number assigned according to vegetation type. C values are derived from the SCS Technical Notes (table 6), for permanent pasture, rangeland, and idle land according to a vegetation coverage of the study are that was downloaded and clipped form the USGS GAP Analysis Vegetation and Land Cover geo-spatial data-set at the AZGENREF library. The watershed comprised 10 vegetation types, predominantly Encinal Mixed Oak and Semidesert Mixed Grass (see Figure 11). The C factor to vegetation type is available in Table 5:
 

The P factor, or conservation practice factor was not relevant to the study area and therefore calculated as value = 1, which does not negatively or positively influence the output of the model.

Once all of the factors are accounted for in GRID based environments, with anything outside the Patagonia Experimental Watershed equaling NODATA, the grids could be multiplied together to get a gross estimate of annual potential soil loss in tons per acre per year. This soil-loss map is available in Figure 12. The total erosion in Patagonia Experimental Watershed per year is 14,129,173,075 tons. The sum of total potential soil loss per watershed is described in Table 6:
 

In watersheds this large in size, most sediment gets deposited within the watershed and only a fraction of soil that is eroded will reach the stream system or watershed outlet. The results derived from the USLE are used for planning purposes to predict the impact of land use on soil erosion and to identify sensitive areas. The determination of areas with potentially low erosion rates is useful if the mitigation strategy is to physically move the toxic materials to sights of low potential harm. It also identifies the critical source areas of pollutants.

SEDMOD

SEDMOD, an acronym for Spatially Explicit Delivery MODel, calculates a Sediment Delivery Ratio (SDR) which can be useful to calculate that amount of eroded material that would be available for transport and that is deposited along hillslopes and streams (Fraser 1999). A SDR is not homogenous across a watershed, instead it varies with changes in watershed area and slope (Ostercamp and Toy, 1997). SEDMOD allows for the calculation of this spatial variation utilizing a GIS. The SDR is multiplied with the predicted amount of erodable soil to calculate Non-Point Source (NPS) sources and sinks within a watershed.

The delivery ratio evaluates deposition that occurs in overland flow before reaching the stream channels (Haan, et al. 1981). Many factors are addressed when calculating this ratio: water availability, texture of eroded material, ground cover, slope shape, gradient and length, surface roughness, and other on-site factors according to the Stiff diagram (Forest Service, 1980). SEDMOD incorporates these parameters in a cell by cell calculation of uniquely specific derivations for changes over space.

Input

Grids representing terrain, soil type, land classification, and soil loss were created at a 30-meter resolution for input into SEDMOD according to specifications (Fraser, 1999). These grids were DEM, SOIL_TEXTURE, ROUGHNESS, and SOIL_LOSS (see Figure 13). SEDMOD also calls for the optional input of SOIL_TRANS (saturated soil transmissitivity) and STREAM (stream network) grids, which were not included in this project. The All the input grids were clipped to the watershed area leaving the surrounding non-watershed cells with a value of NODATA.

The SOIL_TEXTURE grid was created using the STATSGO soil database calculating percent clay. The newly created soils coverage was originally intended for this study, but the documentation of clay percentage was variably unavailable. Within the STATSGO database, table names layer contained two related components: ‘clayl’ and ‘clayh’. These are the minimum and maximum values for the range in clay content for the topsoil layer, expressed as a percentage of the material less than 2 mm. in size (USDA 1995). This was extrapolated in Microsoft Access, averaged in Microsoft Excel, and incorporated into the already existing STATSGO MUID ARC vector polygon coverage feature attribute table and rasterized to GRID format. This data is available in Table 7:
 

The ROUGHNESS grid was derived from of the previously acquired GAP Analysis vegetation and land use data. These vegetation descriptions were used to estimate Manning’s Surface Roughness Coefficients for Overland Flow (Fraser, 1999).

These input grids were then shuttled through a series of Arc Macro Language (AML) scripts in SEDMOD that enable a friendly GUI for calculation, display, and analysis. Secondary grids were derived from the DEM, very similar to the hydrologic modeling described previously (flow accumulation, direction and channel networks). Finally the six grids directly used for calculation of SDR were made. These are available in Figure 14. These are a combination of DEM and site description grids.

Output

The SDR was calculated using SEDMOD (Figure 15). This ranged from 0 –79% of eroded material to be transported under the process of overland flow. This was then multiplied by the soil loss equation prediction to derive net sediment or nonpoint source pollution delivery to the watershed (Figure16) and to calculate the net sediment delivered to the stream channels and the bordering riparian area (Figure 17).

SEDMOD also calculates watershed total potential gross erosion (316,220.7 tons/ acre/ year), estimated sediment delivered to the streams (51,500.1 tons/ year) and finally estimated total delivery to the outlet (16,347 tons/ year).

Analysis

This type of data can be useful in many applications. In this study, critical nonpoint source pollution sources are the focus of analysis. A water quality chemical analysis, including a rigorous water-sampling project has been done in conjunction with the USGS, under the guidance of Floyd Gray, geologist at the Southwest Field Office, in Tucson, Arizona. Results of the aforementioned will be used in comparison with this watershed delineation to complete a paired watershed analysis in the area in an attempt to pinpoint and delineate nonpoint sources of pollution.

The amount of contributing sediment from previously mined areas was analyzed and compared to the pH measured at water sample points along these stream networks to create a synthesis of the mine contribution of sediment and that relationship to acidic conditions.

The coverage of 25 water sample points was converted to grid format. Each sample locality was then entered as a watershed outlet or pour point in ArcView to calculate total contributing drainage area to that point and converted to a 30-meter resolution raster grid. These 25 grids were vectorized and the polygon coverages were used to clip the extents of the mines within the contributing drainage basin (see Figure 18). The sample point locations were queried to discover the amount of sediment (NPS pollutant) delivered to each point according to the derived riparian delivery map created by SEDMOD. This data is available in Table 8:
 

The name and commodity of the mines were teased out using the sample command in GRID. Each mine location grid cell was queried to see the amount of contributing sediment it was supplying to the stream system. These numbers will be calculated to percentages of the sample points’ total sediment. The total mine contribution of sediment to be delivered is available in Table 9:
 

These percentages will then be compared to the change in pH at these points.

The characteristics of acid mine drainage has been described abundantly in this area (Gray 2000, Chaffee 1981, Dean 1982, Chatman 1994, Hyde 1995), however, this study is the first to test the correlation of sediment from mines and acidity in the water. It also pinpoints major contributors of NPS pollutants to the system and identifies areas of high pollution or NPS sinks.

Conclusions

This study offers a lot of advanced technological applications combining a common erosion prediction model with a new approach to spatially derived sediment delivery. GIS is a great platform for this type of study, it facilitates a new approach to watershed management, allowing for a spatial analyses. It limits the mathematical computation for the user, yet equates multiple complex algorithms with the software described. The models used for this study were fairly simple to manipulate within a GIS once all of the necessary input layers had been created and/or acquired.

The products of this study are useful for the project described and can be stored within a digital GIS database for other studies. Other hydrologic models can also be easily applied to this dataset to achieve new goals and management practices. This type of data facilitates management decisions, by creating models of hypothetical situations, money and time can be saved and mistakes can be avoided.

Future research recommendations would include the routing of chemical constituents through the watershed system to the sampling points to account for each mines’ chemical contribution. It would very useful to compare the amount of sediment loaded from mines to the amount of toxic chemicals such as arsenic found per water sample. Additional water sampling points would also help to better characterize the watershed, especially in Providencia Canyon. Finally, high-resolution aerial photos could help to better identify vegetative conditions (at an Association level), which would greatly influence the model outputs.

References

Brooks, K.N., P. F. Ffolliett, H. M. Gregerson and J. L. Thames. 1997. Hydrology and the Management of Watersheds, 2^nd Edition. Iowa State University Press, Ames, Iowa: 438- 439

Chaffee, M.A. Hill, R.H., Sutley, I.J., and Wastterson, J.R., 1981, Regional Geochemical studies in the Patagonia Mountains, Santa Cruz County, Arizona: J. Geochemical Exploration, v. 14, 135-153

Chatman, Mark L., 1994, Mineral Appraisal of Coronado National forest, part 7:US Bureau of Mines Mineral Land Assessment Open File Report 22-94, 6 plates: 117

Dean, Sheila, A.,1982, Acid drainage from abandoned metal mines in the Patagonia Mountains of southern Arizona, Coronado national forest: USDA Forest Service

ERDAS, Inc., 1997. ERDAS Field Guide, 4^th Ed. Atlanta, GA: 307- 342

Forest Service, 1980. ‘An Approach to Water Resources Evaluation of Non-Point Silvicultural Sources’, a procedural handbook. EPA-600/8-80-012, Environmental protection Agency, Washington, D.C.

Fraser, Robert Hungerford, 1999. SEDMOD: A GIS-Based Delivery Model for Diffuse Source Pollutants. Dissertation. Yale University. Connecticut. (p. 76-99).

Gray, Floyd, Laurie Wirt, Kerry Caruthers, Ailiang Gu, Carlos Velez, Douglas M. Hirschberg, and Karen Bolm, 2000. "Transport and Fate of Zn and Cu rich, Low pH Water from Porphyry Copper and Related Deposits, Basin and Range Province, Southeastern Arizona" in Proceedings of the Fifth International Conference on Acid Rock Drainage, May 21-24, 2000, Denver, CO: (1) 411 - 430

Haan, C. T., B. J. Barfield, J. C. Hayes, 1981. Design Hydrology and Sedimentology for Small Catchments. Academic press, Inc, CA: 238-310

Hyde, Peter, 1995, Water Quality at the Trench Camp, Patagonia Mountains, Santa Cruz County, Arizona – April 3-4, 1995: Arizona Department of Environmental Quality Report, Aquifer Protection section: 1-15

Osterkamp, W. R. and T. J. Toy, 1997. Geomorphic Considerations for Erosion Prediction. Environmental Geology 29 (3/4): 152- 157

USDA. U.S. Forest Service, 1977. Anatomy of a Mine from prospect to Production, USDA Forest Service. General Technical Report INT- 35. Intermountain Forest and Range Experiment station Forest Service. US Department of Agriculture Ogden, UT. June 1977

USDA. Natural Resources Conservation Service, 1995. State Soil Geographic (STATSGO) Data Base Data Use Information. USDA NRCS, National Soil Survey Center. Miscellaneous Publication Number 1492.

USDA. SCS, 1976. SCS Technical Notes in Phoenix Arizona, Sept. 1, 1976

USGS. "USGS_DEM." 13 Sept 1999, < http://edcwww.cr.usgs.gov/glis/hyper/guide/usgs_dem (1 June 2000)

Wischmeier, W. H., and D. D. Smith, 1978. Predicting Rainfall Erosion: A Guide to Conservation Planning. Agronomy Handbook No. 537, U.S. Department of Agriculture, Washington, D.C.

Laura M. Brady

Research Assistant

University of Arizona/ U.S. Geological Survey

520 N. Park Ave.

Tucson, AZ 85705

lmbrady@u.arizona.edu