Gully Erosion of Archaeological Sites in Grand Canyon:

Gully Erosion of Archaeological Sites in Grand Canyon:
The Question of Causality

Paul A. Petersen
Geology Department
Utah State University
CEE 6440

1. Background

2. Objective

3. Study area

4. Data sources

5. Processing

6. Results

7. Discussion

8. Conclusions

9. Future work

10. Acknowledgements

11. References

Background

Gully erosion along the Colorado River corridor in Grand Canyon National Park has increased in magnitude and frequency in the last few decades, and as a result many irreplaceable archaeological sites that lie upon and within Holocene terraces are being damaged. Hereford et al. (1993) were the first researchers to suggest that Glen Canyon Dam operation was a potential factor in the increased gully erosion in Grand Canyon. They hypothesized that the effects of flooding at the level of the pre-dam terrace slowed or prevented terraced-based gullies from reaching the river by depositing sand at the mouths of the channels. The lack of flooding and deposition during post-dam time allowed previously terrace-based gullies to reach the Colorado River, reducing effective baselevel by 3 to 4 meters. They claimed that an unusual run of high rainfall and runoff years drove the erosional process, but baselevel controlled the depth of erosion and knickpoint retreat throughout the catchment, implying that gully incision will be intensified relative to pre-dam conditions until the channels adjust to their new, lower baselevel. The geomorphic community (J. Pederson, J. Schmidt, personal communication, 2002), has heavily criticized Hereford et al.’s (1993) baselevel hypothesis. In particular, the baselevel hypothesis is thought to be flawed conceptually. The baselevel hypothesis has not yet been related to the broader geomorphic literature, nor has it been successfully verified or falsified by data. Climate fluctuations and upcatchment controls serve as viable, yet untested controls in the increase of gully erosion in Grand Canyon. My Master’s thesis will place the baselevel hypothesis within the context of previous geomorphic literature and test catchment properties as a primary control for concentrated degradation.

Conceptual setting of gully erosion of Holocene terraces in Grand Canyon (Thompson and Potochnik, 2000). Gullies of interest are eroding through sandy terraces and feature relatively impermeable talus or bedrock upper catchments.

Objectives (top)

The goals of this GIS class exercise are to

1) demonstrate upcatchment control of gully position and extent using an empirical area-slope threshold;

2) model change in erosion potential due to vegetation change

Study Area (top)

Area-slope measurements were taken at 29 gully heads from 7 sites clustered in eastern and western Grand Canyon. High-resolution terrain models were created from 4 photogrammetry sites in western Grand Canyon.

Grand Canyon study sites. All modeling was performed on the western study sites.

The gullies of interest are relatively small, ranging from 20 to 150 meters of thalweg length. Channel widths vary anywhere from 15 cm up to 2 meters. Catchments vary from 0.002 to 0.2 hectares. This small scale of gullies calls for a high resolution of analysis.

Data Sources (top)

Elevation data/Terrain models

Slope and contributing drainage area above each gully head were surveyed in the field during the February, 2002 trip. Also in February, stereopairs were photographed at the four western sites from a helicopter 800 feet above the surface. These black and white photos were then scanned at 12.5 microns, resulting in digital images with 1.98 cm pixels. Photogrammetry was performed on these four stereopair sets, extracting up to 150,000 elevation points for each site using the automatic collection method. These photogrammetry points were edited and merged with survey points of the gully channels to create a robust terrain model. Point density and gully scale warranted interpolation of 10-cm DEMs, using the spline tension technique.

Name	Point distribution	Hillshade
Indian Canyon	Indian_pt	Indian_shade
Arroyo Grande	Arroyo_pt	Arroyo_shade
Granite Park	Granite_pt	Granite_shade
Gorilla Camp	Gorilla_pt	Gorilla_shade

Internet links to elevation point and hillshade maps of four photogrammetry sites.

Vegetation and infiltration

During the October, 2002 research trip, several vegetation transects were conducted at each site using the 8-pin point count method. This estimated percent ground cover for the 5 dominant cover types (rock, bare ground, shrub, grass/forb, cryptobiotic crust) of a given site. Along each transect, an infiltration test was performed at each type of cover (excluding rock), using a tension disk infiltrometer. In all, 59 infiltration tests were performed on bare ground, 34 beneath shrubs, 31 on cryptobiotic crust, and 30 on grass. Since saturated hydraulic conductivity (K_s) represents minimum infiltration rate, these data were then converted to saturated hydraulic conductivity using the method outlined in Reynolds and Elrick (1991).

Procedure (top)

Area-slope thresholds

Contributing drainage area and up-catchment slope of each gully head were plotted in log-log space, yielding an empirical erosion threshold line described by a power function, below which theoretically no incision occurs (Vandaele et al., 1996).

The power equation that defines the erosion threshold in Grand Canyon is y = 0.0137x^-0.4671. Since slope is used as a proxy for runoff, this equation means that for any given contributing drainage area (x), there is a critical slope (y) that must be exceeded to entrain sediment and cause erosion. Above the threshold line, slope is exceeded, indicating that the area is sensitive to gully erosion. The converse is theorized to be true below the threshold line. Thus, by applying a threshold equation to a terrain model in a GIS, one can show areas that are prone to gully erosion. To do this, I used Taudem to calculate slope and contributing drainage area grid of each of the four DEMs. The D¥ algorithm (Tarboton, 1997) was used to best model sheetflow on a hillslope. Gridmath was performed to multiply the area grid by 0.0137 and take it to the –0.4671 power, creating the critical slope grid. Finally, a simple query was run, asking the GIS identify all cells in which the slope (S) exceeds the critical slope (S_cr). Cells that are “true” are considered to be sensitive to gully erosion.

Conceptual diagram showing cells in which the actually slope exceeds the critical slope necessary to cause gullying.

The accuracy of this simple prediction was tested during the October, 2002 field trip by mapping all observed gully features, and comparing their position and extent to the threshold prediction maps.

Erosion index based on ground cover

To keep within a reasonable scope, only one of the four sites, Gorilla Camp, was selected to model changes in erosion potential due to vegetation change. Since cover type strongly influences a given soil’s ability to take in water we analyzed infiltration rates in terms of cover type.

Cover type	Bare ground	Shrub	Crypto	Grass
N	59	34	31	30
Average	0.0197	0.0042	0.0142	0.0079
Median	*0.0077*	*0.0027*	*0.0072*	*0.0037*
Stdev	0.0300	0.0046	0.0144	0.0104
Min	0.0002	0.0000	0.0005	0.003
Max	0.2896	0.0206	0.0619	0.0522
Histogram	bare_k	shrub_k	crypto_k	grass_k

Ground cover at Gorilla Camp was quantified through six vegetation transects: 1) and 2) talus; 3) western gully heads; 4) western gully mouths; 5) eastern gully heads; 6) eastern gully mouths.

	Talus %	W heads %	W mouths %	E heads %	E mouths %
rock	76.25	28.75	32.28	0.00	11.47
bare	0.00	5.63	10.13	68.75	29.30
crypto	20.00	49.37	32.28	23.75	35.03
grass	0.00	0.00	0.63	0.63	1.91
shrub	3.75	16.25	24.68	6.88	22.29

Percent ground cover, data from vegetation transect point counts.

Each cover type was assigned an identification number (rock = 1; bare = 2; crypto = 3; grass = 4; shrub = 5). Five separate grids were clipped from the Gorilla Camp DEM in order to represent and model each transect zone separately.

Vegetation transect endpoints and polygon clip coverages for vegetation grids.

Discrete random number generation was performed in Excel, in order to stochastically represent the percent cover type for each transect in a grid. Cover type grids for each zone were then assembled and merged back together to form one cover type grid for the entire site. Next, each cover type ID, 1 through 5, was assigned the median K value from the population sampled in the field (ie., rock = 0 cm/s; bare ground = 0.0077 cm/s; crypto = 0.0072 cm/s; grass = 0.0037 cm/s; shrub = 0.0027 cm/s).

Hydraulic conductivity map of Gorilla Camp, based on ground cover.

The goal of this exercise was to achieve a relative shear stress for each 10-cm cell in the site. On a hillslope shear stress is primarily a product of slope and water depth. I calculated a water depth for each cell by subtracting the K grid from 0.008, essentially simulating the infiltration-excess of 0.008 cm of rain falling in 1 second. This is a very high intensity, but was necessary to avoid negative numbers in the calculation. Again, the goal of this model was to find an erosion index, not an actual shear stress value (why g and g were left out). The result of subtracting the K grid from 0.008 was the amount in centimeters of water that did not infiltrate in each 10-cm cell. The amount of infiltration excess was then accumulated by using the D¥ flowaccumulation algorithm with the infiltration excess grid as a weight. Finally, the runoff accumulation grid (“depth” for each cell) was multiplied by the slope grid to determine relative shear stress of the event on each cell.

Erosion index map of Gorilla Camp, delineating gully and wash channels. Units are meters (water depth x slope)

In order to model shear stress change at the site due to a vegetation change, the input vegetation grid was altered to increase the percentage of grass cover at the expense of shrub cover.

	Talus %	W heads %	W mouths %	E heads %	E mouths %
rock	76.25	28.75	32.28	0.00	11.47
bare	0.00	5.63	10.13	68.75	29.30
crypto	20.00	49.37	32.28	23.75	35.03
grass	3.13	13.75	21.52	6.25	21.02
shrub	0.63	2.50	3.80	1.25	3.18

Altered percent ground cover, changing only proportions of grass to shrub.

The aforementioned procedures were then performed with the same parameters on the new cover type grid, resulting in a new relative shear stress map. The altered shear stress grid was subtracted from the “natural” shear stress grid, resulting in a difference map showing the results of a vegetation change from shrub to grass on erosion potential.

Results (top)

Gully prediction using area-slope threshold

	# Gullies Observed	# Gullies Predicted	Percentage
Indian Canyon	10	10	100%
Arroyo Grande	19	18	95%
Granite Park	19	18	95%
Gorilla Camp	20	20	100%
Total of 4 sites	68	66	97%

Summary of area-slope prediction analysis. If the position of an observed gully feature was represented on the prediction map, the gully was counted as “predicted.” In the following figures, polygons are archaeological features, yellow lines are cells that exceed the threshold, red lines are observed gullies digitized in.

Threshold map for Indian Canyon

Threshold map for Arroyo Grande

Threshold map for Granite Park

Threshold map for Gorilla Camp

Change in erosion potential due to change in vegetation

Difference map of Gorilla Camp showing the net change of relative shear stress due to a decrease in shrub cover and proportional increase in grass cover.

Shear stress	Shrub:Grass = high	Shrub:Grass = low	Difference map
Average	1.4844	1.4750	0.0095
Stdev	37.4230	37.2017	0.2410
Min	0.0000	0.0000	-2.2826
Max	19302.5723	19162.3359	140.2363

Erosion index and difference map statistics. Erosion index is water depth x slope, so units are technically in meters, although it should be thought of unitless for this exercise.

Discussion (top)

Using an area-slope threshold to predict gully position and extent was successful in Grand Canyon, regardless of effective baselevel. This indicates:

1) The idea of using thresholds for such a prediction analysis is valid in Grand Canyon

2) The high-resolution DEMs created from photogrammetry are topographically correct

3) D¥ is successful at 10-cm scale resolution at showing where sheetflow concentrates on a hillslope

4) 10-cm DEMs are sufficient to predict decimeter-scale gully channels

5) Upcatchment processes can explain the initiation, position, and headward growth of gullies in Grand Canyon. Predicted gullies vary immensely in their order, scale, and effective baselevel, indicating that their causality is controlled by one common element: a sufficient amount of water running over a critically-steep slope. A baselevel change is not need to invoke gully initiation and expansion.

This empirical erosion threshold should be successful at predicting gully position and extent for any site in Grand Canyon with a similar geomorphic setting. The threshold equation cannot be used with confidence in geomorphic and/or climatic settings different from Grand Canyon alluvial terraces.

It should be noted that all sites contain many cells that predict erosion but do not feature gullies (see result maps). This phenomenon has several possible explanations:

It is better to think of these as cells sensitive to gully erosion. If a large enough storm initiates runoff, these cells indicate where the next new gully would erode. Similarly, if the landscape were disturbed through trailing, clearing vegetation, or climate change, new gullies would incise in these predicted areas.
The cells may be areas of deviant substrate or geomorphic setting relative to what typical gullies in Grand Canyon erode through. For example, a dune may have the right slope and drainage area combination to be sensitive to gully erosion in the model, but in reality infiltration is too high to ever sustain runoff. Similarly, a talus slope may have the topography to have predicted gullies, but the bouldery substrate requires too high of a shear stress to be moved by runoff that erodes sand-based gullies.
There may be a flaw in the terrain model itself. There are occasionally anomalous elevation points derived from the photogrammetry or ground survey. These can easily be identified by contouring the site, and then edited out during quality control.

Possible sources for “overprediction.” Artifact scatters are represented by orange polygons.

Keeping these factors in mind will help discern the true erosion potential of a site, and help identify areas most sensitive to landscape disturbance. This will be especially helpful in the fragile desert landscape of Grand Canyon National Park, as managers and archaeologists will be better able to monitor and protect eroding archaeological sites from both gully erosion and human impacts.

Montgomery and Dietrick (1994) note that channel heads are dynamic and will migrate up or down slope in response to changes in the erosion threshold. An increase of the channel initiation threshold can be caused by changes of physical properties and processes through both space and time, and will result in channel contraction. The converse is true if there is a decrease in the channel initiation threshold. One such physical property that could change the threshold is the soil’s infiltration capacity. This ability to take in water is not static either, and can change through time. For example, if a grassland region undergoes climate change and shifts toward shrub communities, infiltration capacity, through both changes in K and roughness, will decrease (Abrahams et al., 1995).

In the groundcover model presented in this project, I attempted to show the increase in runoff described by Abrahams et al. (1995). The model is crude and only takes into account K, but is successful in demonstrating the change in infiltration, runoff, and shear stress in a given site for a given rainfall intensity. Cells consistently experienced less shear stress in the grass-dominated model than in the shrub-dominated model. Essentially this represents an increase in the erosion threshold, and would result in channel infilling and contraction.

In this light, if upcatchment processes control position and extent of gullies in Grand Canyon, then changes in these processes should account for increased erosion rates and network rejuvenation throughout time. Climate does change on the Colorado Plateau on a decadal scale, as Hereford and Webb (1992) showed that the period between 1940 and the late 1970s was an exceptionally dry time, whereas 1979 to the present has been relatively wet. In other words, climate can be the driving force for this particular gullying, especially if coupled with changes in vegetation distribution. In particular, a long period of drought (and subsequent vegetation shift toward shrubs) followed by a very wet period could result in high runoff and erosion (Balling and Wells, 1990). Rogers and Schumm (1991) showed that only a small change in vegetative cover is needed to have an effect on erosion and sediment yield. The shear stress model in this project crudely demonstrates this concept, but much more research needs to be done to validate such a hypothesis in Grand Canyon.

Conclusions (top)

Funding: Grand Canyon Monitoring and Research Center; Geological Society of America; Colorado Scientific Society

Field, lab, and thought assistance: Joel Pederson, David Chandler, Wally McFarlane, Jen Dierker, Jay Norton, Stacy Petersen, Sammie McFarlane, Isaac Larsen, Jesse Allen, Lynn, Thomas, Scott Cragun, Jack Schmidt, Tom Monaco

Abrahams, A.D., Parsons, A.J., and Wainwright, J., 1995, Effects of vegetation change on interrill runoff and erosion, Walnut Gulch, southern Arizona: Geomorphology, v. 13, p. 37-48.

Balling, R.C. and Wells, S.G., 1990, Historical rainfall patterns and arroyo activity within the Zuni River drainage basin, New Mexico, Annals of the Association of American Geographers, v. 80, n. 4, p. 603-617.

Hereford, R., Fairly, H.C., Thompson, K.S., and Balsom, J.R., 1993, Surficial Geology, Geomorphology, and Erosion of Archaeological Sties along the Colorado River, Eastern Grand Canyon, Grand Canyon National Park, Arizona: U.S. Geological Survey Open-File Report 93-517.

Hereford, R. and Webb, R.H., 1992, Historic variation in warm season rainfall on the Colorado Plateau, USA: Climatic Change, v. 22, p. 239-256.

Montgomery, D.R. and Dietrich, W.E., 1994, Landscape dissection and drainage area-slope thresholds, in Kirkby, M.J., ed., Process Models and Theoretical Geomorphology, John Wiley & Sons Ltd, p. 221-246.

Reynolds, W.D. and Elrick, D.E., 1991, Determination of hydraulic conductivity using a tension infiltrometer: Soil Science Society of America Journal, v. 55, n. 3, p. 633-639.

Rogers, R.D. and Schumm, S.A., 1991, The effect of sparse vegetative cover on erosion and sediment yield: Journal of Hydrology, v. 123, p. 19-24.

Tarboton, D., 1997, A new method for the determination of flow directions and contributing areas in grid digital elevation models: Water Resources
Research, v. 33, n. 2, p. 309-319.

Thompson K.S. and Potochnik, A.R., eds., 2000, Development of a Geomorphic Model to Predict Erosion of pre-Dam Colorado River Terraces Containing Archaeological Resources: Cultural Resources Report No. 99-257 prepared for the Grand Canyon Monitoring and Research Center by SWCA, Inc., Environmental Consultants, Flagstaff, AZ, 130 p. plus appendices.

Vandaele, K., Poesen, J., Gover, G., and van Wesemael, B., 1996, Geomorphic threshold conditions for ephemeral gully incision: Geomorphology, v. 16, p. 161-173.

Table of Contents