Influences of Wave Climate and Sea Level on Shoreline Erosion Rates in the Maryland Chesapeake Bay

Abstract

We investigated spatial correlations between wave forcing, sea level fluctuations, and shoreline erosion in the Maryland Chesapeake Bay (CB), in an attempt to identify the most important relationships and their spatial patterns. We implemented the Simulating WAves Nearshore (SWAN) model and a parametric wave model from the USEPA Chesapeake Bay Program (CBP) to simulate wave climate in CB from 1985 to 2005. Calibrated sea level simulations from the CBP hydrodynamic model over the same time period were also acquired. The separate and joint statistics of waves and sea level were investigated for the entire CB. Spatial patterns of sea level during the high wave events most important for erosion were dominated by local north-south winds in the upper Bay and by remote coastal forcing in the lower Bay. We combined wave and sea level data sets with estimates of historical shoreline erosion rates and shoreline characteristics compiled by the State of Maryland at two different spatial resolutions to explore the factors affecting erosion. The results show that wave power is the most significant influence on erosion in the Maryland CB, but that many other local factors are also implicated. Marshy shorelines show a more homogeneous, approximately linear relationship between wave power and erosion rates, whereas bank shorelines are more complex. Marshy shorelines appear to erode faster than bank shorelines, for the same wave power and bank height. A new expression for the rate of shoreline erosion is proposed, building on previous work. The proposed new relationship expresses the mass rate of shoreline erosion as a locally linear function of the difference between applied wave power and a threshold wave power, multiplied by a structure function that depends on the ratio of water depth to bank height.

Introduction

Estuaries and coastal embayments of the Mid-Atlantic region have been significantly impacted by coastal erosion, loss of submerged aquatic vegetation (SAV) (Orth and Moore 1984; Stevenson et al. 1993), loss of marsh, and increasing shoreline hardening (Berman et al. 2004). Erosion leads to property loss, increased levels of suspended sediments and turbidity, and increased nutrient loads, degrading nearshore water quality. Shoreline hardening can be an effective preventive measure for shore erosion (Dean et al. 1982). However, eroding shorelines can also provide ecosystem services, such as beach habitat and beetle habitat, and a source of sediment for SAV beds and marshes (National Academies of Sciences 2007). These sources are cut off when the shoreline is hardened. Furthermore, shoreline hardening permanently modifies the delicate and ecologically important land-water interface.

Increasing rates of sea level rise accompanying climate change will likely lead to increasing rates of shoreline erosion. As sea level rises, unprotected shorelines tend to recede to maintain an equilibrium beach profile (Bruun 1962; Schwartz 1967; Dean 1991), with wave attack as the implicit forcing. Shoreline recession due to translation of an equilibrium profile is generally several orders of magnitude greater than the associated increase in sea level. Equilibrium beach profile laws were developed based on studies of sandy shorelines, and they have served as a basis for several erosion models for sandy beaches, such as the Storm-Induced Beach Erosion (SBEACH) and Generalized Model for Simulating Shoreline Change (GENESIS) models.

The relative role of wave attack is less clear for other shoreline types. For example, Benumof et al. (2000) found that material strength appears to largely determine sea-cliff retreat rates, with wave power as a secondary effect, but experimental results for soft cliffs found that erosion was correlated with oblique wave power (Damgaard and Dong 2004). Wave power was also found to correlate with erosion rates for glacial till bluff in Lake Erie (Kamphuis 1987), marsh shorelines in Rehoboth Bay (Schwimmer 2001), and uniform cohesive bluffs in Lake Ontario (Amin and Davidson-Arnott 1997). A recent study in Hog Island Bay (in Virginia) reinforces the important role of wave power in driving erosion along marsh edges (Mcloughlin et al. 2015). Thus, wave power is likely one of the most important factors for predicting erosion rates of non-beach shores, but the global generality of locally derived relationships is not clear and other factors may also be important.

Past efforts in Chesapeake Bay (CB) to examine relationships between shoreline erosion rates and forcing factors such as wave energy, tidal height, and bank characteristics have met with limited success, though waves are clearly an important factor. Some studies have considered the ratio of silt to sand, bluff height, and cohesive soil strength as predictive variables for erosion rates (Dalrymple 1986; Wilcock et al. 1998), with less attention to the effect of wave attack. Other studies have identified waves as the most important factor out of several possibilities (Wang et al. 1982; Spoeri 1985; Skunda 2000; Perry 2008). Skunda (2000) modeled shoreline instability along the western shore of CB in Virginia (mostly sandy shorelines) and found wave power to be the dominant factor for erosion. Spoeri et al. (1985) and Wang et al. (1982) used data from 107 reaches (2–5 km in length) in the Maryland CB to analyze relationships between erosion rates and wave energy, sediment type, tidal range, rainfall, and 100-year storm surge. Using traditional regression and discriminant analysis, Spoeri et al. (1985) were not able to provide an adequate erosion rate prediction model, but they concluded that wave energy still seems to be primarily responsible for the changes in shoreline erosion. Perry (2008) applied discriminant analysis (CART) and linear models to shoreline data sets in the Maryland CB and found that fetch seems to be the most important factor, again highlighting the role of waves, while geography is the second most important predictor, complicating efforts to develop a single model for Bay-wide predictions. Thus, the search for a better predictor of shoreline erosion rates in the CB remains an important, but elusive goal.

In this study, two modern wind-wave models were implemented in CB for long-term wave climate simulations from 1985 to 2005, and wave variability was compared to separately predicted sea level variability over the same time period. Spatially varying wave and sea level climates were then compared to spatial distributions of historical shoreline erosion rates, shoreline structure information, bank/marsh ratio, and mean bank height in 207 reaches for the Maryland part of CB. The erosion part of the study focused on the Maryland portion of CB because an equivalent detailed shoreline inventory (combining erosion rates with shoreline characteristics) was not available for the Virginia portion of CB.

The datasets assembled for this study are the most comprehensive to date to be applied to the study of shoreline erosion in the Maryland CB. They cover longer shorelines (including the major tributaries), and the climatological forcing is from more advanced numerical models that provide more comprehensive predictions than in previous studies. The ultimate goals of the study were to gain a better understanding of the contributions of a variety of controlling variables to shoreline erosion rates, to establish relatively simple semi-empirical or statistical relationships between erosion rates and these controlling variables, and to formulate an erosion expression incorporating these relationships to provide a framework for future studies and shoreline management efforts.

Study Site

Chesapeake Bay (CB; Fig. 1) is the largest estuary in the USA. It has a relatively shallow average depth of 8.5 m, length of 315 km from the mouth of the Susquehanna River to the Atlantic Ocean, and width that ranges from 5.6 to 56 km (Langland and Cronin 2003). The dendritic shorelines of CB stretch more than 18,800 km (http://www.chesapeakebay.net/discover/bay101/facts). The axis of the Bay runs northeast-southwest in the upper Bay, north-south in the mid-Bay, and northwest-southeast near the Bay mouth. Mean tidal range is lowest in the mid-Bay (∼0.3 m) and increases towards both ends to slightly less than 1 m (Zhong and Li 2006). CB is a partially mixed estuary and has a classic two-layer estuarine circulation, with fresher upper-layer water flowing seaward and saltier lower-layer water flowing landward over subtidal time scales (Pritchard 1956). Ocean swells enter CB between the Virginia Capes, but they are largely dissipated before they reach the mouth of the Potomac River in mid-Bay (Boon et al. 1996). The approximately north-south orientation of CB provides long fetches for generation of wind-forced surface gravity waves, which dominate the upper- and mid-Bay regions (Boon et al. 1996; Lin et al. 2002). For these waves, typical periods (T) are 3–4 s and typical significant wave heights (Hsig) are less than 1 m (Lin et al. 2002). Waves are likely the dominant forcing for sediment transport in shallow nearshore CB waters (Sanford 1994; Boon et al. 1996) and for shoreline erosion (U.S. Army Corps of Engineers 1990).

In the Maryland portion of CB, shorelines consist of banks with heights ranging from 1 m to over 30 m (at Calvert Cliffs) and marshes that are mostly found along the lower eastern shore (Somerset, Wicomico, and Dorchester Counties). Some regions, points, and islands are experiencing severe erosion (>2.4 m/year) on the western shore (Pt. Lookout to St. Jerome, Holland Pt., and Thomas Pt.) and eastern shore (Kent Island, Lowes Pt. to Knapps Narrows, Mills Pt. to Hills Pt., James Island, Oyster Cove to Punch Island Creek, and Barren Island) (Wang et al. 1982). Land loss from shoreline erosion totaled 1.9 × 108 m2 between 1850 and 1950 (Slaughter 1967). Shoreline erosion can lead to nutrient pollution, ecosystem degradation, and huge economic loss (Leatherman et al. 1995). Erosion of unprotected shorelines has been intensified by sea level rise, land subsidence, and increasing rates of shoreline development (Halka et al. 2005).

Methods

This study implemented wave models to produce climatological wind-wave estimates over a 21-year period, combined these with previously calculated model estimates of sea level over the same time period, compiled historical shoreline erosion rates in the Maryland CB, and compiled co-located data on shoreline characteristics. The erosion rate data were compared to the wave and sea level climatologies and to shoreline characteristic data at two spatial resolutions to maximize the value of the analysis. It was anticipated that using 21-year climatologies for waves and sea level might provide a better temporal match to the historical erosion rate estimates, which are derived from shoreline differencing over comparable time scales.

Wave and Sea Level Models

Simulating WAves Nearshore (SWAN) is a third-generation wave model that is based on Eulerian formulation of the discrete spectral balance of wave action density (Booij et al. 1999; Ris et al. 1999). This model simulates random, short-crested waves in coastal regions with arbitrary bathymetry, winds, and current fields. SWAN is driven by boundary forcing and local winds. It accounts for triad wave-wave interaction, shoaling, refraction, white capping, bottom friction, depth-induced breaking, dissipation, and diffraction. Version 40.91 of SWAN was used in this study. The model was configured to run in third-generation mode with 2D propagation and refraction on a boundary-fitted curvilinear Cartesian grid. The model grid and bathymetry were identical to that used by the US Army Corps of Engineers (USACE) for the US EPA Chesapeake Bay Program (CBP) hydrodynamic and water quality models (Fig. 1; Cerco and Noel 2004). Grid resolution was approximately 1 km in the axial direction and 0.4 km in the lateral direction, with a total grid size of 282 × 178 cells. Wave period cutoff limits were 0.001 and 25 s, the peak period was used as a characteristic wave period, the wave growth term from Cavaleri and Malanotte-Rizzoli (1981) and the surface drag coefficient from Wu (1985) were used, and the default JONSWAP coefficient of 0.067 was adopted for bottom friction (Hasselmann et al. 1973). Activated physical processes included white capping, non-linear quadruplet wave interactions, and triad wave-wave interactions. The calculation time step was 10 min, a significant wave height of 0 m and peak period of 0.1 s were used as ocean boundary conditions, and zero wave height and wave period everywhere were used for initial conditions. Ignoring incoming ocean swell at the Bay mouth was necessary because neither ocean data nor ocean wave model results were available for 1985–2005.

Young and Verhagen (1996) developed a semi-empirical method for calculating fetch-limited wave growth that calculates wave heights and periods from water depths, wind inputs, and fetch. They developed fetch-limited and depth-limited shallow-water wave relationships, with empirical parameters determined based on extensive wind-wave measurements at Lake George, Australia. The Engineer Research and Development Center (ERDC) of the USACE developed a parametric wind-wave model for CB based on Young and Verhagen’s formulation (Harris et al. 2012). The most complex calculation in this model is determining smoothed effective overwater fetch at each grid point for each time step, which is accomplished by allowing for cosine-weighted directional averaging to avoid sudden changes in fetch due to slight changes in wind direction. The model was implemented for this study on the same grid, with the same bathymetry and the same winds as the SWAN model. This wave model is called the CBP wave model herein.

Hourly outputs of sea level for the entire CB from 1985 to 2005 were available from the Waterways Experiment Station (CH3D-WES) model (Johnson et al. 1993), referred to as the CBP hydrodynamic model in this study. Sea level is among the best calibrated outputs of the CBP hydrodynamic model (Johnson et al. 1993; Cerco and Noel 2004). Sea level variability at the ocean boundary was directly interpolated from observations for the CBP hydrodynamic model. The wave models accounted for sea level variability by adjusting local water depth using interpolated hourly output from the CBP hydrodynamic model.

The winds used to force both wave models were identical to the winds used to drive the CBP hydrodynamic model. Hourly wind velocity vectors were interpolated onto model grid points from five wind stations: Thomas Point, Patuxent Air Base, Norfolk International Airport, Washington National Airport, and Richmond International Airport (Fig. 1). The resultant wind direction patterns (which varied over the domain) were not adjusted, but wind speeds from the four over-land stations were adjusted by multiplying by 1.5 prior to interpolation to better match data at the over-water Thomas Point station (Johnson et al. 1993; Harris et al. 2012). This adjustment approximately compensated for greater attenuation of overland wind speed in the terrestrial atmospheric boundary layer, as compared to the marine atmospheric boundary layer (Goodrich et al. 1987; Xu et al. 2002).

Predictions of both wave models were compared to the same wave observations as used for validation by Lin et al. (2002). Near-bottom pressure-based wave gauges were deployed from October 10–23, 1995, at Calvert Cliffs and from October 26 to November 9, 1995, at Poplar Island, both in mid-CB. The comparisons (Fig. 2) show that both SWAN and the CBP wave model work reasonably well. Predictions at Poplar Island show a better match than at Calvert Cliffs in terms of both significant wave height (Hsig) and peak period (Tp). At Calvert Cliffs, predictions of Hsig match the data better than Tp predictions from either model. As shown qualitatively in Fig. 2, SWAN did slightly better than the CBP wave model for both Hsig and Tp. In addition, SWAN provides direct estimates of many other wave-associated variables that are potentially important factors for shoreline erosion. Thus, SWAN predictions were used for all model wave data except for fetch, which is only calculated by the CBP wave model.

Comparisons of time series between model outputs and observations at a Poplar Island and b Calvert Cliffs. For each location, panels show (upper to lower) significant wave height (Hsig), peak period, and interpolated winds used to drive the two wave models and the CBP hydrodynamic model

Shoreline Erosion Rate Estimates and Shoreline Characteristics

Shoreline erosion rate estimates were derived for the Maryland portion of CB, including its major tributaries and islands. Some interior ponds, creeks, and heads of tributaries were not included due to limited spatial resolution. The Maryland Geological Survey (MGS) assembled data on shoreline erosion, shoreline structure percentage, bank percentage, and mean bank height into one dataset at “reach” resolution (Fig. 3; Hennessee et al. 2006). Reaches were defined from one point of land to another, then further subdivided if rates of shoreline change varied widely within the reach. The mouths of tributaries, county boundaries, and marked changes in shoreline orientations all influenced reach scopes (Hennessee et al. 2006). MGS demarcated the Maryland CB shoreline into 207 reaches, which were divided almost equally between the eastern shore (100 reaches) and the western shore (107 reaches). Reaches ranged from 1.9 to 87 km in length.

MGS worked with the US Geological Survey (USGS) and Towson University’s Center for Geographic Information Sciences (CGIS) to estimate historical erosion rates for the coastal and estuarine shorelines in Maryland. They used digitized shorelines dating from 1841 to 1995 as inputs into the Digital Shoreline Analysis System (DSAS; Danforth and Thieler 1992), which produced nearly 250,000 shore-normal transects with associated rates of change, including the Atlantic coast, the coastal bays, and the CB and its tributaries (Hennessee et al. 2002; 2003a,b). Transects were spaced approximately 20 m apart. This database is available as a product called Coastal Atlas on the Maryland Department of Natural Resources (DNR) website. In order to calculate a representative shoreline erosion rate for each reach, MGS averaged the most recent erosion rate estimates for all transects within that reach, excluding transects that intersected protected shorelines. Thus, the reach-averaged erosion rates only include unprotected shoreline. In the Coastal Atlas convention, negative rates of change represent shoreline erosion, while positive rates of change indicate accretion. This sign convention is adopted here for consistency, but it must be kept in mind when comparing our results to those of other studies.

The Virginia Institute of Marine Sciences (VIMS) identified land use, bank condition, and man-made structures along the tidal shorelines of navigable waterways in Maryland from 2003 to 2005 (http://ccrm.vims.edu/gis_data_maps/shoreline_inventories/). Observations were made from a slow-moving boat traveling parallel to the shoreline and were organized into a geographically referenced set of shoreline condition data. In addition, the percentage of protected shoreline length in each reach was used as calculated by MGS (Hennessee et al. 2006). MGS used 7.5-min USGS topographic quadrangles to identify marshes and topographic contours along the shorelines. These two features were used to estimate the percentage of bank shorelines vs. marsh shorelines for each reach and the average bank height for each reach (Hennessee et al. 2006).

In order to evaluate relationships between waves, sea levels, and shoreline erosion at finer resolution, erosion rates were also estimated at grid scales (approximately 1 km). Each shoreline transect was manually assigned to the most appropriate model grid, and the average shoreline erosion rate for each grid was calculated. Distinguishing between protected and unprotected transects was not feasible for these calculations. Thus, for reaches with significant fractions of protected shoreline, these grid-averaged erosion rates were biased by the inclusion of transects that had been protected for unknown periods of time. Since marsh grids were previously determined to be mostly unprotected, only marshy shorelines were selected for quantitative analysis at the model grid resolution in this study. Grid cell resolution did not adequately resolve sections of the Bay side (high erosion rates) or the sheltered side (low erosion rates) of a few islands. This caused erosion rates at grid-cell resolution to be underestimated along the Bay sides of Taylors Island, Hoopers Island, Barren Island, and Bloodsworth Island, as well as the upper part of the Bay side of Smith Island (including Martin National Wildlife Refuge), so these cells were not included in the grid-scale analysis.

As a result of the calculations and analyses described above, the assembled data included estimates of wave and sea level forcing factors for the entire CB shoreline every hour for 21 years at grid cell resolution; estimates of average unprotected shoreline erosion rates in Maryland at reach scale; estimates of average bank height, bank fraction, and protected shoreline fraction in Maryland at reach scale; and estimates of marshy shoreline erosion rates in Maryland at grid scale. There were 2217 shoreline grid cells around the entire CB, 1316 shoreline grid cells in Maryland, 207 shoreline reaches in Maryland, and 163 useable marshy shoreline grid cells in Maryland.

Appropriately matching the temporal and spatial resolutions of the various data sets required several additional steps. First, various metrics for waves and sea level were calculated and averaged over the 21-year simulation period at grid cell spatial resolution. This long temporal averaging was done in order to match the decadal averaging of historical erosion rate estimates. In Maryland, grid cells were manually assigned to the closest corresponding reach (Fig. 3) and the associated wave and sea level metrics were then spatially averaged for comparison to reach-averaged erosion rates and shoreline characteristics.

In order to derive the onshore component of wave power, the angle between the shore-normal direction and the incoming wave energy flux (α) was calculated as α = β − θ + 90°, where β is the direction of incoming wave energy propagation calculated from the directional components of the wave energy flux and θis the shoreline orientation (Fig. 4). Nearshore wave direction and onshore wave power were calculated at the center of each shore-adjacent grid cell, at a typical distance offshore of ∼200 m and a typical water depth of ∼1.5 m. SWAN accounts for wave breaking, but typical CB waves are not yet breaking under these conditions. The angle calculations follow the geometric convention that east is 0°, with angle values increasing counterclockwise. cosα > 0 represents offshore-directed waves, cosα < 0 represents onshore directed waves, and cosα = 0 represents alongshore-directed waves. The average onshore wave power was calculated as

Schematic diagram illustrating the calculation of α, the angle between the incoming waves, and the normal to the local coastline

All offshore-directed wave power values were set to 0 in the numerator, since offshore wave power neither contributed to nor subtracted from onshore wave power. Transp_onshore is always negative because of the sign conventions used here, but its absolute value was used for analysis in this study to ensure that correlations between erosion rate and important forcing terms all had the same signs for similar effects.

Time series of fetch were acquired from the CBP wave model; the average wind-weighted fetch at reach i was calculated as

Tidal Range is not the conventional definition of tidal range, but rather the standard deviation of sea level, which is proportional to conventional tidal range—it is referred to as Tidal Range in this study for simplicity. All depths (from the model bathymetry) are relative to mean sea level (MSL) in 1983. Modeled bathymetric steepness was calculated as

where the denominator is the distance offshore at the center of each shoreline model grid. Bank_Height is the elevation of the bank top above MSL. MedianWater60Hsig was calculated as the median value of sea level for the highest 40% of waves at each shoreline grid cell; it is a measure of typical local sea level during high wave events.

For the dataset at reach resolution, shoreline types were identified as “marsh,” “bank,” or “mixed” type based on their bank fraction. A reach was defined as marsh if the bank fraction was ≤0.1 and as bank if the bank fraction was ≥0.9. Because marshes and banks eroded quite differently (shown below), this study focused primarily on reaches that fell into bank (117 reaches) or marsh (27 reaches) categories. Bank shorelines were further subdivided into subregions of main stem, tributary, eastern shore, and western shore. Marsh data at the reach scale had too few points for subdivision, but marsh data at the grid scale were further subdivided into main stem, tributary, eastern shore, and western shore for analysis. In the end, there were 11 groups of data considered separately on the basis of shoreline type or location. Geographic subdivision was largely motivated by conclusions about the importance of geographical location for grouping shoreline erosion rate estimates reported by Perry (2008).

Linear correlation analyses, including calculations of Pearson correlation coefficients and multiple linear regression (MLR), were performed on the full data sets. Non-parametric and non-linear general additive models (GAM; Hastie and Tibshirani 1990) and neural network (NN; Beale et al. 2014) analyses were used selectively to characterize statistical relationships when data sets were too complex or too non-linear for simple linear techniques. Simple calculation of Pearson correlation coefficients among all potentially influential variables turned out to yield the most information for present purposes, with little additional value added by multiple linear regression (MLR), general additive models (GAM), or neural network (NN) analyses.

Erosion rate outliers were excluded using quantile ranges. Data outside the range of Q1 − 3*IQ and Q3 + 3*IQ were identified as outliers, where Q1 is the 25th percentile, Q3 is the 75th percentile, and IQ equals Q3 − Q1. For bank data, the upper bound was relaxed to −0.96 instead of Q3 + 3*IQ = −0.79 in order to preserve more reasonable data for analysis. At reach scales, one bank data point (the Bay side of Hoopers Island and Barren Island) and one marsh data point (Taylors Island) were excluded as outliers. For marsh data at grid cell scale, three outliers were excluded on the Bay side of Smith Island.

Results

Wind Wave and Sea Level Climate, 1985–2005

This section presents the distributions of significant wave height, sea level, and their joint probability along the entire shoreline of CB, including the major tributaries. The spatial distribution of typical sea level during locally high wave events is also shown and discussed. The data for wave height in the lower CB do not include ocean swell, so it should be taken as lower bounds for actual wave heights; ocean swell would be highest at the Bay mouth and decrease northward, were it to be included. The lower Bay data are included here for completeness and because they help to illustrate overall trends and interactions between wind waves and sea level, but they were not used for the erosion analysis.

The spatial distributions of 21-year-averaged Hsig, sea level, and tidal range are shown in Fig. 5. Average Hsig varied between 0 and 0.3 m around the shorelines of CB (Fig. 5a). The smallest waves occurred in tributaries, with Hsig decreasing from tributary mouths towards their heads as expected due to increasingly limited fetch. Larger waves were present along the exposed shorelines of main stem CB, with highest wave heights near the mouth of the Bay on both its western shore (Cape Henry) and eastern shore (Cape Charles). Twenty-one-year averaged Sea Level and Tidal Range (Fig. 5b) showed similar patterns throughout much of the Bay, with a few notable differences. Average sea level was close to 0 at the mouth of CB (recall that sea level is relative to 1983 MSL at the Bay mouth). It increased rapidly with distance into the Bay and then remained relatively flat at approximately 0.075 m through the mid-Bay, increasing again towards the head of the Bay. There were large, distinct peaks in average sea level at the heads of the major tributaries corresponding to large inputs of freshwater. Average tidal range decreased rapidly from the mouth of the Bay and remained relatively flat at approximately 0.25–0.3 m through the mid-Bay, increasing again towards the head of the Bay. Tidal range also exhibited distinct peaks at the heads of the major tributaries, with a slightly different pattern than that of average sea level.

Wave and sea level climate along the unwrapped Chesapeake Bay shoreline. Climate values represent averages over 21 years. The shoreline is unwrapped by starting at Cape Henry and tracing all model shoreline cells, including tributaries, around the circumference of the Bay to Cape Charles. Additional points at the end represent islands

Histograms of Hsig and Sea Level are presented in Fig. 6. These histograms represent the distributions of hourly estimates of Hsig and Sea Level over 21 years and for all shoreline grids in CB. Sea level records at each grid location were detrended over time in order to focus attention on variability over hourly to annual time scales, rather than long-term or large-scale trends. The histogram of Hsig in Fig. 6a shows that 91% of all wave heights were within the lowest bin of 0–0.27 m and that the number of data points decreased exponentially with increasing Hsig. The probability of observing a particular Hsig (P(Hsig)) can be fit to an exponential probability density function (PDF) with R2 > 0.99:

Histograms of aHsig and b detrended Sea Level, from hourly output of the SWAN and CBP hydrodynamic models, respectively, for the entire CB shoreline from 1985 to 2005. Red lines represent best fit probability density distributions given in the text (color figure online)

which is also shown in Fig. 6a. The histogram of detrended Sea Level (Fig. 6b) shows that the variability of sea level was approximately normally distributed. The PDF of sea level is very well described by a normal distribution with R2 > 0.99:

where w stands for detrended sea level; Eq. 5 is also shown in Fig. 6b. Note that the mean of this distribution has a very small negative offset because the data are slightly positively skewed.

Figure 7 shows the joint probability distribution of Hsig and Sea Level. For low wave heights, which are by far the most common (Fig. 6a), sea level was normally distributed with a mean of approximately 0, in agreement with Fig. 6b. However, as Hsig increased the most probable value of Sea Level increased, reaching a value of approximately 1 m above MSL for very high but very infrequent waves.

Joint probability distribution of detrended Sea Level and Hsig for all CB shoreline grid cells at every time step from 1985 to 2005

The sign and magnitude of correlations between high sea level and high waves were spatially non-uniform, however. Figure 8 illustrates this by showing the spatial distribution of MedianWater60Hsig, the median value of sea level for the highest 40% of waves at each shoreline grid cell. To produce Fig. 8, a value of 1 and the color red were assigned to positive values of MedianWater60Hsig while a value of −1 and the color blue were assigned to negative values of MedianWater60Hsig. If counts of positive and negative sea level were approximately the same (within 1% of a 50/50 split), a value of 0 and the color green were assigned. Figure 8 shows that positive MedianWater60Hsig occurred primarily on the western shore in the main stem Bay, the southern sides of tributaries in the lower Bay, and the northern side of tributaries in the upper Bay. Negative MedianWater60Hsig occurred primarily on the eastern shore of the main stem Bay, the northern sides of tributaries in the lower Bay, and the southern sides of tributaries in the upper Bay. These distributions are explained by the different responses of sea level to wind in different parts of the Bay (e.g., Garvine 1985). In the upper Bay, sea level primarily responds by setting up or setting down in response to local south or north winds, respectively. In the lower Bay, sea level primarily responds to set up or set down in the adjacent coastal ocean due to northeast winds or southwest winds, respectively. In both cases, the most frequent combination of large waves and sea level occurs in response to winds that have the longest over-water fetches. For example, in the upper Bay, winds out of the south-southeast increase sea level and produce high waves along the western shore, while winds out of the north-northwest decrease sea level and produce high waves along the eastern shore. In the lower Bay, winds out of the northeast increase sea level and produce high waves along the western shore, while winds out of the southwest decrease sea level and produce high waves along the eastern shore.

Median detrended sea levels for highest 40% of local waves (MedianWater60Hsig) reflect different regional influences of remote or local wind forcing. Red dots indicate grid cells for which sea level tends to be high when waves are high. Blue dots indicate grid cells for which sea level tends to be low when waves are high (color figure online)

Spatial distributions of Maryland shoreline erosion rates at both grid cell resolution and reach resolution are shown in Fig. 9 for comparison to the distributions of onshore wave power, bank height, and bank fraction. These three factors are all expected to significantly influence erosion rates. Erosion rate estimates at grid cell resolution included both protected and unprotected shorelines; including protected shoreline might bias the averages downward (in a positive direction according to the present sign convention). Erosion rate estimates at reach resolution included unprotected shorelines only, so that they represent more of an unbiased natural response. The two estimates generally agreed well except for two high peaks of accretion at grid cell resolution. These peaks, in Baltimore Harbor and at Hart-Miller Island, corresponded to man-made shoreline extensions.

Erosion rates varied widely along the shorelines of the Maryland CB, from 0 to approximately 3 m/year (Fig. 9). Using grid cell 1273 (the mouth of the Susquehanna River) as the dividing line between eastern and western shores, eastern shorelines generally eroded faster than western shorelines. The Bay sides of islands experienced the most severe erosion. The most rapidly eroding islands included Taylors Island, Hoopers Island, Barren Island, Bloodsworth Island, and Smith Island on the eastern shore, and Pooles Island on the upper western shore. Most of these islands have marshy shorelines, except for the Bay sides of Hoopers Island and Barren Island. Onshore wave power was minimal in tributaries, especially near their heads, and was highest along open CB shorelines between the tributaries. Most high erosion rates corresponded to high onshore wave power values, but high onshore wave power alone did not necessarily indicate a high erosion rate. Average bank height, at reach resolution, was generally less than 5 m, except for the 30-m-high banks near Calvert Cliffs on the western shore. Banks constituted a much larger percentage of shoreline than marshes in the Maryland CB. Some marshy shorelines occurred along the upper western shore, but the majority of marshy shorelines occurred along the lower eastern shore and islands.

Though not immediately apparent in Fig. 9, there was a significant relationship between bank percentage and shoreline erosion rate. Mean erosion rates decreased as bank percentages increased (r2 = 0.44; P << 0.05; Fig. 10a) when bank percentages were grouped into five categories. The means of different categories were significantly different from each other as well (ANOVA; 99.9%). The differences of mean erosion rates between 0–10% and 70–90%/90–100% bank, as well as the differences between 90–100% and 10–40%/40–70%, were significant with better than 95% confidence based on Tukey’s test. Banks experienced both erosion and accretion, while marshy shorelines only experienced erosion. Bank height was not significantly correlated with shoreline erosion rate, but mean erosion rates did decrease as bank height increased for bank heights between 0 and 6 m (Fig. 10b). ANOVA analysis indicated that there was 92% confidence that the mean among the three categories in the 0–6-m bank height range was different. Furthermore, Tukey’s test identified a significant difference between the mean erosion rates of bank heights 0–1.5 and 3–6 m, but differences between mean erosion rates among other height combinations was not significant. Two data points representing bank heights higher than 10 m (near Calvert Cliffs on the western shore) showed moderate erosion; these two points were treated as outliers in both the ANOVA and Tukey’s test because there were far fewer data points than in the other categories.

Erosion rate vs. bank percentage at the reach scale in the Maryland CB (left) and erosion rate vs. bank height at the reach scale in the Maryland CB (right)

The different erosion behaviors of marshes and banks are illustrated in Fig. 11, which compares marsh erosion to only low-elevation bank (0–1.5 m) erosion at reach scale, both plotted against onshore wave power. Erosion rates of marshes increased approximately linearly with increasing onshore wave power, beginning at very low wave power. Erosion rates of low banks remained relatively low and constant (except for two sites) until wave power exceeded approximately 0.015 kW/m, when erosion rates began to increase. At low onshore wave power, both low banks and marshes appeared to experience a low background erosion rate of about −0.01 to −0.02 m/year.

Correlations between the listed variables and erosion for all 203 shoreline reaches (with four outliers excluded), aggregating all shoreline types and locations together, were not significant and are not shown in Table 2. This is likely because there were so many different factors affecting erosion of this heterogeneous data set. For example, erosion rates of marshes and low-elevation banks behaved differently with respect to increasing onshore wave power (Fig. 11), indicating that the two shoreline types may need to be treated differently. The plots of bank percentage and bank height against erosion rate in Fig. 10 also show large scatter, presumably due to the multiplicity of influential factors, which led to no significant Pearson correlation between these two factors and erosion in Table 2.

Higher correlations were obtained by breaking the shoreline data down into more geologically and/or geographically related groups. Based on the total number of significant correlations for the different shoreline groups presented in Table 2, the most important distinctions were between marshy and bank-dominated shorelines, and between tributary and main stem bank groups. The main stem bank group generally had less significant correlations simply because there were so few data points. The eastern and western bank groups had sufficient data points but did not have many significant correlations, indicating that this particular grouping did not add much value to the analysis. Stronger correlations were generally obtained for marsh erosion than for bank erosion. Recall that negative correlations in Table 2 indicate a positive relationship; i.e., an increase in the factor of interest led to an increase in erosion rate. This is simply because erosion is defined as a negative change in shoreline position in our analysis.

The most valuable aspects of the significant correlations shown in Table 2 are the insights into erosion dynamics that they offer. The strongest correlations were between marsh erosion and wave characteristics such as significant wave height, wave period, and wave length, all of which were related strongly to each other. In the marsh data, non-directional wave power (transpall) had the highest correlation, but the onshore wave power correlation (transp_onshore) was almost as high. Only four variables were significantly correlated in four or more shoreline categories: frequency spread (FSPR), near-bottom velocity squared (UBOTsq, proportional to wave-induced bottom stress), onshore wave power (transp_onshore), and Tidal Range. The significant positive correlations between FSPR and erosion rate across several shoreline categories indicates an important aspect of waves that favored erosion. Waves with a smaller FSPR (favoring erosion) generally represent a more well-developed narrow-banded wind sea resulting from strong winds blowing across a long fetch, whereas waves with a larger FSPR represent choppy, relatively undeveloped wind seas and lower erosion rates. The positive correlations between tidal range and erosion are an artifact of the spatial distribution of tidal range in the Maryland CB. Tidal range is higher towards the heads of tributaries (Fig. 5), which is where erosion rates are generally lower because of restricted fetches resulting in lower waves (Fig. 5), hence lower erosion rates (Fig. 9). Finally, correlations with UBOTsq most likely represent the influence of long-period, high waves that scour the bottom, removing eroded material from the shoreface. Waves that increase UBOTsq also favor low values of FSPR and high transp_onshore.

Other factors that were significantly correlated with erosion rates were fetch and wind-weighted fetch, which were essentially equal in their influence; this finding is consistent with the results reported by Perry (2008). MedianWater60Hsig had the most significant correlation with erosion along the western Bay, and the sign of the correlation indicated that high sea level accompanying high waves resulted in lower erosion rates. Interestingly, this also is consistent with Perry’s (2008) indication that waves attacking the base of the bank were likely most influential at producing erosion. Finally, in the Bank Stem data group the steepness of nearshore bathymetry (Bath_steepness) was positively correlated to erosion rate; in other words, a steeper nearshore bathymetry accompanied lower erosion rates. This correlation may actually represent a reverse causality—banks that erode slowly may not have the extended shallow nearshore shelf of banks that erode rapidly, and/or may consist of coarser sands that tend to have steeper offshore depth profiles (Dean 1991).

Discussion

The results presented above, while they do not identify any one factor that completely explains the spatial distribution of erosion rates in the Maryland CB, do contribute to a consistent dynamical description that may be summarized as follows. Wave forcing is the dominant forcing factor for shoreline erosion in the Maryland CB. Well-developed wind seas that are formed over long fetches and have large wave heights with long periods, hence high onshore wave power, tend to be the most effective at producing shoreline erosion. Onshore wave power is the most meaningful metric for wave forcing, since it directly represents the physical rate of doing work on the shoreline. Sea levels at or near the base of the shore face during high wave forcing are most effective at causing erosion, since they focus wave power at the location most likely to lead to subsequent shore face collapse. Local shoreline characteristics are also important, however. Marshes erode much faster than banks, other things being equal, which also means that sections of shoreline with higher fractions of marsh erode faster on average. Low banks tend to erode faster than high banks. Finally, there are a number of factors besides the bank height or amount of marsh that were not considered here, but likely contribute significantly to the unexplained variability in the data. Most significant among these are likely to be the sediment texture and cohesion/erosion resistance of the shoreline (Dalrymple 1986; Wilcock et al. 1998; Ravens et al. 2012).

By combining the above summary of empirical shoreline erosion dynamics in the Maryland CB with an understanding of the processes of shoreline erosion from the literature and relationships derived by other investigators, a semi-empirical relationship is proposed below that is consistent with the data presented here and may prove useful. The actual process of shoreline erosion usually proceeds through a sequence of events: waves undercut the bank/marsh base; the bank/marsh collapses; and waves and currents working together transport these materials offshore (Damgaard and Dong 2004; Phillips 1999; Ravens et al. 2012). On sandy beaches, offshore transport alternates seasonally with onshore transport, leading to an annually averaged equilibrium beach profile and no net erosion except as required to maintain the equilibrium profile in the presence of sea level rise (Dean 1991). On muddy or mixed shorelines such as those that dominate the Maryland CB, fine sediments are transported offshore permanently and a true equilibrium shoreline may never occur. In any case, an important part of the shoreline erosion process is transport of sediment mass offshore. Thus, expressing erosion rate in terms of sediment volume or mass loss should be more general than expressing erosion rate in terms of horizontal recession rate, which is the most common definition and the one used here.

Recent studies focusing on marsh edge erosion have made significant progress in this direction. In particular, Marani et al. (2011) used dimensional analysis to derive a relationship of the form

$$ \frac{Ehc}{P}= f\left(\frac{h}{d}\right) $$

(6)

where E is erosion expressed as the linear rate of shoreline recession, h is total bank height, c is a measure of shoreline cohesion, P is the average onshore wave power, f is a function to be determined, and d is the sea level at mean tide. The product Eh is the volumetric erosion rate (VER), which Marani et al. showed to be linearly related to P in a number of marsh erosion data sets, with the product c−1f(h/d) as an empirical proportionality factor.

Figure 12 shows that an expression of this form gives a reasonable fit to the marshy shoreline erosion estimates presented here. In the present study, marsh elevation was not measured directly but was assumed to be constant at 0.5 m. With an assumed constant marsh elevation, there is only a constant multiplicative difference between horizontal shoreline recession rates and volumetric erosion rates. Thus, Eq. 6 becomes

To the extent that the expression inside the square brackets in Eq. 7 is approximately constant for similar marshy shorelines, Eq. 7 predicts an approximately linear relationship between shoreline erosion rate and wave power. Figure 12a shows that an approximately linear relationship between E and P fits our data reasonably, especially at reach resolution. The derived least squares linear fit is given by E (m/year) = −54P (kW/m) − 0.08, with r2 = 0.55. While the data at grid resolution also follow the same general trend as the data at reach resolution, the fit is not significant, possibly because the terms inside the brackets in Eq. 7 are variable even at small scales. The large variability in wave height (hence wave power) from cell to cell shown in Fig. 5a also might have contributed to the lack of significance in the model fit at the smallest scales considered here. McLoughlin et al. (2015) noted that relationships between wave power and marsh erosion at the finest scales in Hog Island Bay, VA, were more variable (and less significant) than relationships averaged over a representative length of shoreline. Scale dependence of shoreline erosion processes is clearly a topic for additional research.

Figure 12a also shows approximate linear relationships from two other marsh erosion studies (Schwimmer 2001; McLoughlin et al. 2015). While the forms of the fits are similar, the proportionality between E and P can be quite different between different studies (up to two orders of magnitude in these cases). The differences between studies may be due to differing cohesion, vegetation, wave climate, marsh elevation, tide range, time/space scales of observation, or differences in methods of calculating wave power and/or erosion rate. For example, many studies estimate wave power by first estimating characteristic wave heights, periods, and directions and then calculating characteristic wave power using a standard formula for monochromatic waves (e.g., Mariotti et al. 2010; Marani et al. 2011). In the present study, we estimated wave power as the average of long-term time series of directional wave energy flux. Differences may arise because one calculation depends on the square of time-averaged wave heights, while the other calculates the time average of instantaneous estimates of squared wave heights.

Recently, Leonardi et al. (2016) have shown that much of the inter-site variability can be removed by normalizing E and P by their respective local average values. Combining data from marsh erosion studies at eight different sites, including the widely varying studies of Schwimmer (2001) and McLoughlin et al. (2015) shown in Fig. 12a, Leonardi et al. (2016) showed that normalizing by local averages collapses all of the observations into a single relationship given by E* = −0.67 P*, where E* and P* are normalized erosion rate and onshore wave power, respectively. Figure 12b shows that normalizing our marsh data in the same way results in a relationship given by E* = −0.88 P* − 0.14; while this is faster erosion than predicted by the Leonardi et al. (2016) relationship, it is within the scatter of their combined data set. The implication is that knowledge of average erosion rate and onshore wave power at a given site will allow prediction of time-varying marsh erosion rates due to variations in wave power, with an empirical proportionality factor given by the ratio of the averages. This is conceptually similar to the empirical proportionality factor in Eq. 7 above, although in the Leonardi et al. (2016) approach there is no specific causality assigned to the value of the proportionality factor.

Equation 6 can be further modified to make it even more general, and to potentially apply to more than just marshy shorelines. The first step is to express erosion rate in terms of sediment mass by multiplying volumetric erosion by the dry bulk density of the bank/marsh sediment, such that the LHS of Eq. 6 now becomes (Ehρdry)c/P. This of course requires an adjustment of the units of c, and an assumption that the mass rate of erosion is the fundamental quantity of interest rather than the volumetric rate of erosion. Variable dry bulk density along the eroding shorelines of CB, in addition to varying bank height, might be responsible for some of the unexplained variability between wave power and linear erosion rate seen in this study. Sediment type and availability for all shoreline reaches in the Maryland CB are unknown, but there are data on sediment characteristics for selected banks and marshes. Analysis of 76 MGS sediment samples from 21 bank sites shows an average composition of 44% sand and gravel, 56% silt and clay, and negligible organic matter. Twenty sediment samples from four marsh sites show 22% sand and gravel, 44% silt and clay, and 34% organic matter (Hennessee et al. 2006). This leads to a factor of 2 difference in dry bulk density between marshes and banks (Table 3), along with possible differences in cohesion. The lower dry bulk density of marsh sediments would cause a linear recession rate for marshes that was twice as fast as erosion for banks of the same height (e.g., Fig. 11), for the same mass erosion rate.

Table 3

Average bulk density of banks and marshes from selected sites in Maryland, from Hill et al. (2003)

Shoreline type

Samples (N)

Locations (N)

Dry bulk density (g/cm3)

Bank

66

20

1.38

Marsh

22

5

0.62

Another important potential difference between erosion of banks and marshes may be a critical wave power threshold for erosion of banks. This may be indicated by the bank erosion behavior at low onshore wave powers shown in Fig. 11, as pointed out in the “Results” section above. A critical threshold for erosion has been incorporated into a model of soft cliff erosion proposed by Hackney et al. (2013), who suggested that accumulated wave power above a critical threshold is a better erosion predictor than wave power alone. Leonardi et al. (2016) point out that the most important consequence of the linear relationship (with no wave power threshold) in their proposed marsh erosion equation is that marshes should erode more in response to moderate wave activity than to extreme wave activity. This is simply because moderate waves are much more common than extreme waves. The steep negative exponential wave height probability distribution for CB waves shown in Fig. 6, combined with the linear relationships between CB wave power and marsh erosion rates shown in Fig. 12, supports this suggestion for CB marshes as well. However, if banks only begin significant erosion above some wave power threshold, then bank erosion would be much more affected by rare high wave events. A threshold dependence could be incorporated by replacing P by (P − Pcr) in Eqs. 6 and 7.

A final modification to Eq. 6 proposed here is to change the dependence of erosion rate on the ratio of bank height to average water depth, f(h/d) in Eq. 6, to a dependence on the ratio of time-varying water depth to bank height, g(D/h), where D is the wave-averaged time-dependent water depth and h is the bank height. The primary advantages of this approach are that it allows water depth to pass through 0 without a singularity in the value of the governing ratio, and that it allows erosion to vary temporally as a function of time-varying water depth. Halka and Sanford (2014) proposed such a functional dependence to ensure that sea levels that were either too high or too low would result in significantly less erosion through the influence of the structure function g. g must approach 0 for D/h < 0 because the water depth is below the base of the bank/marsh and the waves break offshore; g must approach 0 as D/h becomes greater than 1 and the waves pass over the top of the bank/marsh edge without breaking; and g ≈ 1 in the intermediate range 0 < D/h < 1. Tonelli et al. (2010) found a similar sea level dependence to be very important in their model of marsh edge erosion. This behavior is consistent with the inverse correlation between high sea level (as MedianWater60Hsig) and erosion rate shown in Table 2.

Combining all of these concepts, an expression for the mass rate of shoreline erosion is proposed here as a modification of the Marani et al. (2011) relationship (Eqs. 6 and 7):

where LHS is the mass rate of erosion, α = c−1 is the erodibility (inverse of cohesion), and Pcr is the critical wave power threshold for bank erosion. In practice, α would be treated as an empirical local constant of proportionality (an erosion rate constant) and Pcr as an empirical critical wave power. It is likely that Pcr = 0 for marsh edge erosion, but not for bank erosion. Equation 8 is directly analogous to one of the most widely used expressions for bottom sediment erosion (e.g., Sanford and Maa 2001). Full evaluation of Eq. 8 is beyond the scope of this paper, but it would be a worthy topic for future research. It should also be noted that Eq. 8 is only applicable to shorelines dominated by wave-forced erosion, since it does not account for the potential influence of riverine, tidal, or wind-forced currents. Such currents are generally too weak to cause shoreline erosion on their own in most nearshore environments of Chesapeake Bay (e.g., Halka and Sanford 2014), but current-induced suspension is implicitly assumed to work in concert with wave-forced nearshore erosion (e.g., Sanford 2008) to transport eroded sediments offshore.

Finally, it is important to acknowledge that this paper, like most other shoreline erosion studies, has focused on just one of the two components of shoreline erosion. The term “erosion rate” used in this study refers to the fastland erosion rate, which only accounts for soil eroded from above the waterline. Nearshore erosion, which occurs from the waterline to the base of wave action, must also occur for the full shore profile to migrate landwards. In much of the MD CB, total mass erosion rate is roughly a factor of 2 higher than the fastland erosion rate because of the addition of nearshore subaqueous material (Halka et al. 2005). Clearly, the actual ratio of nearshore erosion to fastland erosion depends on the geometry of the full shore profile, which will vary from site to site, but further discussion of this effect is beyond the scope of this paper.

Notes

Acknowledgements

The authors gratefully acknowledge the support of the US Army Corps of Engineers (ERDC), through contract number W912HZ-13-2-0002, as well as the support of Maryland Sea Grant through award number NA10OAR4170072. We thank our colleagues at the Maryland Geological Survey and US Army Corps of Engineers, who helped us gain access to computer codes, databases, and reports containing much of the information upon which this paper is based. We also think two anonymous reviewers for their helpful suggestions. This paper is dedicated to the memory of Evamaria W. Koch, friend, colleague, and fellow estuarine nearshore enthusiast. UMCES Contribution No. 5348.

Harris, C.K., J.P. Rinehimer, and S.-C. Kim. 2012. Estimates of Bed Stresses within a Model of Chesapeake Bay, Proceedings of the Twelfth International Conference on Estuarine and Coastal Modeling. Washington: American Society of Civil Engineers.Google Scholar

Leonardi, N., N.K. Ganju, and S. Fagherazzi. 2016. A linear relationship between wave power and erosion determines salt-marsh resilience to violent storms and hurricanes. Proceedings of the National Academy of Sciences 113: 64–68.CrossRefGoogle Scholar