San Juan River Historical Ecology Assessment Changes in Channel Characteristics and Riparian Vegetation September 2015 U.S. Bureau of Reclamation Agreement Number R09AP40004 Prepared by Steven Bassett The Nature Conservancy in New Mexico 212 East Marcy Street, Suite 200 Santa Fe, NM 87501 Introduction The San Juan River is the third largest tributary of the Colorado River. The river’s natural flow regime is modified by Navajo Dam and other smaller dams and diversions. The lower reaches of the San Juan River are submerged under Lake Powell, the reservoir behind Glen Canyon Dam. Portions of the lower San Juan River are classified as critical habitat for the Colorado pikeminnow (Ptychocheilus lucius) and razorback sucker (Xyrauchen texanus). These species were listed as endangered by the U.S. Fish and Wildlife Service (USFWS) in 1967 and 1991 respectively. In 1992 the San Juan River Basin Recovery Implementation Program (SJRIP) was established to recover the Colorado pikeminnow and the razorback sucker while allowing water development and other river management activities to continue. Actions pursued by the SJRIP include stocking the endangered fish, removing non-native fish that compete with and prey upon the native fish, providing fish passage around diversion dams, providing flows that mimic the natural hydrograph, installing screens and weirs to reduce entrainment, and more recently, modifying habitat by removing non-native vegetation and opening up secondary channels and backwaters that have been overgrown with non-native vegetation. Backwaters, embayments, and secondary channels provide the low-velocity habitats that are used by larval and small-bodied fish, including Colorado pikeminnow and razorback sucker, and are important in the recruitment process (Farrington et al. 2014). Loss of low-velocity habitat is thought to have contributed to the decline of the listed species (Holden 1999). Simplification of the river channel reduces the area of low-velocity habitats (Van Steeter & Pitlick 1998). Channel simplification can have many causes but in the study area is thought to be primarily caused by changes in hydrologic and sediment regimes due to regulation by Navajo Dam, and subsequent encroachment and bank armoring by non-native vegetation, specifically Russian olive (Elaeagnus angustifolia) (Holden 1999) and to a lesser extent saltcedar (Tamarix spp.) (Brooks et al. 2000). My study quantifies channel simplification and other changes in river planform morphology and stream bank vegetation using historical aerial photographs. Imagery from the 1930s, 1970s, and 2010s was used to map the active channel and bankside vegetation in each period. Digitized features were then compared to quantify changes. Changes in planform morphology and riparian vegetation likely occurred before the mid-1930s but these antecedent changes occurred before aerial imagery is available so are not included in this analysis. Study Area The San Juan River originates in the San Juan Mountains and flows east to west through Colorado, New Mexico, and Utah. The highest elevations of the basin are dominated by spruce-fir, mixed conifer and ponderosa pine forest but the majority of the basin is desert grasslands. The basin is 2 underlain with Mancos shale and oil and gas development is widespread. Cities in the basin include Farmington, Bloomfield, and Shiprock, New Mexico, and Durango, Colorado (Figure 1). Surface water development includes over 20 agricultural, municipal, and industrial diversions on the main stem and tributaries (NMOSE 2015). Navajo Dam was completed in 1962 as part of the United States Bureau of Reclamation’s Colorado River Storage Project and provides irrigation for about 70,000 acres of farmland (Glaser 1998) and 30 MW of hydroelectricity to the City of Farmington (EMNRD 2005). Geomorphological reaches have been delineated below Navajo Dam (Bliesner & Lamarra 2000). These reaches provide a common reference for research and facilitate comparison of studies. Reach delineations were produced by conducting multivariate analysis of geomorphic characteristics to identify areas with similar characteristics. Reaches 1 and 2 are below the confluence with Chinle Creek and are canyon bound. Reach 1 is heavily influenced by fluctuations in the surface elevation of Lake Powell. Reach 3, between Chinle Creek and Aneth, Utah, has a low gradient and is sand- dominated. Reach 4 is transitional between sand-dominated lower reaches and cobble-dominated upper reaches. Reach 5 is characterized by a wide valley and is primarily cobble. Reach 6 has extensive anthropogenic modification including four diversion dams. Reaches 7 and 8 are heavily influenced by releases from Navajo Dam and banks have been stabilized to control lateral movement of the channel (Bliesner & Lamarra 2000). The reaches included in this study were 3, 4, 5, and 6. Figure 1. Map of the study area. 3 Methods Analysis of historical aerial photographs is a common technique for quantifying landscape-scale habitat changes (Bliesner & Lamarra 2000; Cadol et al. 2010; Jones et al. 2011; Block 2014). Since the 1930s, aerial photographs have been collected regularly over the continental United States. In the San Juan Basin, this imagery is available beginning with black and white photographs taken in 1934, with imagery available for most of the study reach at least once a decade since. These photographs are stored in various archives and libraries but recently much of this imagery has been scanned and published electronically. For this study, imagery from 1934–1935, 1975–1979, and 2011 was analyzed (Table 1). Additional imagery is available for many other years but was not included in this analysis. The aerial photographs used in this analysis were selected because they offer complete coverage of the study area and allow documentation of pre-dam, post-dam, and current conditions. Table 1. Aerial imagery used in this analysis. Year Series Extent (VM) Description Scale 1934 Soil Conservation Service (via USFWS) 131–150 B&W1 1:30,000 1935 Soil Conservation Service(via UNM EDAC2) 187–249 B&W 1:31,680 1935 NAV 1935 (UTDNR3) 148–178 B&W 1:31,680 1975 USGS series VDXK0 184–235 B&W 1:30,800 1978 USGS series VEPO0 281–249 B&W 1:32,000 1979 USGS series VETD0 131–189 B&W 1:24,000 2011 USDA National Agriculture Imagery Program 131–249 Color + IR4 1m 1 Black and white panchromatic imagery 2 University of New Mexico Earth Data Analysis Center 3 Utah Department of Natural Resources 4 Natural color and color infrared imagery To prepare aerial photographs for analysis they were scanned at high resolution, collar information was cropped away, and adjacent images taken on the same day were merged using the Adobe Photoshop photomerge utility (Adobe Systems 2010). The utility uses a proprietary feature- matching algorithm to reference images to each other. The proprietary photomerge algorithm is similar to the Scale-Invariant Feature Transform algorithm described in Lowe (2004). Merged images were then georeferenced to current aerial photographs using manually placed control points for features visible in both the historical and current reference imagery. High resolution imagery from the National Agricultural Imagery Program (NAIP) flown in 2011 was used to reference the historical imagery. The NAIP imagery has a horizontal accuracy of 6m with a 95% confidence interval (USDA 2013). Any error present in the reference imagery is minimized during analysis because all photographs are aligned to the same imagery. I used a spline transformation in ArcMap to align the historical imagery with the current imagery (Esri 2014). For images with low distortion, as few as one control point per valley mile was used; in 4 other images, as many as ten control points per valley mile were used. Many control points were used for each image with a much higher density of points near the river. Spline transformation distributes distortion between points and is most accurate at each control point (Esri 2015). The spline transformation was particularly suited for use with the historical photographs where image distortion is common. Areas of high relief or where photographs were taken obliquely or off of vertical are difficult to georeference using traditional polynomial transformation, but the spline transformation correctly aligns the images to ground control points. Figure 2. Sample imagery from each time period used in analysis. Valley miles 203 and 204 are shown. Imagery Interpretation and Digitization Georeferenced images where interpreted and digitized at a scale of 1:5,000. The left bank, right bank, and islands were digitized with islands stored as polygon features and banks stored as line features. Left and right banks were digitized to correspond to the edges of the active channel as defined by the U.S. Army Corps of Engineers (Lichvar & McColley 2008). While the active channel is relatively stable compared to wetted area, streamflow at the time of photo and quality of the aerial imagery affected how the channel was interpreted and digitized. Each photo series was digitized by an analyst and then reviewed by a second analyst. In areas where there was disagreement, a consensus was reached. Some of the aerial photographs were low quality which made interpretation difficult. Specular reflection and tone differences between images, as well as damage and scanning errors (scratches, dirt, and fold lines), were common but could usually be avoided by using other photographs of the same area. Islands were digitized but sand or gravel bars were not. The active channel definition used during imagery interpretation characterizes exposed sand and cobble bars as within the active channel and more permanent vegetated islands that have upland characteristics as outside the active channel (Lichvar & McColley 2008). Vegetation, weathering such as erosion caused by precipitation, and other visible disturbances were indicators used to differentiate islands and bars. Many smaller exposed bars were present only in Reach 3 in the 1930s imagery and appeared to be part of a heavily braided system; these bars were not digitized as islands. 5 There were several instances of relief displacement in the 1930s imagery where cliffs and canyon walls obscure parts of the river. In the 1930s imagery, high relief areas with oblique imagery obscured 1.85 miles of channel edge or 0.66% of both banks. In these areas, the approximate location of the bank was estimated and digitized. Additional verification of bank location was not completed due to the short portions of the bank that were obscured. The edge of the active channel was not obscured in the photographs from the 1970s and 2011. Valley Delineation and Linear Referencing River miles have historically been used for linear referencing in studies on the San Juan River. The most common delineation used by the SJRIP begins with river mile zero just downstream of the Clay Hills boat ramp. In 2003 the SJRIP Biology Committee decided to continue using these mile delineations even though the delineated miles were not one mile in length (Gottleib 2003). There have been significant changes in channel sinuosity that have created discrepancies between current river distance and the standard river miles. Legacy river miles are useful for defining any arbitrary point along the river but are poor for summarizing and comparing linear and areal measures (Novak & Walker 2013). To enable accurate comparisons between reaches through time and to allow percent change measures, valley miles were delineated which address the issues presented by using river miles. Other recent studies have used valley miles for summarizing results along the San Juan River (e.g. Mangano 2013). A crosswalk between legacy river miles and current valley miles is attached to this report (Appendix A). Valley miles were mapped continuously from the historical confluence with the Colorado River to the top of the study area. To define valley miles, the valley floor was delineated using methods developed by the Vermont Agency of Natural Resources (VANR 2009). Digital elevation data, slope data, and aerial imagery were used to identify the edges of the valley floor. Areas where tributary valleys intersect the main valley floor were only included in the mapped valley floor where they appeared to be influenced by flood flows on the main stem. The valley centerline was derived from the valley floor as the line equidistant from both valley walls (VANR 2009). The valley centerline was then used to create a linear reference framework for the study area. Lines were projected perpendicular to the valley centerline at one mile increments and were used to create valley floor polygons referenced to the valley mile. Planform Channel Analysis To compare changes in the active channel, linear bank features and polygonal island features were converted into a polygon representation of the active channel. Valley miles and reach delineations were used to summarize the area of the active channel and islands. Island count, island area, active channel area and valley floor area were summarized for each valley mile and reach, and the whole study area. The complete coverage for the study area allowed calculation of areal measures that relate to fish habitat. Traditional measures of complexity calculated from river centerline and sample transects 6 such as braiding intensity and sinuosity are possible with the imagery, but were not completed during this analysis in favor of 100% sampling of channel and island features. Vegetation Sampling Bank armoring by vegetation is considered to be a significant control on erosion and bank avulsion (Cadol et al. 2010). To quantify the vegetation located on the banks of the river, continuous transects were placed parallel to the digitized active channel on both banks and islands. Sample points were placed at 37.4 meter intervals along the transects throughout the entire study reach (Figure 3). This point density was selected to provide a minimum of 150 sample points per valley mile. Presence or absence of vegetation was recorded for each sample point at a scale of 1:5,000. This protocol was repeated for each series of aerial photos. The sample point technique overestimates true canopy cover due to gaps in vegetation smaller than the photo resolution and shadows cast by the canopy (Frescino & Moisen 2012). The cover estimates derived from my analysis should not be considered analogous to true canopy cover but do serve as an index to changes in vegetation cover over time. Figure 3. An example reach including valley mile 203 with vegetation sampling points color coded by vegetation presence. Trees and large shrubs were easily identifiable from the aerial photos. Stands of cottonwood, Russian olive and tamarisk were classified as riparian vegetation. Willows were difficult to identify in the older imagery and may have been incorrectly classified as non-vegetation. In 1994 willow accounted for 6% of all riparian vegetation cover (Bliesner & Lamarra 2000) suggesting that any error in willow classification would have a relatively small effect on estimates of riparian vegetation cover in this analysis unless there was much more willow in the 1930s than in the 1990s. Results Throughout the study area the planform channel has changed dramatically since the 1930s; there has been a large decrease in channel area, island count, and total island area in the study area (Figure 4). Channel area decreased by 77.3% from 14,277 acres in the 1930s to 3,235 acres in 2011. A 34.7% decrease in island count between the 1930s and 1970s did not dramatically decrease total 7 island area due to a corresponding increase in average island size. Island area stayed relatively constant between the 1930s and 1970s but declined from 4,494 acres in the 1970s to 2,725 acres in 2011, a 39.4% decrease. Figure 4. Planform channel characteristics for the study area in the 1930s, 1970s, and 2010s. The geomorphic reaches identified by Bliesner and Lamarra (2000) exhibit distinct differences in channel characteristics and change (Figure 5). During the study period Reach 3 experienced the largest decrease in channel area (83%). The other study reaches demonstrate the same trend but with less severity (decreases between 58% and 76%). Between the 1930s and 2011, island count increased by over 200% in Reach 3, stayed relatively constant in Reach 4, and decreased dramatically in Reach 5 and 6 (by 67% and 61% respectively). Between the 1930s and 2010s, average island area increased in Reach 5 but decreased in all other reaches (Figure 5c). The largest decrease in average island area occurred in Reach 3 (12.1 acres to 2.8 acres, 76%). Despite the large decrease in average island size, total island area decreased only slightly in Reach 3 because of the increase in island count over the study period. Figure 5. Planform channel characteristics summarized for each reach. 8 Reach-level summaries conceal the high variation along the river. No two valley miles are alike. Despite similarities between each reach’s river miles, there are many river miles that do not follow the patterns seen in other parts of the reach (Figure 6). Channel area, when expressed as a percent of valley floor illustrates how the channel area has changed over time relative to maximum possible river area (Figure 6d). The percent of valley floor occupied by the river channel is a measure of channel confinement by the valley walls. Channel confinement by valley walls was high in the 1930s throughout reaches 3 and 4 but confinement by valley walls has since decreased as channel area decreased. Valley miles 135, 136, 233, and 234 have a consistently low level of confinement but this is likely the result of incorrect delineation of the valley walls. The decreasing amount of confinement by valley walls suggests that the river is being confined by other mechanisms. For the reach and valley-mile summaries, islands are split at divisions between analysis units. This is apparent in island count tabulations because islands that span a unit division will be double- counted if the island counts for adjacent units are added together. When analyzing the datasets that accompany this report, island count for individual valley miles or reaches should not be aggregated or added because islands that extend between multiple valley miles or reaches will be counted twice. 9 Figure 6. Planform channel characteristics summarized for each valley mile. 10