GIS analysis of spatial variability of contaminated watershed components in a historically mined region, Arizona

By Laura M. Brady^[1], Floyd Gray², Craig A. Wissler³, and D. Phillip Guertin³

 

 

 

2000

 

 

 

 

This report is preliminary and has not been reviewed for conformity with U.S. Geological Survey editorial standards or with the North American Stratigraphic Code.  Any use of trade, product or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

 

U. S. DEPARTMENT OF THE INTERIOR

U.S. GEOLOGICAL SURVEY

  Abstract:

GIS can be used to integrate and accurately map field studies, information from remotely sensed data, watershed models, and monitor the dispersion of pyrite rich mine waste and tailings.  The study area included historic small to medium scale mining in the Patagonia and southern Santa Rita Mountains District, 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, most of the area is in the Coronado National Forest. 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 possible ways the pattern could be used to develop mitigation strategies are suggested.

Introduction

Patagonia and the southern Santa Rita Mountains area, located in southern Arizona (see fig.1) was mined intermittently from the 1600’s to the mid- 1960’s primarily for silver, lead, copper and zinc as well as gold. What is the movement of water through these mined areas?  What factors effect that flow? Several of the mines have been identified as 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 nearby resident 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 (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 used 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 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 (see USGS website "USGS_DEM." 13 Sept 1999, http://edcwww.cr.usgs.gov/glis/hyper/guide/usgs_dem (downloaded 1 June 2000)). 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.

Table 1: Delineated sub-watershed size
 

 

The watersheds that were effectively "poured" into the DEM were labeled according to their corresponding canyons (see fig.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 commands in ARC/INFO’s GRID module (the FLOWDIRECTION and FLOWACCUMULATION). 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 as shown in  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 and others, 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:

A = annual soil loss in tons per acre per year

R = rainfall erositivity factor

K = soil erodibility factor

L = slope length factor

S = slope gradient factor

C = cover management factor

P = erosion control practice factor

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,307 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, Ariz. 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 fig. 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 percent Comoro, 25 percent Riveroad, and 20 percent 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.

Table 2: STATSGO derived K factors.
 

 

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 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 fig.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 fig.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 fig.8) and is also available as a soil series map to interested others (see fig.9).

Table 3: SCS Technical Notes K Factor values.
 

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:

S = (0.43 + 0.30s + 0.043s^2)/ 6.613

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:

LS = L * S (10,000/ (10,000 + s^2))

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 fig.11).

Table 4: Assigned C Factor to vegetation type.