Future sea-level rise drives rocky intertidal habitat loss and benthic community change

Sea-level rise projections and potential impacts to the rocky intertidal zone

Sea-level rise is predicted to alter habitat availability and modify community structure in many marine ecosystems (Hoegh-Guldberg & Bruno, 2010; Nicholls & Cazenave, 2010). Long-term monitoring of sea-level shows considerable variability among locations, but an overall rising trend with recently increased rates of change at many locations (Church & White, 2006; Mcleod et al., 2010; Nicholls & Cazenave, 2010; Rahmstorf, 2010; National Research Council, 2012; Church et al., 2013; Cazenave et al., 2014; Chen et al., 2017). While estimates remain uncertain, recent studies project up to 2.5 m of SLR within the next century (Church et al., 2013; Sweet et al., 2017). Such a large magnitude and rapid shift poses a substantial risk to the integrity of coastal ecosystems, yet the extent to which SLR will modify the physical and ecological structure of rocky coastlines remains mostly unknown (Denny & Paine, 1998; Thompson, Crowe & Hawkins, 2002; Harley et al., 2006; Helmuth et al., 2006).

Certain characteristics of rocky intertidal systems make them particularly vulnerable to SLR driven habitat loss (Denny & Paine, 1998; Thompson, Crowe & Hawkins, 2002; Harley et al., 2006; Helmuth et al., 2006). When backed by steep cliffs or anthropogenic structures, as is the case on the majority of coastlines globally (Emery & Kuhn, 1982), many rocky shores are expected to experience “coastal squeeze”—a general narrowing of the extent of the intertidal zone (Doody, 2013; Pontee, 2013) and a steepening of the coastal profile (Vaselli et al., 2008). Whether SLR will lead to losses or gains in habitat at a given site is highly dependent upon geomorphology, rates of erosion, underlying stratigraphy, and alteration in sedimentation processes at local and regional scales (Emery & Kuhn, 1980; Young & Ashford, 2006; Masselink et al., 2020; Schaefer et al., 2020). Studies investigating durable coastlines with low erosion rates have suggested that SLR may cause substantial rocky intertidal habitat loss. Estimates of habitat loss range from 10–27% and 26–50% with 0.30 m and 1.90 m of SLR, respectively, in northern Scotland (Jackson & Mcilvenny, 2011), to 10% and 57% with 1.0 and 2.0 m of SLR, respectively, in Oregon, USA (Hollenbeck, Olsen & Haig, 2014). Thorner, Kumar & Smith (2014) found that with 0.3–1.0 m of SLR on five rocky headlands in New South Wales, Australia, impacts will be variable but will largely result in substantial habitat loss. Most recently, Schaefer et al. (2020) suggested that along 50 km of rocky shoreline in eastern Australia, SLR will generally reduce habitat availability, though impacts will show high spatial variability, with gently sloped shores that lack the ability to migrate landward experiencing the greatest losses of habitat. Because habitat loss is among the greatest threats to global biodiversity (Brooks et al., 2002; Mantyka-Pringle, Martin & Rhodes, 2012), it is critical to evaluate the current state of rocky intertidal ecosystems and assess how SLR may affect these important communities in coming years.

The unique ecological characteristics of rocky intertidal systems may also make them particularly vulnerable to changes in ecological structure and function as a result of SLR. The rocky intertidal zone is characterized by patterns of ecological zonation that manifest as distinct, and sometimes strict, bands along the tidal elevation gradient. These bands are generated by a variety of spatially variable and density-dependent biological mechanisms, such as competition (Connell, 1961; Menge, 1976), mutualism (Menge, 1995; Bertness & Leonard, 1997), predation (Paine, 1969; Paine, 1974; Paine, 1980; Menge, 1976), and differences in physiological tolerances to extreme environmental conditions (Connell, 1972; Paine, 1974). Regardless of the potential for the concomitant provisioning of suitable space via erosion, SLR will cause an upward shift in this banding as the current intertidal zone is submerged. It is unclear whether the changing intertidal zone will be suitable for colonization by intertidal species, as habitat characteristics, such as slope (Vaselli et al., 2008) and substrate type (Thorner, Kumar & Smith, 2014), and physical environmental conditions, such as wave stress and exposure (Harley et al., 2006), may also change. Thus, coastal squeeze and the rapid upward shift in intertidal area will likely impact the abundance, distribution, and competitive interactions of rocky intertidal species.

Approaches to quantifying climate change impacts in the rocky intertidal zone

Rocky intertidal habitats are among the most extensively studied marine ecosystems, and are well known for their rich and varied history of influential research (Ricketts & Calvin, 1939; Connell, 1961; Lewis, 1964; Paine, 1966; Paine, 1974; Paine, 1994; Dayton, 1971; Menge, 1976; Sousa, 1979; Paine & Levin, 1981; Underwood & Denley, 1984; Gaines & Roughgarden, 1985; Roughgarden, Gaines & Possingham, 1988; Underwood, Chapman & Connell, 2000, to list a few). Despite much advancement in the field, a combination of factors has constrained field sampling designs (Andrew & Mapstone, 1987). Many studies have been tradition bound, selecting sampling units a priori rather than empirically evaluating the appropriate grain size for the particular question of interest (Andrew & Mapstone, 1987; Wiens, 1989; Underwood, 2000; Denny et al., 2004). Logistical challenges, such as difficulty of access and a brief sampling window constrained to periods of emersion, have necessitated labor-intensive in-situ surveys, which are limited in spatial extent (Blakeway et al., 2003). While remote sensing approaches have been used to collect data on larger spatial scales, technological limitations have somewhat restricted the resolution and deployment of these approaches. Such limitations have previously been a barrier to determining how SLR will modify ecological patterns and processes in rocky intertidal communities at ecologically relevant spatial scales (Runting, Wilson & Rhodes, 2013).

Recently, remote-sensing approaches and geographic information system (GIS) analysis software have been used to provide high-resolution, landscape-scale ecological information in the rocky intertidal zone (Guichard, Bourget & Agnard, 2000; Blakeway et al., 2003; Bryson et al., 2013; Thorner, Kumar & Smith, 2014; Gomes et al., 2018; Konar & Iken, 2018; reviewed by Garza, 2019). Importantly, when combined with detailed in-situ information, these approaches now provide researchers the opportunity to investigate ecological dynamics at scales commensurate with the landscape processes affecting these communities (Murfitt et al., 2017; Gomes et al., 2018; Garza, 2019).

We investigated the potential ecological impacts of future SLR on rocky intertidal ecosystems, using a multi-scale approach at two marine reserves in San Diego, CA, USA as a case study. Light Detection and Ranging (LiDAR) data were used to estimate current intertidal habitat extent at four study sites and site-level habitat area changes under a range of SLR scenarios. Using high-resolution in-situ imaging, we mapped intertidal habitats and quantified the percent cover of dominant benthic space occupiers and the density of sessile and mobile invertebrates across tidal elevations in 180 m² plots at each site. We then used a model-based approach to investigate future SLR driven changes in the areal extent of dominant benthic space occupiers and the numerical abundance of focal rocky intertidal invertebrate taxa. This work takes a critical step toward determining the future impact of SLR on rocky intertidal communities at ecologically relevant scales, and provides a framework for future monitoring and experimental efforts.

Survey locations and sites

Two rocky intertidal study locations in San Diego County were chosen: the Scripps Coastal Reserve (SCR) and the Cabrillo National Monument (CNM) (Fig. 1). Within these locations, one site was sampled at SCR (SCR 0), and three sites were sampled at CNM (CNM 1-3). These locations are recognized for their ecological and economic importance, and are both designated as marine protected areas (MPAs) under California state legislation. Regular long-term ecological monitoring by the Multi-Agency Rocky Intertidal Network (MARINe) has occurred at one site at SCR since 1997, and at three sites at CNM since 1990 (UCSC, 2020). The SCR site was studied under a permit granted by the Scripps Coastal Reserve Manager (application #: 33783), and CNM sites were studied under a permit granted by the US Department of the Interior National Park Service, Cabrillo (permit #: CABR-2016-SCI-0007).

As previously described by MARINe, the intertidal habitat at SCR is comprised of a gently sloping boulder-field backed by a large metamorphic dike running south by southwest through the site. The dike provides a distinct upper limit to the intertidal zone, which is immediately adjacent to steep cliffs that rise vertically to an elevation of approximately 75 m. CNM2 and CNM3 are characterized by broad and gently sloping sandstone benches, while CNM1 is comprised of large scattered boulders atop the same sandstone base found at the other CNM sites. At CNM, the intertidal zone is backed by low cliffs, which rise steeply to approximately 5–10 m elevation, then more gradually to a higher ridge. Cliff stratigraphy in both locations is comprised of lithified mudstone, siltstone, shale, and sandstone, and these cliffs regularly experience failure and variable rates of net yearly erosion (Emery & Kuhn, 1980; Benumof et al., 2000; Young & Ashford, 2006).

SLR scenarios and intertidal area estimation

In the Fifth Assessment Report, the United Nations Intergovernmental Panel on Climate Change (IPCC) projected a rise in global sea-level of between 0.26 m and 0.98 m by 2100 (Church et al., 2013). More recently the National Oceanic and Atmospheric Administration (NOAA) released projections that include SLR extremes of up to 2.5 m on US coastlines (Sweet et al., 2017). Sea-level rise projections for California generally fall within the range of global projections. Cayan et al. (2008) forecasted SLR on the coast of California of 0.11–0.72 m by 2070–2099. The National Research Council (NRC) projected 0.42–1.67 m of SLR by 2100 for the coast of California south of Cape Mendocino (National Research Council, 2012). For this study, SLR scenarios from 0.10 to 2.0 m were analyzed in 10.0 cm increments (twenty scenarios) to cover the generally accepted potential range of SLR for the California coast in the next century. By analyzing SLR in increments, we also avoided using projections for specific dates, and thus our analysis is not time specific and is more flexible to the uncertainty of the projections.

The total area of the intertidal zone at each site currently, and under each SLR scenario, was estimated from the 2009–2011 California Coastal Conservancy Coastal LiDAR Project: Hydro Flattened Bare Earth Digital Elevation Model (“NOAA, 2020”), which is hereafter referred to as the site-scale DEM or SDEM. WGS1984 and MLLW were selected as horizontal and vertical datums to allow later embedding of orthophotomosaics and direct referencing to common tidal datums. Tidal datums from NOAA tide stations nearest our survey sites were selected as the initial vertical boundaries of each site (Table S1, Fig. 1). Mean Lower Low Water (MLLW, 0 m) was used as the lower boundary of the intertidal zone because it is a commonly used datum and the penetration limitations of LiDAR do not allow accurate data below MLLW. Highest astronomical tide (HAT) was chosen as the upper boundary of the intertidal zone to encompass the zone that is only covered during the highest high tides and the spray zone. Straight lines running perpendicular to the elevation contours through the northernmost and southernmost long-term monitoring plots at each site were established as fixed upcoast and downcoast boundaries (Fig. 1). The DEM Surface Tools Extension for ArcGIS (ESRI 2016. ArcGIS Desktop: Release 10.5 Redlands, CA, USA) was used to provide surface area estimates on a grid by grid basis for the SDEM (Jenness, 2004). The Surface Tools Extension accounts for vertical relief values in area estimates by using adjusted surface area values for each SDEM grid cell rather than strictly planimetric area (Jenness, 2004).

As future erosion rates are difficult to estimate, we chose to follow a passive flooding approach. In each SLR scenario at each site, total intertidal area was estimated by adjusting the vertical boundaries to the future intertidal extent (current datums + SLR), then summing surface area of SDEM grid cells within the future intertidal extent. This allowed migration of the realized available habitat upward and shoreward under each SLR scenario, so that predicted SLR was applied to the existing coastal topography of our sites in our modeling (e.g., Vaselli et al., 2008). Thus, we prevented artificial constraint on the upper limit of the intertidal zone during future scenarios. These surface area values were also subset to allow analyses by intertidal zone (e.g., lower, middle, and upper, Table S1).

Large-area imaging of intertidal survey plots and image processing

In order to characterize the composition of the biological communities at these sites, we used an imaging approach referred to as large-area imaging, which produces spatially accurate, detailed maps of the benthos from digital imagery collected in field plots (Murfitt et al., 2017; Edwards et al., 2017). The general large-area imaging workflow involves the use of commercially available Structure-from-Motion photogrammetry (SfM) software (Agisoft PhotoScan Professional V.1.3 software (Agisoft LLC 2014, St. Petersburg, Russia)). This software allows creation of composite 3-dimensional reconstructions, also known as 3D models, from raw imagery. We then generate plot-scale Digital Elevation Models (PDEM) and orthorectified 2-dimensional top-down views, or orthophotomosaics, of the 3D models (Fig. 2). These 2D models are then scaled and uploaded into ArcGIS for ecological analysis and extraction of landscape-scale metrics.

At each site, imaging plots were positioned to span the maximal distance across the intertidal zone while also overlapping a subset of the long-term study sites. To establish plot locations, ten random coordinates were first generated in ArcGIS within the upper zone of the SDEM at each site. Coordinates were then ground-truthed to ensure they occurred on natural substrate within representative upper intertidal habitat. One coordinate was then randomly selected from the coordinates meeting these criteria as the upper, upcoast corner for the plot at each site. Rectangular 6.0 m × 30.0 m plots were then established running perpendicular to the shoreline along the elevation gradient. Plots of this size were deemed sufficient for accomplishing project goals while also balancing image and data processing capabilities.

Imagery was collected at all sites between December 2016 and January 2017, concurrent with the season’s lowest tides and with existing long-term monitoring surveys. A GoPro Hero 5 camera (GoPro Inc., San Mateo, CA) was mounted to a frame on a handheld transect line and passed between two surveyors across the plot, with adjacent passes separated by 0.5 m. Images were collected every 0.5 s using a linear field of view setting with an equivalent focal length of 24–49 mm. The camera was held approximately 1.0 m above the substrate to maximize overlap of images while also ensuring sufficient image resolution for benthic community identifications. Ten scale bars of known length (0.5 m) were deployed throughout the plot as x, y, z spatial references, and ground control point (GCP) coordinates were collected at the upcoast end of each scale bar. Images were collected over a 9.0 m × 33.0 m area to ensure that the target plot was imaged with sufficient coverage (1.5 m buffer around perimeter) to avoid areas missing data within the study plot.

The details of 3D model and orthophotomosaic creation have been described in detail elsewhere (Burns et al., 2015; Naughton et al., 2015; Murfitt et al., 2017). Briefly, Agisoft was used to first align imagery, estimate camera positions and scene geometry, and produce a sparse cloud of points extracted from the imagery. These sparse clouds were then optimized to correct model geometry and minimize alignment error. They were then assigned a coordinate system and scale using GCP coordinates and lengths from the scale bars within each plot. Textured dense point clouds were then produced and subsequently meshed to provide PDEMs (1.0 cm nominal post spacing). Finally, top down orthophotomosaics (0.3–0.4 mm/pixel resolution) were produced to allow visual identification and quantification of benthic organisms (Fig. 2). This resolution allowed identification of organisms of sizes greater than or equal to approximately 1.0 cm in diameter, though taxonomic resolution varied depending on morphology. Orthophotomosaics and PDEMs were exported from Agisoft in a raster format and uploaded into ArcGIS for subsequent data extraction.

Biological data extraction

SLR impacts to areal extent and invertebrate abundances

In order to assess how changes in habitat area as a result of SLR may influence the rocky intertidal community at our sites, the total areal extent (m²) of dominant benthic space occupiers and the abundance of focal invertebrate taxa (#) was estimated for each site under each SLR scenario. These estimates were produced by extrapolating the estimates of percent cover and density from the orthophotomosaics and PDEMs to the more spatially expansive SDEMs. The total areal extent (m²) of each dominant benthic space occupier in each SLR scenario at each site was calculated by multiplying the percent cover of each species in each bin by the total area in the respective bin as calculated from the SDEM, and then summing across bins. This metric provided an assessment of how SLR will influence absolute areal extent of each taxon and allowed comparison of relative changes in area. The abundance for each focal invertebrate taxon in each SLR scenario at each site was calculated by multiplying the current density of each taxon in each bin by the total area in the respective bin as calculated from the SDEM, and then summing across bins.

Change in intertidal habitat area and zonation with SLR

We estimated habitat area within current and future intertidal elevation ranges at our study sites and found that SLR will substantially reduce total intertidal habitat area (m²) (Fig. 3). Following a SLR trajectory consistent with the observed trend in San Diego, CA (approximately 20.0 cm by 2100), we estimate that total intertidal habitat area loss will be on average 29.88% (±3.78, SE) across study sites. Under the IPCC upper-end global projection of 1.0 m by 2100 (Church et al., 2013), this value will reach 77.72% (±4.65). Under the NRC upper-end projection for California (National Research Council, 2012) of 1.7 m, this value will rise to 85.32% (±2.33). Habitat loss will be greatest for the lower and middle intertidal zones, which currently occupy a broad and gently sloping intertidal shelf that will rapidly become subtidal as sea-levels rise (Fig. 3). Under scenarios greater than 0.2 m the lower intertidal zone will nearly always experience the greatest proportional habitat area loss, followed by the middle, then upper zones (Fig. 3). As a result, we expect that the proportional contribution of each zone to total intertidal area will shift, with the contribution of the lower intertidal zone diminishing, and that of the middle zone and the upper zone increasing (Fig. 3, Fig. S1).

Change in areal extent of dominant benthic space occupiers with SLR

With little exception, we found that the SLR scenarios examined here will result in reductions in areal extent (m²) of dominant benthic space occupiers. Importantly, our data suggest that species will experience changes in areal extent of different magnitudes, resulting in changes in community structure as the relative abundance of species shifts (Fig. 4). The most pronounced changes will likely occur in the first 0.5 m of SLR for all species. At 0.5 m of SLR, we estimate a mean decrease in total areal extent across all species and study sites of 56.95% (±2.40). Changes in areal extent will be most pronounced for those taxa that primarily occupy lower and middle intertidal habitats, such as articulated coralline algae, some brown algae, red foliose algae, turf algae, and surfgrass (see Table S2 for species identifications and classifications). For example, we estimate that the total areal extent of articulated coralline algae at our study sites will decrease by an average of 83.74% (±4.72 SE; range 70.10–91.77) under the IPCC upper-end projection. Areal extent is expected to change less dramatically for taxa occupying primarily upper intertidal habitat, such as mussels and barnacles. For example, we estimate that the total areal extent of barnacles (Balanus/Chthamalus spp.) will decrease by an average of 52.99% (±12.01; range 20.38–77.84) across our study sites under the IPCC upper-end projection.

Change in invertebrate abundance with SLR

Our results suggest that there will be nearly ubiquitous declines in overall numerical abundance of sessile and mobile invertebrates at our study sites with future SLR (Fig. 5). Lower and middle intertidal taxa will likely exhibit greater population declines than upper intertidal taxa, and the largest declines will be observed in the first 0.5–1.0 meter of SLR. For example, our results suggest that the abundance of green anemones (Anthopleura sola/elegantissima) will decrease by an average of 64.37% (±8.66) and 76.20% (±15.60) across our study sites at 0.5 and 1.0 m of SLR, respectively. In contrast, we estimate smaller declines in the abundance of upper intertidal periwinkles (Littorina spp.) with 26.22% (±14.9) and 51.19% (±0.17) and upper intertidal gooseneck barnacles (Pollicipes polymerus (Sowerby, 1883)) with 29.10% (±5.60) and 47.14% (±7.54) across study sites at 0.5 and 1.0 m of SLR, respectively. Overall our results suggest mean decreases in total abundance of all invertebrates at our study sites of 55.82% (±4.25) and 66.92% (±4.09) under the IPCC and NRC projections, respectively. These results, however, do not take into account the possibility of these species becoming increasingly subtidal.