Effects of dams and geomorphic context on riparian forests of the Elwha River, Washington

Abstract

Understanding how dams affect the shifting habitat mosaic of river bottomlands is key for protecting the many ecological functions and related goods and services that riparian forests provide and for informing approaches to riparian ecosystem restoration. We examined the downstream effects of two large dams on patterns of forest composition, structure, and dynamics within different geomorphic contexts and compared them to upstream reference conditions along the Elwha River, Washington, USA. Patterns of riparian vegetation in river segments downstream of the dams were driven largely by channel and bottomland geomorphic responses to a dramatically reduced sediment supply. The river segment upstream of both dams was the most geomorphically dynamic, whereas the segment between the dams was the least dynamic due to substantial channel armoring, and the segment downstream of both dams was intermediate due to some local sediment supply. These geomorphic differences were linked to altered characteristics of the shifting habitat mosaic, including older forest age structure and fewer young Populus balsamifera subsp. trichocarpa stands in the relatively static segment between the dams compared to more extensive early‐successional forests (dominated by Alnus rubra and Salix spp.) and pioneer seedling recruitment upstream of the dams. Species composition of later‐successional forest communities varied among river segments as well, with greater Pseudotsuga menziesii and Tsuga heterophylla abundance upstream of both dams, Acer spp. abundance between the dams, and P. balsamifera subsp. trichocarpa and Thuja plicata abundance below both dams. Riparian forest responses to the recent removal of the two dams on the Elwha River will depend largely on channel and geomorphic adjustments to the release, transport, and deposition of the large volume of sediment formerly stored in the reservoirs, together with changes in large wood dynamics.

Introduction

Dam emplacement and operation can trigger changes to flow regime (Nilsson et al. , Poff and Zimmerman ), sediment flux (Vorosmarty et al. ), water‐transported plant propagules (Nilsson et al. ), large woody debris (Angradi et al. , Collins et al. ), and associated aquatic and floodplain ecosystems (Bunn and Arthington , Nilsson and Svedmark ). Floodplain vegetation responses to dams can vary among regions and river segments (Friedman et al. ), but some general outcomes have been documented. Flow regulation and decreased sediment loads can contribute to channel narrowing and simplification (Peipoch et al. ) and sometimes a transient increase in riparian forests that colonize former channels (Friedman et al. , Lobera et al. ). Channel simplification and reduced channel migration or widening rates can contribute to reduced vegetation recruitment, overall older age structures, and increases in late‐successional forests or more xeric vegetation (Merritt and Cooper , Kloehn et al. , Johnson et al. , Magdalena et al. ). Reduced low flows can lead to mortality of mesic tree seedlings and mature trees in semi‐arid regions, favoring non‐native species in some cases (Rood and Mahoney , Stromberg et al. , Vale et al. ), whereas increased low flows can contribute to vegetation expansion (Shafroth et al. ).

An understanding of the downstream effects of dams on riparian vegetation can be used to test and potentially extend the application of general concepts in river science. For example, riparian landscapes have been described as a shifting habitat mosaic, a unifying concept that articulates spatial and temporal habitat dynamics as a function of key biophysical drivers and feedbacks (Stanford et al. ). The mosaic of riparian habitats, composed of various fluvial geomorphic surfaces and associated vegetation, is typically altered and simplified downstream of dams (Friedman et al. , Petts and Gurnell , Graf , Egger et al. ). These effects of dams provide insights into the fluvial and ecological processes that sustain habitat mosaics in riparian ecosystems. Understanding how dams influence the shifting habitat mosaic is important to protect the many ecological functions and related goods and services that riparian forests provide (Naiman et al. ) and can inform efforts to restore riparian ecosystems (Webb and Erskine , Richardson et al. , Gonzalez et al. ).

In the Pacific Coastal Ecoregion of western North America (PCE), riparian forest development and succession have been described as dynamic patch mosaics (analogous to the shifting habitat mosaic) that are largely dependent on flow‐mediated sediment and wood dynamics (Latterell et al. , van Pelt et al. , Collins et al. ), soil formation, and nutrient cycling (Naiman et al. ), within the context of valley form and channel pattern (Beechie et al. ). At the site or river reach scale, these processes and feedbacks often result in simultaneous changes to fluvial geomorphology and riparian forests (van Pelt et al. ). Although the literature from the PCE has elucidated key processes associated with forest dynamics on unregulated rivers, few studies have examined the effects of dams on riparian forests in this region (Braatne et al. , Kloehn et al. ), though geomorphic changes associated with dams are better documented (e.g., Draut et al. , Gendaszek et al. ).

We studied the riparian forests of the Elwha River in Washington, USA, at a time when two large dams (32 and 64 m high) had been in place for 90 and 76 yr, respectively. Effects of the dams on fluvial processes and geomorphology along the river have been described in several earlier studies, revealing slower rates of channel change along the dam‐affected downstream river segments than upstream of the dams, though the channel still exhibited considerable mobility near the river mouth (Draut et al. ). The river segment between the dams changed considerably as a result of dam construction and operation, including a somewhat incised and armored channel and relatively static channel configuration (Pohl , Kloehn et al. ), representing a possible shift from an island‐braided to straight channel pattern in some reaches (Beechie et al. ).

We sought to answer four primary questions and address associated hypotheses related to Elwha River riparian forests and the effects of dams on these forests:

1. How does riparian forest composition and structure change during ecological succession? We hypothesized that overstory forest species composition and structure (i.e., cover, basal area, stem density) along the Elwha River would differ predictably with successional stage, as indicated by stand age and by proximity to the channel and associated fluvial disturbance.
2. How do riparian forest patches vary among different geomorphic settings? We hypothesized that forest age, structure, and species composition would differ predictably with geomorphic surface development from bars to floodplains and terraces, reflecting a shifting habitat mosaic driven by fluvial disturbance and ecological succession.
3. How do different geomorphic settings vary among different river segments (dammed, undammed)? We hypothesized that geomorphic surfaces downstream of the Elwha River dams would be older than those upstream, due to reduced fluvial dynamics downstream of the dams.
4. How do riparian forest patches vary among different river segments (dammed, undammed) with respect to species composition, stand structure, and age? We hypothesized that riparian forest patches downstream of the dams would have an older stand age distribution, altered stand structure, and later‐successional species composition than those upstream, because the dams disrupted the fluvial dynamics that drive the shifting habitat mosaic in riparian ecosystems.

In addition to these questions aimed at furthering our understanding of the downstream effects of dams on riparian forests, our results provide baseline data to which future studies of riparian forest dynamics on the Elwha River can be compared. This baseline data context is particularly novel on the Elwha River, where between September 2011 and November 2014 both dams were removed for ecosystem restoration (Warrick et al. ).

Methods

Study area

The Elwha River is located on the northern side of the Olympic Mountains in northwestern Washington State and drains approximately 833 km² of mountainous, forested terrain primarily in Olympic National Park along its 72 km length before reaching the Strait of Juan de Fuca (Fig. ). From its headwaters to its mouth, the Elwha River drops nearly 1370 m, and annual precipitation declines from ~600 to 100 cm along this elevation gradient (Duda et al. ). Average instantaneous discharge is 43 m³/s, and the 2‐, 10‐, and 100‐year recurrence interval floods are approximately 400, 752, and 1240 m³/s, respectively (Duda et al. , USGS streamflow gaging station 12045500, Elwha River at McDonald Bridge). Flood flows are typically associated with two different hydroclimatic processes. The largest annual peaks result from winter rains, primarily in November to January, and tend to be of relatively short duration. Smaller, longer‐duration peaks occur in association with snowmelt in late spring and early summer. Our study sites were located in alluvial valleys of the lower 29 km of the Elwha River, characterized by an island‐braided channel pattern and gravel‐bed channels (Beechie et al. , East et al. ).

Map of the Elwha River study area. Transects were located 31.7, 30.6, 30.4, 29.0, 28.3 (upper segment), 21.2, 19.5, 18.5, 17.1, 15.8 (middle segment), 6.0, 4.5, 3.4, 2.6, and 1.8 (lower segment) river km upstream of the river mouth.

Two dams constructed on the Elwha River to generate power for timber and paper mill operations had relatively minor effects on flow regime, but substantially altered sediment supply, river connectivity, and channel processes downstream. The 32 m high Elwha Dam, completed in 1913, was 7.9 km from the mouth of the Elwha River and created a 4 km long reservoir (Lake Aldwell). Glines Canyon Dam, a 64 m high structure completed in 1927, was located 21.6 km from the mouth and created a 4.5 km long reservoir (Lake Mills; USDOI BOR ). Power generation operations at Glines Canyon Dam resulted in significant flow fluctuations from 1927 to 1975. From 1975 to 2000, flow was largely “run‐of‐river” (i.e., inflows and outflows were roughly equal). After 2000, flow was still largely run‐of‐river, except during September and October when low flows were augmented to benefit native fish. The operations of Elwha Dam were always predominantly run‐of‐river. Peak flows were never significantly affected by either dam (Duda et al. ). A more dramatic effect of the two dams was the restriction of downstream sediment and wood transport and upstream anadromous fish movement. An estimated 21 ± 4 million m³ of sediment was trapped behind the two reservoirs prior to dam removal (Randle et al. ). Changes to sediment flux and flow regime contributed to significant armoring and stabilization of the active channel downstream of both dams, with the exception of the downstream‐most 4 km of the river, which exhibited some channel mobility (Pohl , Draut et al. ).

We measured riparian forests and described bottomland geomorphology in three segments of the Elwha River: an unregulated segment above both dams (hereafter the “upper” segment), a regulated segment between the two dams (“middle” segment), and a regulated segment below both dams (“lower” segment; Fig. ). The upper segment was located in an unconstrained valley between two bedrock canyons, ~28–32 river km above the Elwha River mouth, and had an island‐braided channel pattern. In the middle segment, our study sites were ~15–21 river km above the mouth, where the Elwha River flows through a largely unconstrained, alluvial valley, with either a single, wandering or anabranching gravel‐bed channel. The lower segment was characterized by an anabranching gravel‐bed channel and our study sites were within a ~5‐km stretch, beginning downstream of a narrow bedrock canyon ~1 km from the Elwha Dam and ending ~2 km upstream of the river mouth. Channel gradient decreased from upstream to downstream, approximately 0.015, 0.007, and 0.003 in the upper, middle, and lower segments, respectively (Kloehn et al. ).

Field methods

Within each river segment, we selected five cross‐valley transects to represent a range of common geomorphic surfaces and forest communities, based on field reconnaissance and examination of aerial photography. We distributed transects throughout the upstream to downstream extent of the three segments, but spacing between transects was not systematic and ranged from 0.4 to 1.7 km depending on factors such as private property or other access constraints. Transects extended between valley walls or high terraces that were >8 m above the channel except for the transects at river km 1.8 and 2.6, which extended to a levee that functionally bounded the east side of the bottomland. Transect cross‐valley widths were 319, 322, 390, 412, and 437 m in the upper segment, 319, 401, 534, 550, and 557 m in the middle segment, and 430, 476, 819, 1027, and 1174 m in the lower segment.

To characterize the overall geomorphic structure of the transects and river segments, we first divided the transects into distinct patch‐types based on dominant overstory vegetation and geomorphic position and measured the length of each patch along the transect line. We then used these patch‐type transect lengths, plot geomorphic surface‐age classes (see below), and topographic data to estimate the proportions of the transects occupied by general channel and landform types, including primary channels, perennial floodplain channels, non‐perennial floodplain channels, bars, floodplains, and fluvial terraces. We categorized the patches that contained plots based on geomorphic surface‐age classes and categorized the patches that did not contain plots based on topography relative to neighboring patches.

To measure plant community attributes in the forested areas on each transect, we placed 133 plots at random along the transect lines within the forested patches. Plots ranged in size from 20 m² (2 × 10 m in homogeneous patches with high tree densities) to 250 m² (10 × 25 m in patches with mature forest). Plots were rectangular and were perpendicular to the transect line. We added more than one plot per patch, depending on patch width (distance along transect): one or two plots in patches that extended 5–50 m along the transect line, two plots in 51‐ to 100‐m patches, three plots in 101‐ to 150‐m patches, and four plots in >150‐m patches. Because we were most interested in areas directly influenced by the river and hence more likely to be influenced by dam removal, we did not place plots in patches on surfaces with elevations >8 m above the primary channel water surface, and limited the number of plots on relatively high surfaces to one or two on the portion of the surface nearest the primary channel.

For each plot, we determined the elevation above and distance along the transect to the water surface of the nearest primary channel (disturbance channel; hereafter, elevation_dist and distance_dist), and, where applicable, the nearest secondary channel (groundwater channel; hereafter, elevation_gw and distance_gw). For our study, a primary channel was one that carried significant flow during high‐flow events and thus was likely to contribute to fluvial disturbance. Secondary channels were small, perennial, often groundwater‐fed, and not likely to be associated with fluvial disturbance of nearby plots, but likely to indicate the approximate elevation of alluvial groundwater. We surveyed locations and elevations of the plots using a Leica TCR705 total station during the 2003 growing season for the transects at river km 18.5, 19.5, and 21.2 and during the 2004 growing season for the remaining transects. We tied the relative elevation of surveyed points along the transects to a final water surface survey during a low‐flow period (11.3 m³/s) in August 2005.

We estimated overstory tree cover in the plots, by species, and converted cover values to the midpoints of Daubenmire cover classes for analysis: 0.1–5% = 2.5%, 5–25% = 15.5%, 25–50% = 37.5%, 50–75% = 62.5%, 75–95% = 85%, and 95–100% = 97.5% (Daubenmire ). Further, we measured diameter at breast height (dbh) for all individuals with dbh > 5 cm and used those data to calculate basal area and stem density per hectare for each species. We calculated overstory species importance as the sum of relativized cover, stem density, and basal area in each plot. In addition, we estimated cover of tree species in the understory (saplings and seedlings, hereafter “understory trees”) in two 2 × 4 m subplots in each plot; subplots were averaged for analyses. Salix scouleriana and Salix sitchensis were sometimes difficult to distinguish, so we lumped them for all analyses. We collected most plant community composition data during the 2004 growing season, except for the transects at river km 18.5, 19.5, and 21.2 which we sampled during the 2003 growing season.

We estimated stand age for every plot by collecting tree cores from the largest individuals of characteristic species within or adjacent to the plots. We measured stem diameter at the core height and at breast height for each cored tree. We collected cores from a total of 247 trees near or within 98 plots, mounted and sanded the cores, and counted the annual rings under a microscope. For cores that did not include the pith (i.e., off center), we estimated additional years that were missing. We also added years to the annual ring counts to account for time required for an individual of a given species to grow from the ground surface to the coring height, using published values (King , Harrington and Curtis , Fried et al. , DeBell , Foiles et al. ). These age estimates were likely accurate to ±3 yr, similar to reports in other studies of western riparian vegetation (Samuelson and Rood , Braatne et al. ). For plots from which cores were not collected, we estimated age based on diameter–age linear regressions that we developed for five of the tree species (Appendix S1: Fig. S1). R² values were >0.49 for all models, and P‐values <0.012. The oldest tree cored in a plot or estimated by regression was used to estimate the establishment year of woody vegetation in the plot.

Next, we assigned the vegetation plots to one of six geomorphic surface‐age classes: (1) open pioneer bar (1–10 years old and <50% woody vegetative cover; N = 2); (2) woody pioneer bar (1–10 years old and >50% cover; N = 9); (3) developing floodplain (11–30 years old; N = 29); (4) established floodplain (31–50 years old; N = 19); (5) transitional fluvial terrace (51–90 years old; N = 39); and (6) mature fluvial terrace (91–200 years old; N = 35). This classification largely followed the model developed for the Queets River (Latterell et al. , van Pelt et al. ) and assumed that other aspects of these surface‐age classes were similarly correlated with age in our study area (but see Results).

To evaluate tree seedling recruitment on younger surfaces, we measured seedling densities on bars with no overstory cover adjacent to the main channel (hereafter “recruitment areas”) during the 2004 growing season. We counted and identified tree seedlings (plants <1 m tall) to species in 0.1‐ or 1.0‐m² plots (smaller plots along transects with higher seedling densities) placed at 1‐m intervals from the water's edge to the edge of established vegetation on each transect. For each transect, we calculated the mean stem density, the length of the recruitment area along the transect line, and a total seedling stem count in a 1 m wide strip through the recruitment area (the product of the mean stem density and the length of the recruitment zone).

Statistical methods

We used non‐metric multidimensional scaling (NMDS; Kruskal , Mather ) with Sørenson (Bray–Curtis) distances in PC‐ORD version 5.0 (McCune and Mefford ) to ordinate log‐transformed overstory species importance. We excluded species present in <2% of plots, and plots with no overstory, from the analysis. We used random starting configurations with 250 runs with real data, examined NMDS scree plots to select between one and six dimensions for each ordination, and ran 500 iterations using the best starting configuration for the final solution. We used Varimax rotation to maximize correspondence between the ordination axes and species importance.

We tested for correlations between the Varimax‐rotated ordination axes and overstory species importance, as well as correlations between the ordination axes and other plot characteristics using Kendall's tau (τ) analyses in PC‐ORD. To avoid discussion of weak correlations, we reported only correlations with |τ| > 0.3. Further, we examined causal relationships among abiotic plot characteristics (stand age, elevation_dist, elevation_gw, distance_dist, distance_gw), as well as between those plot characteristics and several vegetation measures (total overstory cover, stem density, basal area, and understory tree cover), using single‐factor linear regression analysis (REG procedure, SAS 9.3; SAS Institute, Cary, North Carolina, USA). We used corrected Akaike's information criterion scores to choose between linear and quadratic models. To correct non‐normality and unequal variances, we log‐transformed total stem density, basal area, and understory tree cover for regression analyses.

We separated the plots into different overstory community types using hierarchical cluster analysis of log‐transformed overstory species importance, with Sørenson (Bray–Curtis) distances and flexible β = −0.25 as the linkage method in PC‐ORD. We excluded species present in <2% of plots, and plots with no overstory, from the analysis and placed plots with no overstory into an a priori “no overstory” community type. We stopped merging clusters when all resulting clusters represented clear and ecologically interesting differences in overstory composition. We identified misclassified plots in the resulting clusters using non‐parametric discriminant analysis (DISCRIM procedure, SAS 9.3; SAS Institute), using the cross‐validation feature and the value for k nearest neighbors (k = 8) that corrected the most obvious misclassifications while minimizing the total number of misclassifications (15 plots).

We compared log‐transformed overstory species importance and log‐transformed, relativized understory tree species cover among overstory community types, geomorphic surface‐age classes, and river segments using non‐parametric, multi‐response permutation procedures (MRPP; Mielke and Berry ) in PC‐ORD, with Sørenson (Bray–Curtis) distances and weights = n/sum(n). We excluded species present in <2% of plots, and plots with no species present, from the analyses. For the segment comparisons, we analyzed transect means rather than plot data because transects were the scale of replication for segment comparisons and MRPP could not account for the nested design. We used indicator species analysis (ISA) to identify species indicative of each group (Dufrene and Legendre ) and tested the significance of the ISA indicator values (hereafter “ISA IV”) with a Monte Carlo method using 4999 randomizations in PC‐ORD. We named the community types based on species with ISA IV > 15, in order of mean species importance across the plots within each cluster.

We compared other biotic and abiotic plot characteristics among community types, geomorphic surface‐age classes, and river segments using mixed linear model analyses of variance (ANOVA; MIXED procedure, SAS 9.3; SAS Institute). For all comparisons, the models included river segment as a fixed effect and transect (nested within segment) and patch (nested within segment and transect) as random effects, with degrees of freedom calculated using the Kenward–Roger method (Kenward and Roger ). We included community types and geomorphic surface‐age classes as additional fixed effects in separate analyses. To avoid pseudoreplication, when two plots with the same geomorphic surface‐age class and community type, and with similar ages and elevations above the channel, occurred within the same small patch (<50 m), we deleted one of the plots for river segment comparisons (14 plots). When main effects were significant, we conducted pairwise comparisons using Tukey‐adjusted least squared means tests. To meet model assumptions, we cube‐root‐transformed elevation_gw, distance_dist, and distance_gw and log‐transformed total basal area, total stem density, and total understory tree cover for both community type and geomorphic surface‐age class comparisons, and cube‐root‐transformed stand age for geomorphic surface‐age class comparisons. Transformations could not correct heteroscedasticity among community types in NMDS axis 3 scores, stand age, elevation_dist, total overstory cover, basal area, and stem density, or among geomorphic surface‐age classes in NMDS axis 1 and 2 scores, stand age, elevation_dist, elevation_gw, and total overstory cover, so we added a group statement to those models to calculate separate variances for each type or class (Littell et al. ).

We compared the frequencies of the overstory community types among geomorphic surface‐age classes with Pearson chi‐square tests (FREQ procedure, SAS 9.3; SAS Institute). For post hoc comparisons of individual community types, we accounted for variation in sample size among geomorphic surface‐age classes by calculating expected frequencies as the total number of plots containing a particular community type, multiplied by the percentage of total plots containing each geomorphic surface‐age class. Because expected counts in some cells were low, we used exact statistics instead of asymptotic methods.

We compared the proportions of the transects composed of each general landform and channel type among the three river segments with one‐way ANOVAs (MIXED procedure, SAS 9.3; SAS Institute). To meet model assumptions, we arcsine‐square‐root‐transformed pioneer bar and floodplain proportions.

In addition, we compared the frequency of the geomorphic surface‐age classes and overstory community types among the three river segments by comparing the proportions of the plots on each transect that contained each geomorphic surface‐age class or community type with one‐way ANOVAs (MIXED procedure, SAS 9.3; SAS Institute). To meet model assumptions, we log‐transformed the proportions for analysis.

Finally, we compared the lengths of seedling recruitment areas along the transect lines, transect means of seedling stem densities, and transect means of total seedling stem counts (length × density) among river segments, with one‐way ANOVAs (MIXED procedure, SAS 9.3; SAS Institute). To meet model assumptions, we square‐root‐transformed transect lengths of seedling recruitment areas and log‐transformed seedling stem densities and total stem counts for analysis.

Results

Elwha River riparian forest species

Generally, the species composition of Elwha River riparian forests was similar to other systems in the Pacific Coastal Ecoregion (Franklin and Dyrness , Naiman et al. , Acker et al. ), with mixtures of coniferous taxa such as Pseudotsuga menziesii (Mirb.) Franco (Douglas fir), Abies grandis (Douglas ex D. Don) Lindl (grand fir), and Thuja plicata Donn ex D. Don (western redcedar), as well as deciduous taxa such as Acer macrophyllum Pursh (bigleaf maple), Acer circinatum Pursh (vine maple), Alnus rubra Bong. (red alder), Salix scouleriana Barratt ex Hook (Scouler willow), Salix sitchensis Sanson ex Bong. (Sitka willow), and Populus balsamifera L. subsp. trichocarpa (Torr. & A. Gray ex Hook) Brayshaw (black cottonwood; Appendix S1: Table S1). Two regionally dominant trees, Tsuga heterophylla (Raf) Sarg. (western hemlock) and Picea sitchensis (Bong.) Carrière (Sitka spruce), occurred infrequently at our study sites relative to other parts of the PCE. Numerous shrub species also grow in the Elwha riparian zone, including Holodiscus discolor (Pursh) Maxim. (ocean spray), Symphoricarpos albus (L.) S.F. Blake (snowberry), Rubus spectabilis Pursh (salmon berry), and other Rubus L. species. Throughout the text, we refer to overstory tree species by genus, with Acer_m and Acer_c for the two Acer species. Nomenclature follows the USDA plants database (http://plants.usda.gov).

Riparian forest succession—relationships between forest structure, stand age, and channel proximity

Not surprisingly, elevation_dist and elevation_gw were strongly correlated across plots (τ = 0.81), as were distance_dist and distance_gw (τ = 0.73). In regression analyses, elevation above the channel increased with stand age (Appendix S1: Fig. S2). Total basal area also increased with stand age, whereas total stem density decreased, and total overstory cover was highest at intermediate stand ages. Further, total basal area increased with elevation above the channel, whereas total stem density was lowest at intermediate elevations. Total understory tree cover was highest at intermediate distances from the channel (Appendix S1: Fig. S3). There were no significant trends in stand age, elevation above the channel, total overstory cover, total basal area, or total stem density with distance from the channel, and no significant trends in total understory tree cover with stand age, elevation above the channel, total overstory cover, total basal area, or total stem density (P > 0.05; data not shown).

Riparian forest succession—NMDS ordination of overstory species composition

Non‐metric multidimensional scaling ordination of overstory species importance (Fig. ) yielded three axes, describing 38% (axis 1), 34% (axis 2), and 19% (axis 3) of the variance (total R² = 92%). Final stability was 0.001 and final stress was 12.0. Axis 1 was positively correlated with Pseudotsuga and Abies and negatively correlated with Alnus (Table ). Axis 2 was positively correlated with Acer_m, Acer_c, and Thuja and negatively correlated with Salix and Alnus. Axis 3 was positively correlated with Populus, Thuja, Acer_m, and Abies and negatively correlated with Alnus.

Varimax‐rotated scores for axis 1 and axis 2 from non‐metric multidimensional scaling ordination of log‐transformed overstory species importance. Overlays of 10 overstory community types and five geomorphic surface‐age classes are also shown. Joint plots indicate correlations between the ordination axes, species importance, and plot characteristics with Kendall's |τ| > 0.3. ED, elevationdist; EW, elevationgw; Au, understory Alnus cover; Su, understory Salix cover.

Correlations between forest overstory NMDS axis scores and plot measurements, including overstory and understory forest characteristics and proximity to the channel (Kendall's τ)

Plot measurements	Axis 1	Axis 2	Axis 3
Overstory Salix importance	−0.03	−0.52	−0.29
Overstory Alnus importance	−0.54	−0.36	−0.58
Overstory Populus importance	−0.16	−0.08	0.51
Overstory Pseudotsuga importance	0.66	0.11	0.17
Overstory Abies importance	0.47	0.07	0.31
Overstory Tsuga importance	0.23	0.10	0.10
Overstory Thuja importance	0.11	0.40	0.37
Overstory Acerm importance	0.22	0.71	0.36
Overstory Acerc importance	0.18	0.47	0.14
Overstory Picea importance	−0.12	0.09	−0.08
Overstory stand age	0.31	0.56	0.35
Overstory total cover	0.03	0.20	0.25
Overstory total basal area	0.11	0.25	0.27
Overstory total stem density	−0.16	−0.28	−0.22
Understory Salix cover	−0.11	−0.32	−0.10
Understory Alnus cover	−0.13	−0.31	−0.24
Understory Populus cover	−0.06	−0.06	0.09
Understory Pseudotsuga cover	0.04	−0.08	−0.15
Understory Abies cover	0.16	−0.02	0.07
Understory Tsuga cover	0.11	−0.11	−0.06
Understory Thuja cover	0.04	0.09	0.09
Understory Acerm cover	0.10	0.17	0.13
Understory Acerc cover	0.12	0.27	0.20
Understory total cover	0.08	0.12	0.06
Channel proximity: elevationdist	0.33	0.28	0.16
Channel proximity: elevationgw	0.34	0.39	0.25
Channel proximity: distancedist	0.07	0.14	0.05
Channel proximity: distancegw	0.04	0.08	0.04

Notes

Strong correlations (|τ| > 0.3) are shown in bold. NMDS, non‐metric multidimensional scaling.

Stand age was positively correlated with all three axes (Table ). Elevation_dist and elevation_gw also were positively correlated with axis 1, and elevation_gw was positively correlated with axis 2. Understory Alnus and Salix cover were negatively correlated with axis 2, but no other understory tree species were strongly correlated with any axis. Distance from the channel, total overstory cover, basal area, stem density, and total understory tree cover were not strongly correlated with any axis.

Riparian forest succession—overstory community types

Cluster analysis of overstory species importance divided the plots into ten overstory community types (in addition to the a priori “no overstory” community type), which retained ~55% of the information in the overstory importance data. All species except Picea were significant indicators of at least one community type (Appendix S1: Table S1). Based on species with ISA IV > 15, the community types can be described as (1) “no overstory” (N = 10), (2) Alnus–Salix (N = 22), (3) Alnus (N = 18), (4) Alnus–Populus (N = 24), (5) Pseudotsuga–Alnus (N = 13), (6) Acer_m–Populus–Thuja (N = 10), (7) Pseudotsuga–Abies (N = 11), (8) Pseudotsuga–Tsuga (N = 5), (9) Pseudotsuga–Acer_c (N = 6), (10) Thuja–Acer_m–Acer_c (N = 5), and (11) Acer_m–Acer_c (N = 9). Overstory species‐specific importance differed significantly between all community types, except for between Pseudotsuga–Tsuga and both Pseudotsuga–Acer_c and Thuja–Acer_m–Acer_c communities (Fig. a; MRPP, A = 0.55, P < 0.0001; pairwise comparisons, Bonferroni‐adjusted α = 0.001).

Mean relativized overstory tree species importance and understory tree species cover for (a, b) overstory community types, (c, d) geomorphic surface‐age classes, and (e, f) river segments. Overstory species importance was calculated as the sum of relativized overstory cover, basal area, and stem density. For overstory community types, none, “no overstory”; Al‐Sa, Alnus–Salix; Al, Alnus; Al‐Po, Alnus–Populus; Ps‐Al, Pseudotsuga–Alnus; Ac‐Po, Acerm–Populus–Thuja; Ps‐Ab, Pseudotsuga–Abies; Ps‐Ts, Pseudotsuga–Tsuga; Ps‐Ac, Pseudotsuga–Acerc; Th‐Ac, Thuja–Acerm–Acerc; Ac, Acerm–Acerc. For geomorphic surface‐age classes, OB, open pioneer bars; WB, woody pioneer bars; DF, developing floodplains; EF, established floodplains; TT, transitional fluvial terraces; MT, mature fluvial terraces. Plots with no overstory trees were excluded from the overstory means and multi‐response permutation procedures (MRPP) analyses, and plots with no tree saplings or seedlings in the understory were excluded from the understory means and MRPP analyses.

Scores for all three NMDS axes differed among community types, with higher axis 1 scores for Pseudotsuga–Tsuga, Pseudotsuga–Acer_c, Pseudotsuga–Abies, and Pseudotsuga–Alnus communities, higher axis 2 scores for Thuja–Acer_m–Acer_c, Pseudotsuga–Acer_c, and Acer_m–Acer_c communities, and higher axis 3 scores for Acer_m–Populus–Thuja communities (Figs. , a‐c).

Mean non‐metric multidimensional scaling (NMDS) ordination scores and other plot characteristics that differed significantly among the 11 overstory community types: (a) NMDS axis 1 (F9,111 = 76.4, P < 0.0001), (b) NMDS axis 2 (F9,99 = 29.9, P < 0.0001), (c) NMDS axis 3 (F9,96 = 20.9, P < 0.0001), (d) stand age (F10,19 = 20.1, P < 0.0001), (e) elevationdist (F10,16 = 10.2, P < 0.0001), (f) elevationgw (F10,65 = 11.4, P < 0.0001), (g) total overstory cover (F9,23 = 4.6, P = 0.001), (h) total basal area (F9,23 = 4.4, P = 0.002), and (i) total stem density (F9,23 = 4.8, P = 0.001). Error bars are 1 SEM. Within each graph, different letters indicate significant differences among community types (Tukey‐adjusted least squared means, P < 0.05).

The Pseudotsuga–Tsuga, Pseudotsuga–Acer_c, Acer_m–Acer_c, Acer_m–Populus–Thuja, Thuja–Acer_m–Acer_c, and Pseudotsuga–Abies communities occurred in the oldest stands, followed by the Pseudotsuga–Alnus, Alnus–Populus, and Alnus communities (Fig. d). The Alnus–Salix and “no overstory” communities were youngest. The Pseudotsuga–Tsuga communities occurred at the highest elevations above the channel, followed by the Pseudotsuga–Acer_c, Thuja–Acer_m–Acer_c, Pseudotsuga–Abies, and Acer_m–Populus–Thuja communities (Fig. e, f). The Pseudotsuga–Alnus, Alnus–Populus, Alnus, and Alnus–Salix communities occurred at the lowest elevations. Total overstory cover was lower in Alnus–Salix, Alnus, and Pseudotsuga–Tsuga communities than in Acer_m–Populus–Thuja and Acer_m–Acer_c communities (Fig. g). Alnus–Salix communities also had lower total basal area and higher total stem density than many other communities (Fig. h, i). Plot distances from the channel and total understory tree cover did not differ significantly among community types (data not shown).

Understory tree species composition differed among overstory community types, with significantly different species composition in “no overstory” communities than in Alnus, Pseudotsuga–Alnus, Acer_m–Populus–Thuja, Pseudotsuga–Abies, Pseudotsuga–Acer_c, and Acer_m–Acer_c communities, as well as significantly different species composition in Alnus–Salix communities than in Pseudotsuga–Acer_c and Acer_m–Acer_c communities (Fig. b; MRPP, A = 0.12, P < 0.0001; pairwise comparisons, Bonferroni‐adjusted α = 0.0009). Understory Salix and Populus were significant indicators of the “no overstory” community type (ISA, IV = 23, P = 0.02; and IV = 38, P = 0.0008, respectively), understory Alnus was a significant indicator of the Alnus–Salix community type (ISA, IV = 19, P = 0.04), understory Pseudotsuga was a significant indicator of the Pseudotsuga–Tsuga community type (ISA, IV = 20, P = 0.04), and understory Acer_c was a significant indicator of the Pseudotsuga–Acer_c community type (ISA, IV = 23, P = 0.02).

Geomorphic surface‐age classes and forest composition and structure

Mature fluvial terraces occurred at higher elevations above the channel than woody pioneer bars and developing floodplains for elevation_dist, and higher than woody pioneer bars, floodplains, and transitional fluvial terraces for elevation_gw (Fig. d, e). Total basal area was higher on mature fluvial terraces than on floodplains, higher on transitional fluvial terraces than on developing floodplains, higher on established floodplains than on pioneer bars, and higher on developing floodplains than on open pioneer bars (Fig. g). Total overstory cover was higher on mature fluvial terraces than on woody pioneer bars and developing floodplains (Fig. f). In contrast, total stem density was higher on developing floodplains than on mature fluvial terraces and open pioneer bars (Fig. h). Plot distances from the channel and total understory tree cover did not differ significantly among geomorphic surface‐age classes (data not shown). The high variability and small sample size of open pioneer bars (n = 2) limited the utility of statistical comparisons between this and the other geomorphic surface‐age types.

Mean plot characteristics and non‐metric multidimensional scaling (NMDS) ordination scores that differed significantly among the six geomorphic surface‐age classes: OB, open pioneer bars; WB, woody pioneer bars; DF, developing floodplains; EF, established floodplains; TT, transitional fluvial terraces; MT, mature fluvial terraces. (a) NMDS axis 1 (F4,54 = 13.6, P < 0.0001), (b) NMDS axis 2 (F4,61 = 40.8, P < 0.0001), (c) NMDS axis 3 (F4, 104 = 8.7, P < 0.0001), (d) elevationdist (F5,66 = 3.0, P = 0.02), (e) elevationgw (F5,48 = 4.3, P = 0.003), (f) total overstory cover (F5,55 = 4.2, P = 0.003), (g) total basal area (F5,90 = 11.6, P < 0.0001), and (h) total stem density (F5,95 = 4.3, P = 0.002). Error bars are 1 SEM. Within each graph, different letters indicate significant differences among geomorphic surface‐age classes (Tukey‐adjusted least squared means, P < 0.05).

Overstory species‐specific importance differed significantly between most geomorphic surface‐age classes, but not between woody pioneer bars and developing floodplains, nor between established floodplains and transitional fluvial terraces (Fig. c; MRPP, A = 0.20, P < 0.0001; pairwise comparisons, Bonferroni‐adjusted α = 0.005). Salix was a significant indicator of woody pioneer bars (ISA, IV = 55, P = 0.002), whereas Acer_m and Acer_c were significant indicators of mature fluvial terraces (ISA, IV = 46, P = 0.0002; and IV = 44, P = 0.002, respectively). Scores for all three NMDS axes also differed among geomorphic surface‐age classes, indicating higher Pseudotsuga, Abies, Thuja, Acer_m, Acer_c, and Populus importance and lower Alnus and Salix importance on fluvial terraces than on developing floodplains and woody pioneer bars, and particularly high Acer_m and Acer_c, importance and low Salix importance on mature fluvial terraces (Figs. , a‐c).

The different geomorphic surface‐age classes also were associated with different overstory community types (χ² = 184.7, df = 50, P < 0.0001; Fig. ). “No overstory” communities occurred on the youngest geomorphic surface‐age classes, followed by Alnus–Salix, Alnus, and Alnus–Populus communities. The remaining community types occurred almost exclusively on fluvial terraces, resulting in increasing diversity of overstory community types with increasing geomorphic surface age (Fig. ). Pseudotsuga–Alnus communities occurred mainly on transitional fluvial terraces, Pseudotsuga–Abies and Acer_m–Populus–Thuja communities occurred on established floodplains and transitional and mature fluvial terraces, Pseudotsuga–Tsuga, Pseudotsuga–Acer_c, and Thuja–Acer_m–Acer_c communities occurred only on transitional and mature fluvial terraces, and Acer_m–Acer_c communities occurred only on mature fluvial terraces.

Proportion of plots containing each overstory community type on each geomorphic surface‐age class. OB, open pioneer bars; WB, woody pioneer bars; DF, developing floodplains; EF, established floodplains; TT, transitional fluvial terraces; MT, mature fluvial terraces. Numbers at the tops of the bars indicate the total sample size for each geomorphic surface‐age class. Asterisks mark frequencies with >0.5 relative deviations from expected frequencies, for community types with significant differences between expected and observed frequencies: “no overstory,” χ2 = 57.6, df = 5, P < 0.0001, from left to right, relative deviations = 12.3, 6.4, and 0.4; Alnus–Salix, χ2 = 38.5, df = 5, P < 0.0001, relative deviations = 1.7, 2.1, and −0.1; Pseudotsuga–Alnus, χ2 = 14.8, df = 5, P = 0.03, relative deviations = −0.7, 1.6, and −0.4; Acerm–Acerc, χ2 = 25.2, df = 5, P = 0.009, relative deviation = 2.8.

Understory tree species composition differed among geomorphic surface‐age classes, with significantly different species composition on woody pioneer bars and developing floodplains than on transitional and mature fluvial terraces (Fig. d; MRPP, A = 0.09, P < 0.0001; pairwise comparisons, Bonferroni‐adjusted α = 0.003). Understory Salix and Alnus were significant indicators of woody pioneer bars (ISA, IV = 29, P = 0.045; and IV = 28, P = 0.04, respectively).

River segment comparisons—geomorphic settings

Over the full lengths of the transects extending to an elevation of >8 m above the channel (or to the levee for two transects in the lower segment), bars comprised a greater proportion of the transects in the upper segment than in the middle and lower segments (39% ± 6% compared to 2% ± 1% and 6% ± 2%, respectively; Fig. a). Correspondingly, fluvial terraces comprised a lesser proportion of the transects in the upper segment than in the middle and lower segments (27% ± 10% compared to 77% ± 3% and 65% ± 7%, respectively). Although, on average, floodplains also comprised a greater proportion of the transects in the upper segment than in the middle segment, this difference was not significant. Proportions of the transects composed of channels did not differ significantly among river segments. Total transect distances were greater in the lower segment than in the upper segment (Table ). However, this difference did not explain the differences in proportional abundance of bars or fluvial terraces, as including total transect distance as a covariate did not alter the results of river segment comparisons for any landform or channel type.

(a) Mean proportions of total transect distances composed of different general landform and channel types for each river segment, (b) mean proportions of plots containing different overstory community types within each segment, and (c) mean proportions of plots containing different geomorphic surface‐age classes within each segment. Different letters indicate significant differences among segments (Tukey‐adjusted LS means, P < 0.05) in (a) pioneer bars (F2,12 = 29.9, P < 0.0001) and terraces (F2,12 = 9.2, P = 0.004), and (b) Alnus–Salix communities (F2,12 = 6.2, P = 0.01), Alnus–Populus communities (F2,12 = 6.9, P = 0.01), Acerm–Populus–Thuja communities (F2,12 = 4.7, P = 0.03), Pseudotsuga–Tsuga communities (F2,12 = 5.8, P = 0.02), Thuja–Acerm–Acerc communities (F2,12 = 6.9, P = 0.01), and Acerm–Acerc communities (F2,12 = 8.6, P = 0.005).

River segment means and ANOVA comparisons

Plot measurements	River segment (mean ± 1 SEM)	df	F	P
Lower	Middle	Upper
Stand age (yr)	56 ± 6ab	86 ± 5a	50 ± 7b	12	5.2	0.02
Elevationdist (m)	1.64 ± 0.12	1.87 ± 0.14	2.79 ± 0.29	11	4.1	0.047
Elevationgw (m)	1.61 ± 0.12	1.94 ± 0.14	2.42 ± 0.31	11	1.0	ns
Distancedist (m)	100.7 ± 13.0	70.6 ± 11.2	92.9 ± 12.7	13	0.4	ns
Distancegw (m)	86.3 ± 12.9	43.2 ± 5.2	77.9 ± 12.8	13	1.2	ns
NMDS axis 1	−0.21 ± 0.06b	0.01 ± 0.08ab	0.33 ± 0.16a	12	5.3	0.02
NMDS axis 2	0.02 ± 0.08	0.18 ± 0.09	−0.21 ± 0.10	13	2.3	ns
NMDS axis 3	0.17 ± 0.07a	0.03 ± 0.08ab	−0.33 ± 0.08b	14	6.5	0.01
Total overstory cover (%)	108.5 ± 6.4a	111.2 ± 6.4a	73.3 ± 6.5b	12	5.1	0.03
Total basal area (m2/ha)	63.7 ± 15.5	59.5 ± 6.5	56.4 ± 10.7	92	0.1	ns
Total stem density (/ha)	1186 ± 147	882 ± 160	1566 ± 361	87	2.1	ns
Total understory tree cover (%)	7.1 ± 1.9	5.9 ± 1.6	5.2 ± 1.1	87	0.3	ns
Total transect distance (m)	785 ± 147a	472 ± 48ab	376 ± 24b	12	6.1	0.01

Within a row, means marked with different superscript letters are significantly different (Tukey‐adjusted least squared means, P < 0.05). NMDS, non‐metric multidimensional scaling; ns, not significant.

⁴Denominator degrees of freedom. Numerator df = 2 for all comparisons.

River segment comparisons—forest composition and structure

Stands were younger in plots in the upper segment than in the middle segment (Table ). In the lower segment, mean stand ages on most transects were similar to stand ages in the upper segment; however, stand ages on the uppermost transect of the lower segment were similar to those in the middle segment (Appendix S1: Fig. S4). In contrast, elevation_dist was higher in the upper segment than in the lower segment, although the post hoc means comparison was not significant (Tukey, P = 0.06). Total overstory cover was lower in the upper segment than in the middle and lower segments. Elevation_gw, distance from the channel, total basal area, total stem density, and total understory tree cover did not differ significantly among river segments.

Both overstory and understory species importance differed between the upper segment and the middle and lower segments (Fig. e, f; MRPP, A = 0.18, P = 0.003 and A = 0.11, P = 0.006, respectively; pairwise comparisons, Bonferroni‐adjusted α = 0.016). Overstory Salix and Tsuga and understory Pseudotsuga were significant indicators of the upper segment (ISA, IV = 56, P = 0.03; IV = 57, P = 0.02; and IV = 61, P = 0.04, respectively), whereas overstory Populus and Thuja were significant indicators of the lower segment (ISA, IV = 46, P = 0.006; and IV = 64, P = 0.02, respectively). Non‐metric multidimensional scaling axis 1 scores were higher in the upper segment than in the lower segment, whereas NMDS axis 3 scores were lower in the upper segment than in the lower segment, indicating higher Pseudotsuga importance in the upper segment and higher Populus and Thuja importance in the lower segment (Table ).

Alnus–Salix communities occurred more frequently in the upper segment than in the middle segment, and Pseudotsuga–Tsuga communities occurred only in the upper segment (Fig. b). In contrast, Acer_m–Acer_c communities did not occur in the upper segment and occurred most frequently in the middle segment. Alnus–Populus and Acer_m–Populus–Thuja communities also did not occur in the upper segment, but occurred most frequently in the lower segment. Thuja–Acer_m–Acer_c communities occurred only in the lower segment. Relative frequencies of other overstory community types did not differ significantly among river segments (Fig. ).

River segment comparisons—tree seedlings on channel margins and bars

Transects at river km 5.6 (lower segment), 18.5, 19.5, and 21.2 (middle segment) included no seedling recruitment areas (i.e., unvegetated channel margins and gravel bars). Among transects that did support new seedlings, recruitment areas were more extensive in the upper segment than in the middle segment (Fig. a). Therefore, total seedling stem counts in the recruitment areas were higher in the upper segment than in the middle segment (Fig. c), even though mean seedling densities did not differ significantly among river segments (Fig. b). Species‐specific total seedling stem counts in the upper segment also differed significantly from the middle and lower segments (MRPP, A = 0.17, P = 0.01), although species‐specific stem densities did not. Pseudotsuga and Populus seedlings were significant indicators of recruitment areas in the upper segment compared to the middle and lower segments (ISA, IV = 87, P = 0.0004; and IV = 48, P = 0.03, respectively).

Tree seedling recruitment on bare channel margins on the three river segments: (a) mean lengths of seedling recruitment areas along the transect lines, defined as areas adjacent to the main channel, on shallow slopes, with no overstory vegetation (F2,12 = 12.6, P = 0.001); (b) mean seedling stem density; and (c) mean total seedling counts in a 1 m wide strip through each recruitment area (F2,12 = 8.7, P = 0.005), estimated by multiplying the mean seedling stem density by the length of each recruitment area. Error bars are one SEM. Different letters indicate significantly different means (Tukey‐adjusted least squared means, P < 0.05).

Discussion

Riparian forest succession across geomorphic settings

Changes in overstory vegetation with stand age and geomorphic position represent toposequences and reflect co‐development of stand composition, structure, and geomorphology through time, which has been well described within the PCE (Balian and Naiman , van Pelt et al. , Naiman et al. ) and globally (Corenblit et al. ). As we expected, on the Elwha River, overstory stem density, basal area, cover, and dominant species varied predictably with stand age and geomorphic position. Older stands tended to occur on surfaces at higher elevations above the channel (i.e., fluvial terraces) and tended to have higher total overstory cover, tree basal areas, and stem densities than younger stands on woody pioneer bars, and developing floodplains. The forest overstory in older stands on fluvial terraces tended to contain more Pseudotsuga, Abies, Tsuga, Thuja, and especially Acer_m and Acer_c, and less Salix and Alnus, than in younger stands on woody pioneer bars and developing floodplains (Figs. a, c, ). Tree recruitment in the understory mirrored these differences, with more understory Acer_c in older stands, and more understory Salix, Alnus, and Populus in younger stands on woody pioneer bars (Fig. b, d). The spatial distributions of overstory community types, geomorphic surface‐age classes, and riparian forest characteristics did not follow smooth, across‐valley lateral gradients; very few forest characteristics were correlated with distance from the channel. These vegetation patterns correspond to the concept of riparian ecosystems as a shifting habitat mosaic, in which patchy fluvial disturbance and floodplain turnover sustain an array of successional stages within the valley bottom (Stanford et al. ).

Linked geomorphic and vegetative responses to dams

Operation of two large dams on the Elwha River for 76–90 yr influenced the character of the riparian forest mosaic downstream. Differences in riparian vegetation between the three river segments were conditioned by differences in bottomland geomorphology likely related to effects of the dams (Kloehn et al. , East et al. ). Previous studies of fluvial processes and geomorphology on the same river segments showed that the dams induced channel armoring and reduced channel mobility and floodplain turnover, particularly in the middle segment and to a lesser extent in the lower segment, probably as a result of dramatically reduced sediment supply (Pohl , Kloehn et al. , Draut et al. ).

In the PCE, greater channel dynamics result in more frequent creation of bar surfaces and associated Alnus–Salix communities (O'Connor et al. , van Pelt et al. ). Accordingly, we found more pioneer bars, more extensive seedling recruitment areas, more riparian tree seedlings, and more overstory Salix in the upper segment, intermediate levels in the lower segment, and very low levels in the middle segment—consistent with the relative dynamism of the three stretches of the river. Reduction of geomorphic dynamism along dammed rivers has been observed to reduce woody pioneer plant recruitment on rivers around the world (Gonzalez et al. , Dixon et al. , Yan et al. ).

Reduced rates of bar creation or floodplain turnover also can lead to an overall older age structure of riparian forests downstream of dams (Merritt and Cooper , Gonzalez et al. , Johnson et al. ), including on the Elwha River. Kloehn et al. () reported a significant increase in older forest stands in the lower and middle segments, relative to the upper segment and another reference segment on the Quinault River. Similarly, we observed an older age structure in the middle segment, where Glines Canyon Dam apparently stabilized the channel and isolated floodplain forests, which then functionally became transitional fluvial terraces (sensu Latterell et al. ). Many of the older stands in the middle segment established prior to the completion of Glines Canyon Dam (62% compared to 28% and 16% of plots in the upper and lower segments, respectively), reflecting the near lack of floodplain turnover on this segment since dam emplacement. The resulting habitat mosaic in the middle segment still included an array of community types and geomorphic surfaces, but was dominated by fluvial terraces and established stands of late‐successional species.

While Kloehn et al. () observed older stands in both the lower and middle segments, we noted a mix of stand ages in the lower segment. The upstream‐most transect in the lower segment, nearest to Elwha Dam, had an older stand age structure, similar to those in the middle segment, whereas our downstream‐most transects included numerous floodplain stands (<51 yr) with an age structure more similar to the upper segment. This pattern fits with the serial discontinuity concept, which posits that dams interrupt longitudinal gradients in rivers and that those gradients recover with increasing distance downstream of dams (Ward and Stanford , Stanford and Ward ). Sediment capture by the two dams, combined with the narrow, less erodible valley just downstream of Elwha Dam, restricted sediment supply and channel mobility in upstream parts of the lower segment (Pohl ), thus reducing floodplain turnover and maintaining an older forest age structure. In contrast, increasing sediment supply from local erosion with distance downstream from Elwha Dam allowed the channel to remain relatively mobile in the wider valley near the river delta (Draut et al. ), resulting in floodplain turnover and younger riparian forests. We did not observe a similar pattern of decreasing dam effects with distance downstream in the middle segment, where sediment supply and valley morphology were consistent throughout.

In addition to effects of reduced sediment supply on channel mobility below the dams, another factor that likely influenced Elwha River geomorphic and vegetative responses to dam operations was the relative lack of large woody debris recruitment, log jam formation, and associated fluvial processes, including bar and island formation. In the PCE, large wood jams and subsequent sediment deposition patterns are significant drivers of bottomland geomorphology (Abbe and Montgomery , Gurnell et al. ). Although our work did not address this aspect of the system response, it was almost certainly a contributing factor and is the subject of ongoing research along the Elwha River (V. Leung, personal communication).

Other factors affecting riparian forest composition

Our work was aimed at relating changes in riparian forests to fluvial dynamism and hydrogeomorphic changes associated with dams. However, differences in numerous environmental factors, including climate, soils, fire, herbivory, and land management, also may have affected successional trajectories, species composition, and the character of the shifting habitat mosaic along the Elwha River (Mouw et al. , Kleindl et al. ). Besides the dams, human activities along the Elwha River have altered fire frequency (due to suppression efforts), large ungulate population size and behavior (due to predator control), and abundance of large conifers and woody debris (due to selective timber extraction), each with potential effects on riparian forest age structure and species composition.

The species composition of young geomorphic surface‐age classes on the Elwha River was consistent with other rivers in the region (Naiman et al. ), but the oldest stands we studied included relatively low frequency and abundance of some of the regional, zonal dominant conifers, particularly Picea and Tsuga, which occurred in only 2% and 13% of plots, respectively. Picea and Tsuga dominate mature fluvial terraces along rivers in the western part of the Olympic Peninsula (Fonda , van Pelt et al. ), but were also uncommon on the Elwha River in young (<50 years old) stands examined by Stolnack and Naiman (). Instead, Pseudotsuga and Abies were the most common and abundant conifers in our plots, occurring in ~35% of plots across a range of geomorphic surface‐age classes, similar to the stands examined by Stolnack and Naiman (). Picea and Tsuga require moister conditions than Pseudotsuga and may be hindered by the relatively dry climate of the Elwha River compared to the moister, western side of the Olympic Peninsula. Overstory Thuja also occurred in 20% of plots, mainly in stands >90 years old but sometimes in relatively young stands 34–55 years old, unlike the findings of Stolnack and Naiman (), because they did not examine forests along the lower segment where almost all of the Thuja we observed were located. Higher soil calcium and magnesium in the lower segment may have contributed to higher Thuja abundance (Hawkins and Robbins ; Perry et al., in press).

Recurring fires affect the distribution, age structure, and species composition of Pacific Northwest forests, including on the Olympic Peninsula, favoring Pseudotsuga over later‐successional Tsuga (Huff ). More frequent fires associated with the drier conditions of the northeastern side of the Olympic Peninsula (www.nps.gov/olym/learn/management/fire-history.htm) may have contributed to higher Pseudotsuga abundance, and lower Tsuga abundance, along the Elwha River than along rivers in the western side of the Olympic Peninsula. A study of fire history over a 440‐year period in our upper segment and surrounding uplands estimated a mean fire return interval of 127 yr (Wendel and Zabowski ). Some large historical fires extended across the bottomland, the most recent of which occurred in 1898. Fire suppression since the 1920s has reduced fire frequency and area burned along the Elwha River, although a period when large fires were rare also occurred in the 18th and 19th centuries (Wendel and Zabowski ). Fire can be a driver of the shifting habitat mosaic in some riparian ecosystems (Kleindl et al. ). Fewer and smaller recent fires on the Elwha River may have resulted in older forest age structures and later‐successional species composition along all three river segments compared to historical conditions.

Ungulate population size and behavior, particularly of elk, is thought to have changed beginning in the 1920s when wolves were all but eliminated from the Olympic Peninsula (Ratti et al. ). Increased population size and behavioral changes such as more local and sedentary herd movements have been associated with substantial declines in palatable riparian trees, including Populus balsamifera subsp. trichocarpa, Acer macrophyllum, and Tsuga heterophylla in the PCE and various Populus spp. elsewhere in western North America (Beschta and Ripple , , ). Ungulate herbivory, in combination with dam effects on channel mobility and floodplain turnover, may explain the differences in Populus abundance among river segments along the Elwha River. Populus seedlings and saplings were most abundant in the upper segment, where suitable surfaces for seedling recruitment were common. However, both overstory and understory Populus trees were relatively rare. Populus saplings in the upper segment were often browsed, and browsed individuals typically did not exceed 2 m in height, suggesting that ungulate herbivory may have limited Populus establishment to maturity along our upper segment transects. Populus did, however, contribute significantly to overstory basal area in a portion of the upper segment in another study (Acker et al. ). In contrast, in the middle segment, we observed almost no Populus seedlings or saplings, but this appeared to be a result of lack of gravel bars, surfaces that typically support Populus germination and establishment (Polzin and Rood ), due to reduced channel mobility downstream of Glines Canyon Dam. Of Populus overstory trees on the middle segment, 85% were older than ~40 yr (>30 cm dbh), suggesting that ungulate herbivory may have affected Populus less in the mid‐1900s, when Populus were establishing on pre‐dam surfaces abandoned after dam emplacement. Overstory Populus were most abundant in the lower segment, perhaps as a result of lower elk herbivory downstream of the Olympic National Park boundary, coupled with moderate channel mobility and availability of young geomorphic surfaces for Populus seedling establishment. Populus seedlings and saplings were present on pioneer bars and in the understory in the lower segment, and 63% of Populus overstory trees were younger than ~40 yr. Clonal expansion can also contribute to new Populus stands (Polzin and Rood ); however, we did not observe clear indications that this was common on the Elwha River.

Logging of the Elwha Valley and adjacent landscape was a primary land use in the late 1800s and much of the 1900s (Sadin et al. ), continuing until 1958 within the Olympic National Park boundary (Lien ), and until more recently downstream of the park. Both clearcuts and selective logging were common (Lien ). Pseudotsuga were widely harvested along the Elwha River throughout the late 1800s and 1900s, Tsuga were harvested for paper mills beginning in the 1920s, Picea were harvested for airplane production during World War I, and large woody debris was salvaged from river margins and Lake Mills in the mid‐1900s (Lien , Sadin et al. ). Forest clearcuts can increase flood frequency and magnitude in the Pacific Northwest (Zhang and Wei ). At the same time, logging can result in long‐term decreases in large woody debris (Murphy and Koski , Bilby and Ward ). Thus, changes to streamflow hydrology and wood budgets associated with logging may have influenced fluvial dynamics and creation of new geomorphic surfaces along the Elwha River. In addition, conifer logging can increase dominance by hardwoods in Pacific Northwest forests (Berg ), sometimes reducing conifer abundance for >100 yr (Banner and LePage ). Logging of Picea and Tsuga may have contributed to their relative scarcity along the Elwha River. Selective logging of conifers, coupled with lower channel mobility and floodplain turnover, also may have contributed to the high abundance of Acer_m and Acer_c in the middle segment.

Implications for dam removal

Dam removal is a river restoration approach that has increased dramatically in recent decades due mostly to aging infrastructure (O'Connor et al. ). Bottomland vegetation restoration is seldom the primary rationale (Shafroth et al. , Orr and Stanley ), but dam removal may at least partially restore riparian forest dynamics. As with riparian forest responses to dam operations, responses to the recent removal of the Glines Canyon and Elwha dams will depend largely on channel and geomorphic responses to the release, transport, and deposition of the large volume of sediment formerly stored in the reservoirs, as well as changes in large woody debris dynamics. A simple return to pre‐dam vegetation dynamics is not necessarily to be expected, because system adjustments in response to the substantial pulse of sediment and large wood could take decades, and because the system will be responding from a condition that had already been altered by decades of damming, fire suppression, increased ungulate herbivory, and logging. Sediment erosion and transport since dam removal have been fairly rapid thus far (Warrick et al. ), however, and increased rates of channel migration should lead to more sites (such as new bars) available for young forest establishment in the middle and lower segments than were present when the river was dammed (Kloehn et al. ), perhaps leading to forest dynamics more similar to the upper segment. Many new bars will likely be unstable for several years while sediment and wood fluxes are still high, but eventually some surfaces should facilitate new forest development, as sediment flux declines and some wood jams stabilize. In the first year following dam removal, most sediment was transported to the Strait of Juan de Fuca, with relatively little deposition on floodplains (Warrick et al. ). Some new bar creation was observed in the middle and lower segments, with pioneer woody plants germinating on many of these surfaces (East et al. ). Significant volumes of sediment are still available for transport as of this writing, and more sediment deposition on floodplains and in floodplain channels is expected in future years, which could also drive vegetation responses. Understory plant communities could be buried, with compositional shifts occurring where a different mix of species respond favorably to burial. Deep sediment deposits could contribute to mortality of some woody plants as well, potentially resulting in patches with somewhat novel plant communities (Shafroth et al. ).

Acknowledgments

Financial support for this work was provided through the U.S. Geological Survey's Park Oriented Biological Support Program. Thanks go to staff at the National Park Service, Olympic National Park, for facilitating access to study sites in the middle and upper river segments, especially Jerry Freilich and Steve Acker. For access to the lower segment sites, thanks go to the Lower Elwha Klallam Tribe, particularly Mike McHenry and Mel Elofson. Several people helped with field data collection including Ann Marie Casey, Robin Jenkinson, Greg Auble, and Sonny Sampson. S. Acker and two anonymous reviewers provided constructive suggestions for improving the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.