Linking Marine Ecosystem Response to Shoreline Armor Removal and Large Dam Removals in the Elwha River and Nearshore, Washington, USA

Large dams are well documented to have wide-scale ecosystem effects, and, as such, large dam removals (DRs) have become an essential ecosystem restoration tool that can have direct influence on the marine nearshore zone (Ezcurra et al., 2019; Maavara et al., 2020). However, DR plans often give little consideration to nearshore ecosystem function, which can lead to omission of basic nearshore functional elements that are critical to ecosystem restoration associated with large DR and, ultimately, incomplete ecosystem restoration (Shaffer et al., 2008, 2009, 2017a,b).

Similarly, ecosystem impairment from shoreline armor, including rock or concrete, installed to protect shorelines from erosion and ecosystem restoration through shoreline armor removal (AR), has been well documented (Dethier et al., 2016; Heerhartz et al., 2014, 2016; Lee et al., 2018; Toft et al., 2013). However, no information exists on the relationship of large DRs, shoreline armor and removal, and associated shoreline restoration of armored and unarmored beaches affected by in-river dams and DRs. Globally, shoreline armoring continues to occur, whereas large DRs, often at the same time, continue to be implemented with the goal of restoring coastal watershed ecosystems (Boucher, 2022; Dam Removal Europe, 2022). Therefore, understanding the relationship of the shoreline and riverine ecosystem response to restoration actions provides new and especially valuable information for nearshore ecosystem function and management, including large DR restoration priorities.

Located on the north central Olympic Peninsula in Washington State, the Elwha River had two hydroelectric dams installed in the main river channel more than 100 years ago. Most of the Elwha River is located within Olympic National Park. The Elwha primary drift cell comprises approximately 21 km of shoreline that extends from the west end of Freshwater Bay east to the tip of Ediz Hook (Shaffer et al., 2008). The Elwha River DR project, which included the removal of the 32 m high Elwha and 64 m high Glines Canyon dams that subsequently liberated approximately 18 Mt (∼9 million m³) of sediment to the watershed, was the world’s largest completed DR restoration project to date (Warrick et al., 2019). The DR began in 2011 and ended in 2014.

Approximately one-half the liberated sediment was delivered to the coast within 5 years of DR and resulted in abrupt changes to the coastal system (Warrick et al., 2019). During these dramatic transitions in sediment supply and delivery to the coastal zone, the unarmored reach of the Elwha nearshore responded almost immediately (Shaffer et al., 2017b). In contrast, most of the armored Elwha nearshore continued to be highly erosional, likely because of remaining alterations in the lower river and shoreline armoring placed along the east delta, feeder bluffs, and spit of the Elwha drift cell over the last 100 years (Parks, 2015; Shaffer et al., 2008, 2017b). To address this persistent ecosystem impediment, a large-scale shoreline restoration project was undertaken. From 2016–18, the armored portion of the east Elwha Delta reach was restored through the removal of 4700 m³ of large riprap (shoreline armor) from more than 650 m of the Elwha nearshore shoreline. This paper compares and quantifies the nearshore ecosystem response to restoration actions of DR, shoreline AR, and the cumulative results of both restoration actions.

Standard nearshore physical metrics of beach profile topography, beach sediment grain size, large woody debris (LWD) composition and distribution, and ecological metrics of beach wrack composition, including invertebrate species composition and abundance, were used to define nearshore physical restoration action response to individual actions of DR and shoreline restoration (Rich et al. 2014; Toft, Litle, and Adams, 2015). The metric responses of the two sites that experienced DR (one with armoring, one without armoring) were then compared with the control site (no DR, no AR) and analyzed with each restoration action to define the relative role of DRs, shoreline AR, and the two actions combined to achieve nearshore ecosystem restoration. These metrics were used to test two hypotheses: (1) that DR alone does not restore the ecological metrics when shoreline impediments of shoreline armoring remain and (2) when shoreline armor is removed, DR shorelines physically and ecologically fully restore relative to unarmored shorelines.

The study area is located along the nearshore of the central Strait of Juan de Fuca, including the Elwha and Dungeness River drift cells, on the northern Olympic Peninsula in Washington State (Figure 1). In the early 1900s, two dams were constructed on the Elwha River, and the shoreline of the eastern Elwha Delta and drift cell was heavily armored repeatedly until the early 2000s, nearly eliminating river conveyance of sediment, which previously sustained the beaches in the drift cell (Shaffer et al., 2008; Ward et al., 2008). Impacts resulting from in river dams and shoreline armoring included increased shoreline erosion, decreased habitat extent and complexity, and decreased connectivity of the shoreline and the lower river estuary due to sediment starvation and extensive shoreline armoring (Parks, 2015; Parks, Shaffer, and Barry, 2013). A century later, the two in-river dams were removed with the intent of a complete watershed scale ecosystem recovery; however, the shoreline armoring impediments remained along much of the Elwha drift cell following DRs (Shaffer et al., 2017a).

This study included three treatments and three phases (Table 1). Two of the three treatments were restoration actions that were within the Elwha DR drift cell (one with and one without armoring). The restoration treatment shorelines of the Elwha shoreline are described in further detail in Parks (2015), Parks, Shaffer and Barry (2013), and Warrick et al. (2019). The third treatment is a control site on the Dungeness drift cell, well outside the Elwha drift cell, which served as a comparative treatment unaffected by either the dam or shoreline armor (Figure 1). These nearshore treatments were classified as follows.

This treatment experienced only DR response and is the unarmored area of the west shoreline located in the Freshwater Bay area of the Elwha drift cell (Figure 1). It is an unarmored embayed 2000-mm shoreline to the west of the Elwha River mouth. This stretch of shoreline experienced persistent erosion because of sediment starvation caused by the Elwha and Glines in-river dams and significant sediment delivery during DR (Parks, 2015, Shaffer et al., 2017a,b; Warrick et al., 2019). This segment of the Elwha nearshore was monitored for beach profiles and LWD before, during, and after the DR project and also for beach profiles, LWD, beach wrack, and sediment composition after the shoreline AR project.

This treatment experienced both dam and armor removal (DAR) response. This is the armored area of the east Elwha shoreline also known as the east delta. This reach of shoreline experienced persistent erosion due to sediment starvation from in-river dams and significant shoreline armor over the last century, including during and after DRs (Warrick et al., 2019). To address this persistent ecosystem impediment, a large-scale shoreline restoration project was undertaken. The armored portion of the east Elwha Delta reach was restored through the removal of 4700 m³ of riprap from more than 650 m of the Elwha nearshore shoreline. In this study, the nearshore ecosystem response to restoration actions of DR, shoreline AR, and the cumulative results of both restoration actions are defined.

This treatment experienced neither the restoration action of DR nor shoreline AR. The control located in the Dungeness Bay region of the Dungeness drift cell is 31 km east of the Elwha drift cell and included a 600-m long shoreline that has had no armoring and no alterations over the last 20 years. The treatment has the same NW orientation and experiences similar fetch as the Elwha drift cell study treatments but is somewhat protected to the north by Dungeness Spit. It is a county park and is protected from development. The control is ecologically dissimilar in some ways from the other study treatment locations (i.e. shorter, somewhat protected from dominant wind and fetch by a spit, in a different drift cell). This paper primarily focuses on the differences in relative changes of the treatments related to restoration actions and phases, so the control treatment is an appropriate comparative area to represent physical and ecological variability in a region where no restoration actions are occurring. This area was monitored for beach topography, sediment composition, LWD, beach wrack composition, and invertebrate composition both before and after the AR restoration project.

Standard nearshore physical metrics of beach profile topography, beach sediment grain-size, LWD composition and distribution, and ecological metrics of beach wrack composition were used, including invertebrate species composition and abundance (Toft, Litle, and Adams, 2015), to define nearshore physical restoration action response to individual actions of DR and shoreline restoration. The metric responses of the two sites that experienced DR (one with armoring, one without armoring) were compared with the control site (no DR, no AR) and analyzed with each restoration action to define the relative role of DRs, shoreline AR, and the two actions combined to achieve nearshore ecosystem restoration. These metrics were used to test the two hypotheses: (1) DR alone does not restore the ecological metrics when shoreline impediments of shoreline armoring remain, and (2) when shoreline armor is removed, DR shorelines physically and ecologically fully restore relative to unarmored shorelines.

The study was divided into three phases: pre-dam removal (PreDR); post-dam removal and pre-armor removal PostDRPreAR; and post-dam removal and post-armor removal (PostDRPostAR; Table 1).

Beach profile data included beach shore face topography (elevation in meters) and beach shore face, sediment mean, intermediate-axis, grain-sized diameter in millimeters. Beach topography was measured normal to the alongshore direction in a series of topographic profiles extending from the upper-beach berm or bluff toe, 30 m down to the low-tide water line, and were spaced approximately 30 m apart parallel to one another. Transects were truncated to 30 m in length to normalize transect length across the study area. Measurements of vertical elevation of the beach surface were recorded every meter along the topographic profile.

Topographic beach profiles were measured within the study area from a variety of sources. These included aerial fixed-wing LIDAR digital elevation models (DEMs; Dewberry, Inc., 2016; Entrix, Inc., 2009; Geomatics Data Solutions, 2015; Quantum Spatial, 2019; Woolpert, Inc., 2013, 2014, 2015), real-time kinematic GPS (RTK-GPS) data collected by the U.S. Geological Survey in 2009 and 2016 (Stevens et al., 2016; Warrick et al., 2007), and RTK-GPS data collected by the authors in 2017, 2018, and 2019. Beach profiles from 2009, 2012, 2014, and 2015 were extracted from LIDAR-based DEMs using the ARC GIS 10.6.1 3D extension tool (ESRI, 2020).

Beach topography was measured using foot surveys in 2017, 2018, and 2019, with an Arrow Gold™ RTK-GPS mounted on either a standard 2-m survey pole or backpack and connected to the Washington State Reference Network to receive real-time horizontal and vertical position corrections. Horizontal survey data were collected in World Geodetic System 1984 coordinates in meters, and vertical data were collected using the National Vertical Datum of 1988 meters. Horizontal and vertical accuracy of RTK-GPS measurements was assessed by occupying survey benchmarks within the sample area. Overall root-mean-square accuracy of DEMs and RTK-GPS measurements ranged from ±5.9 to ±10 cm.

The beach wrack sampling method was used, which was developed by Heerhartz et al. (2014) and is employed by the Puget Sound Ecosystem Monitoring Program described in Toft et al. (2017). For both biological and sediment sampling, two 50-m transects were placed parallel to the shoreline along the high-tide line with fresh wrack deposition at each of the three beaches. Transect locations were selected randomly, and at 10 random points along each transect, a 30 cm × 30 cm quadrat was placed for sampling. To determine beach wrack composition, within each quadrat a visual estimate of the percentage of composition of algae, eelgrass, terrestrial plant material, and human debris was conducted for each point. A minimum of three observers quantified beach wrack composition per quadrat, with the same person making observations at all quadrats. Observations of all observers were averaged for each quadrat.

To visually estimate sediment grain size, two random points along the transect were randomly selected, and the 30 cm × 30 cm quadrat was again placed on the transect line. The percentage of each sediment substrate category was estimated at the surface and then again after digging down 5 cm into the sediment (Table 2). These two percentage estimates were then averaged to gain a total percentage of each substrate category. Sediment substrate categories were identified as per Heerharz et al. (2014): cobble (>6 cm), pebble (4–6 cm), granule (2–4 mm), sand (“gritty” up to 2 mm), and silt/clay (smooth; Table 2).

Beach wrack invertebrates were collected from five randomly selected points of the 10 points used for composition estimates. At each point, a 15-cm diameter corer made from PVC pipe was pushed from the surface of the wrack to the beach surface. The wrack samples were bagged, labeled with sample location information (i.e. sample number, GPS coordinates, location along transect), and frozen until processing. In the laboratory, invertebrate samples were identified by trained invertebrate specialists to lowest taxa using standard dissecting and microscope techniques and counted per sample.

The LWD monitoring protocols developed by Rich et al. (2014) were used. A series of 18-m-wide transects were established at permanent locations parallel to the shore and then sampled from the bank to the waterline. Three transects each were sampled at the DR and DAR treatment; two transects were sampled at the control treatment. Transects were sampled once a year in either June or July 2016–19.

Data were collected within a 4-hour period around the low tide of the day. Length and width of all LWD greater than 0.5 m in length and 10 cm in diameter within the transect line were recorded. All LWD pieces less than 0.5 m in length and diameter less than 10 cm were estimated for percentage of coverage of total transect area for each transect.

A PERMANOVA was used to examine the effect of phase and treatment and also whether an interaction between phase and treatment occurs on beach wrack composition and beach wrack sediment categories. The adonis2 function from the vegan package (version 2.6-4) was used, and data were displayed with a nonmetric multidimensional scaling plot and Bray-Curtis dissimilarities (Oksanen et al., 2022). Beach wrack sediment was analyzed separately for the top layer and 5-cm depth, as well as the averaged total sediment categories.

Generalized additive models were used to model beach wrack invertebrate total abundance, Hill-Richness, and Hill-Shannon metrics. Hill-Richness and Hill-Shannon metrics were calculated using the hillR package (version 0.5.2; Li, 2018). All models included phase, treatment, an interaction between phase and treatment as fixed effects, and year as a random effect. Each metric (invertebrate total abundance, Hill-Richness, and Hill-Shannon) had two models conducted: total average percentage of wrack cover and average percentage of algae cover. The Akaike Information Criterion model selection was conducted for each metric to determine which model best fit each metric. Negative binomial distributions were used to model invertebrate total abundance, whereas a scat distribution was used for Hill-Richness and a gamma distribution for the Hill-Shannon models. The final comparison between the DAR and DR treatments was determined using a Wilcoxon rank sum test.

A generalized linear mixed model with a negative binomial distribution was used to examine the effect of phase and treatment and to determine whether an interaction between phase and treatment on LWD counts occurred. Additionally, linear mixed effects models were used to examine the effect of phase and treatment and to determine whether an interaction occurs between phase and treatment on log-transformed LWD length and diameter. All three models included year as a random effect. Pairwise contrasts were examined for each model using the glht function from the multcomp package (version 1.4-25; Hothorn, Bretz, and Westfall, 2008).

Assumption checks for models were conducted visually, and all data analyses were conducted in R (version 4.4.0; R Core Team, 2024).

Beach profiles of all treatments changed with phase. Lower elevations of the profiles showed the most change. The DR treatment with no armor changed positively overall by an average of 0.3 m PostDR and 0.4 m after the AR period (even though AR occurred outside this part of the study area); it had a net change of 0.7 m relative to PreDR, with the majority of the change occurring within the waterward 15–30 m horizontal distance of the beach profile (Figure 2).

In contrast, the DAR treatment site continued to erode after DR and before DAR, losing on average 0.5 m after the dams were removed but before armor was removed. After AR, this shoreline aggraded, gaining 0.7 m in elevation relative to PreAR profile. Net gain of the armored shoreline after armor removal was therefore 0.2 m relative to preDR (Figure 2).

Beach wrack vegetation communities differed based on treatment (PERMANOVA; F_2,231 = 30.9; p = 0.001), phase (PERMANOVA; F_1,231 = 8.4; p = 0.001), and an interaction between treatment and phase (PERMANOVA; F_2,231 = 10.4; p = 0.001; Figures 3 and 4). PostDrPreAR phases were more positively correlated with algae and human debris, whereas PostDRPostAR phases were more positively correlated with eelgrass and terrestrial plants in the wrack composition (Figure 3). Control and DAR treatments were also more positively correlated with terrestrial and eelgrass than DR treatments (Figure 3). The DR treatments tended to have higher wrack cover than DAR treatments in both phases, similar total wrack cover as the control during the PostDRPreAR, and more than the control during PostDRPostAR (Figure 4).

Beach wrack sediment communities in the total sediment samples differed based on treatment (PERMANOVA; F_2,66 = 30.0; p = 0.001), phase (PERMANOVA; F_1,66 = 4.8; p = 0.015), and an interaction between treatment and phase (PERMANOVA; F_2,66 = 5.1; p = 0.002; Figures 5 and 6). PostDRPreAR phase for the DAR and DR treatments was more positively correlated with sand, whereas the control treatment is more positively correlated with cobble and pebbles (Figure 5). In the PostDRPostAR phase, the DAR was more positively correlated with sand, whereas the DR and control treatments were more positively correlated with gravel and pebbles (Figure 5). Sand was the most common sediment type, especially during the PostDRPreAR phase (Figure 6). The top layer and 5-cm-depth sediment samples showed similar results to the total sediment samples (Figures 6, S1, and S2; Tables S2 and S3).

The most commonly observed taxa were Amphipoda and Oligochaeta, followed by Acari and Nematodes (Figures 7 and 8). Abundance of most abundant taxa were lower at the DAR site and increased yearly after AR (Figures 7 and 8). Individual taxa responses to the AR event were varied and were different among restoration actions.

The best model for invertebrate total abundance included total percentage of algae, whereas best models for invertebrate richness and Shannon diversity metrics included total wrack percentage (Table S6). For total abundance models, no main effect of phase occurred; however, a main effect of treatment was found, with the DAR treatment having lower invertebrate abundance than the control treatment. Additionally, an interaction between phase and treatment occurred (Table 3, Figure 9).

For the best Hill-Richness model, higher richness at the control treatment was found, as compared with the DAR treatment, and no other treatments or phases were significantly different (Table 3, Figure 9). A main effect of total wrack cover on Hill-Richness occurred, with richness declining to the lowest, around 30% wrack cover, before increasing with wrack cover (Figure S3a).

For the best Hill-Shannon model, higher evenness was found at the DR treatment compared with the control; however, no other treatments or phases were significantly different (Table 3, Figure 9). A main effect of total wrack cover was found on the Hill-Shannon metric, with evenness increasing to the maximum around 30% wrack cover before decreasing with wrack cover (Figure S3b).

The LWD count increased in both the PostDRPreAR and PostDRPostAR compared with PreDR; however, no difference occurred between both PostDR phases (Table 4, Figure 10). The DR had the highest LWD count, followed by DAR, with the control treatment having the lowest LWD count regardless of phase (Table 4, Figure 10).

On average there were longer LWD logs in the PreDR phase compared with the PostDRPostAR as well as longer logs in the DR treatment compared with the DAR treatment (Table 4, Figure 10). All other phase and treatments did not differ from each other significantly (Table 4). Additionally, a significant interaction occurred between treatment and phase (Table 4).

Finally, for LWD mean log diameter, the PreDR phase had larger diameter logs than both PostDR phases, although the PostDR phases did not differ from each other (Table 4, Figure 10). No difference occurred between the control treatment and the DAR or DR treatment; however, larger diameter logs in the DR treatment occurred, as compared with the DAR treatment (Table 4, Figure 10). A significant interaction also occurred between phase and treatment (Table 4).

This work provides a novel study of the role that shoreline armoring plays in DR shorelines and illustrates that in the Elwha system, both dam and shoreline ARs contribute to restore a degraded shoreline. The study also indicates that when dams are removed but shoreline armor remains, armored shorelines within the drift cell continue to function at an impaired state until armor is removed. After the armor is removed, the physical and ecological metrics of the beach respond quickly with restoration trends consistent with, but delayed relative to, changes observed at unarmored shorelines.

Earlier work has documented that shoreline armoring and in-river dams work synergistically to impair shorelines. Parks (2015) documented that when the nearshore is sediment starved because of in-river dams, armored feeder bluff shorelines are coarser, and grain size is larger than unarmored feeder bluff shorelines. This study illustrated that armored shorelines that experienced DR responded and rapidly became persistently broad and finer grained—but only after armoring was removed. This response was not only seen at the AR treatment locations but also farther along the shoreline (Figures 2 and 11; Parks, 2022), illustrating that shoreline armor and its removal can have both a local and drift cell–wide ecosystem effect. This is a key consideration for broader geographic implications for both shoreline armoring and shoreline restoration through AR.

When shoreline armor remains down drift of large sediment restoration events, the armoring will continue to impede ecosystem (functions) restoration, even after significant sediment additions and farther afield than the shoreline armored treatment. When the shoreline armoring is removed, ecological response is varied, indicating shoreline ecological response may take longer than the timeline of this study. Finally, these changes to the physical system are also being reflected to higher trophic systems. CWI et al. (2022) documented that the distribution of surf smelt spawn expanded exactly concurrent with or mirrored sediment changes associated with dam and then shoreline AR along the drift cell.

Shoreline armoring effects are valuable ecological factors to shorelines (Dethier et al., 2016; Heerhartz et al., 2014, 2016; Lee et al., 2018) but to date have been overlooked when quantifying hydrodynamic response of shorelines affected by DR (Warrick et al., 2019; Zurbuchen et al., 2020). Our work confirms the value of including shoreline AR into coastal planning for DRs and the significance of proper geographic scale in understanding and optimizing the ecological response to them. Lack of replication of all work in the Elwha system is a strong motivator in replicating these studies along other DR shorelines.

Nearshore ecosystem function is complex and synergistic, and it varies with restoration actions. Prior work in the Elwha nearshore has demonstrated that intact, mature nearshore habitats are functionally more stable and resilient than newly restored habitats after large-scale DRs (Shaffer, Munsch, and Juanes, 2018). Further, nearshore ecosystem function and restoration have critical seasonal and inter-annual temporal dimensions (Shaffer et al., 2012), and newly created habitats may take decades to establish, stabilize, and become functionally mature. Finally, many of the species driving restoration action have complex and long life histories that depend on multiple nearshore habitats seasonally that collectively make ecosystem restoration a decades’ long endeavor. It is therefore clearly crucial to prioritize restoring impaired coastal zones at the proper ecosystem scale (drift cell) and temporal scales.

Both large dam and shoreline armoring have significant and interactive effects to nearshore ecosystem forming and functional processes. Restoration actions to remedy these impacts are key for ecosystem recovery; however, DR alone will not achieve complete nearshore ecosystem restoration when shoreline impediments remain along the DR shoreline. Shoreline AR can remedy this gap; it results in significant, dramatic, and synergistic restoration both at and far distal to the AR treatment with ecological restoration changes that are translated across many facets of nearshore ecosystem function.

A number of college collaborators (including faculty and students) and CWI interns, staff, and volunteers (including Jenise Bauman, Wesley Greentree, Kirsten Simonsen, Tara McBride, Sara Schoenemann, Katrina Campbell, David Harvey, Lindsey Howard, Tony Thompson, Breyanna Waldsmith, Seren Weber, and Curtis Welker) assisted with the analytical and field elements of this work. The Lower Elwha Klallam Tribe, Pam Lowry, the Dudley family, and Malcom and Phoebe Moore provided collaboration. Bruch and Bruch Construction and 2 Grade LLC provided partnership and removed the armor from the shoreline. Professional videographers Laura James and John Gussman provided filming services. Funding for this work was provided by Patagonia, Inc.; USFWS grant numbers F16AP000157, F17AC00393, and F18AP00149; and Washington Recreation and Conservation Office grant numbers 15-1045 and 16-2089 through the Salmon Recovery Funding Board and the Estuary and Salmon Restoration Program. Thank you all.

An extended summary of this research was originally published in Shaffer et al. (2024).