Quantification of Swell Energy and Its Impact on Wetlands in a Deltaic Estuary

Abstract

The extensive coastal wetlands in Mississippi River Delta represent the seventh largest deltaic floodplain in the world, contributing to many services that sustain the economies of the region. Subsidence, sea level rise, saltwater intrusion, wave action from storms, and sediment depletion have contributed to chronic wetland losses, converting vegetated lands into open waters and increasing wind fetch. Among these factors listed, wave energy plays the largest role in marsh edge erosion in an open bay environment. Degrading barrier islands along the shoreline of this delta allow swell energy to enter protected bay areas, contributing to marsh edge erosion. Locally generated wind waves within enlarged bays also contribute to wetland loss. Quantifying the roles of swell and wind waves in marsh edge erosion is essential to any ecosystem restoration design. In this study, a numerical model is implemented to describe the wave climate of combined swell and wind waves in a deltaic estuary. Terrebonne Bay was chosen as the study area because it has experienced one of the largest reductions in barrier islands and wetland loss rates among Louisiana estuaries. A continuous wave measurement in upper Terrebonne Bay was obtained over the course of a year. A spectral wave model is used to hindcast the wave climate in the estuary. The model results are compared against the in situ wave measurement. The wave power is partitioned into swell and wind sea at different locations in Terrebonne Bay using the model results. An extensive analysis on a valid effective wave power range that directly impacts the marsh edge is performed and presented. Insight into the temporal and spatial variability of wave power is gained. Through differentiating swell and wind sea energies around the bay, improvements of long-term wave power computation for shoreline retreat prediction are made. It is found that the swell energy becomes the primary driver of marsh edge retreat in the southwest part of Terrebonne Bay as the barrier islands are degrading.

Introduction

The Mississippi River Delta consists of extensive coastal wetlands within deltaic estuaries on Louisiana’s coastline. Coastal wetlands provide significant defense mechanism against hurricane storm surges, limiting their propagation inland, reducing threats to infrastructure and the lives of both people and wildlife (Hu et al. 2015; Wamsley et al. 2010). Such services of coastal wetlands are currently under great strain from a variety of coastal proccesses including relative sea level rise, saltwater intrusion, storms, and sediment deprivation. Wetland loss is not only a result of reduced sediment delivery to the detlaic floodplains, but is also caused by the erosive power of waves along the more exposed edge of wetlands in more exposed bays. Chronic wetland losses have converted vegetated lands into open waters and increased wind fetch. This increase in wind fetch generates larger waves, which in turn leads to greater erosive forces that impact the soil stability of wetlands.

The deltaic estuaries along Louisiana’s coastline were formed over the last 6000 years by the Mississippi River naturally changing course, depositing large amounts of sediment across the Louisiana coast known as delta switching (Coleman et al. 1998). Such deltaic processes resulted in wetland dominated estuaries with minor portion of the estuary made of bay water area. However, the need for navigation and flood protection increased the construction of levees along the river reducing sediment delivery to surrounding wetlands. Such interruptions of deltaic processes lead to degrading wetlands as a direct result of the anthropogenic riverine management strategies, severely limiting sediment supply delivered to deltaic estuaries. As a result, wetland areas declined, bay area increased, and barrier islands continued to degrade (Twilley et al. 2016). While the increased wind fetch leads to greater locally generated wind wave energy, the reduction in barrier islands increases swell energy penetrating into the estuaries. The increase in swell energy contributes to erosion processes along exposed marsh edges that are also subject to the impact of locally generated wind sea waves. On open coasts, wind sea and ocean swell waves often occur simultaneously. Most of the wave studies in estuaries have been focused on locally generated wind waves (e.g., Chen et al. 2005; Mariotti and Canestrelli 2017) and ocean swell has been less investigated in estuaries. Therefore, improving the understanding of swell contribution to the total wave energy budget in degrading deltaic estuaries and its impact on marsh edge retreat rates is paramount for mitigation of future land loss and storm impact on coastal communities.

Louisiana has lost approximately 60 km² of wetlands each year from 1932 to 2010 (Couvillion et al. 2011). Due to high rates of relative sea level rise (subsidence rates in most of coastal Louisiana is higher than eustatic sea level rise), the reduction of sediment input to deltaic floodplains has focused on wetlands drowning as major cause of wetland loss in deltaic estuaries (Boesch and Turner 1984; Lee and Gosselink 1988). Less attention has been paid to erosion along marsh edge boundaries resulting from wave action as a causative agent of coastal wetland loss. In an analysis performed by Penland et al. (2000), it was determined that 26% of the wetland loss in the Mississippi River deltaic estuaries can be attributed to erosion due to wind waves. For wetland restoration and coastal management, it is crucial to understand how meteorological conditions and hydrodynamic processes are functioning and interacting with each other in a coastal setting with unique hydrological processes and geomorphologic features.

Over the past couple of decades, many types of numerical models for coastal wave processes have been developed. The spatial scales of these models range from deep water, phase-averaged wave models to nearshore, phase-resolving models. The phase-resolving models can solve the Boussinesq-type equations, e.g., Chen et al. (2000), but are restricted to regions of limited size since they require such fine resolution both spatially and temporally. For larger scale applications, such as an estuary, spectral phase-averaged models are the most appropriate (Battjes 1994; Qin et al. 2005). Deep-ocean or shelf-sea waves can be predicted well using third generation wave models based on the energy or action balance equation, such as the WAve Model (WAM) (WAMDI1988) and WAVEWATCH model of Tolman (1991). However, these models cannot be realistically applied to nearshore regions of estuaries, tidal inlets, barrier islands, or tidal flats with spatial scales less than 20–30 km and water depths less than 2–3 m (Booij et al. 1999). Booij et al. (1999) developed the SWAN model to transition the shallow water formulations from deep water processes to shallow water through the incorporation of (1) a shallow-water phase speed in the expression of wind input, (2) a depth-dependent scaling of quadruplet wave-wave interactions, (3) reformulation of the whitecapping in terms of wave number rather than frequency, (4) adding bottom dissipation, (5) depth-induced wave breaking, and (6) triad wave-wave interactions. Chen et al. (2005) were among the firsts to model wind waves in an estuary considering the effects of water level change and currents, although they neglect the effects of swell. In the present study, through a detailed numerical model simulation, locally generated wind sea power and penetrating swell energy are quantified at numerous locations near the exposed marsh boundaries in Terrebonne Bay, Louisiana to better understand their temporal and spatial variability within the bay.

It has been well documented that wave energy flux, also known as wave power, serves as a good proxy for determining retreat rates on marsh shorelines. Schwimmer (2001), Marani et al. (2011), and McLoughlin et al. (2015) have all contributed to correlating the wave power to marsh edge retreat. In previous studies (Marani et al. 2011; McLoughlin et al. 2015), parametric models such as Young and Verhagen (1996) were used to calculate wave power within the boundaries of estuaries to compare with salt marsh edge retreat rates. While this approach may be applicable to smaller spatial scale studies (i.e., on the order of hundreds of meters of shoreline), it is not suitable for estuary-scale studies, such as this one, due to the complex shoreline geometry found throughout. Moreover, limitations for using the Young and Verhagen (1996) equations to hindcast waves include fetch-averaged depths, quasi-steady wind inputs, and the assumption that the wave is traveling directly along the fetch. In this study, average wave power values calculated using the SWAN spectral wave model are compared to marsh edge erosion rates established by Allison et al. (2015) at multiple locations in Terrebonne Bay.

The objective of this study is threefold: (1) to implement and validate the SWAN wave model to accurately predict the wave climate in Terrebonne Bay, a shallow water estuary with complex geomorphology; (2) to determine the effective range of wave power that directly impacts the wetlands in this estuary; and (3) to quantify the magnitude of swell energy along the shoreline in Terrebonne Bay and compare it to locally generated wind wave energy that affects marsh edge erosion rates. Advancements to existing studies include improving the accuracy of wave power estimates by using the SWAN model, tested against field measurements, to account for spatially varying fetch and bathymetry along with temporally varying wind speeds and directions as well as swell energy from the Gulf of Mexico. Incorporating the effects of swell can more accurately represent the erosive wave power near the marsh edge boundary, especially near the entrance of the bay. The methodology of quantifying and partitioning the swell and wind wave effective power can be used to study marsh edge erosion in other estuaries along the northern Gulf Coast and beyond.

Methods

Study Area

The study site is in Terrebonne Bay (hereafter, TB), which is located west of the Mississippi River in the northern Gulf of Mexico just west of Barataria Bay (Fig. 1). TB is chosen for this research because it has experienced one of the largest wetland loss rates among Louisiana estuaries. Originally, TB was formed as old delta lobes of the Mississippi River, the Teche and Lafourche deltas, subsided and deteriorated (Wang and Wang 1993). During the period of 1932–2010, the Terrebonne Bay Basin lost 3082 km², making up 25% of total wetland lost in coastal Louisiana during this time (Couvillion et al. 2011). Today, the bay receives virtually no fluvial inflow, depriving it of a sediment influx that it vitally needs.

Fig. 1

Terrebonne Bay, Louisiana (LA) study site showing bathymetry (NAVD88), SWAN domain coverage and orientation, and model input source locations

Terrebonne Bay is on average about 40 km wide and 25 km long and has an average depth of 1.7 m with a maximum depth of 3.0 m. It is surrounded by wetlands and bounded offshore in the south by barrier islands that are approximately 10 km long and have an 11 km gap between them. Locally generated wind waves dominate the wave climate in TB; however, some swell enters the bay from the Gulf of Mexico. Further, TB experiences a diurnal, microtidal range varying from 0.1 to 0.2 m during equatorial tides and from 0.3 to 0.6 m during tropic tides (Leonard and Luther 1995).

Wave Measurements

Detailed wave measurements were collected for the duration of a full year in TB (Parker 2014) with an OSSI (Ocean Sensors Systems Inc., OSSI-010-003C) pressure transducer (or wave gauge) that records a time series of water pressure at a deployment depth. These instruments were selected because of their ability to record data at high sampling rates for a long duration of time. This allows for a continuous collection of wave measurements over the entire study period in TB. For this study, the bottom-mounted wave gauge was deployed at 29°13′21.07″ N, 90°36′24.56″ W and set to record at a rate of 10 Hz, for 20-min bursts in every 30 min, resulting in 48 data sets per day. The high sampling frequency of the instrument allows the short wind waves that occur during this study period to be resolved. The instrument was replaced at the same location approximately every 3 months to obtain a continuous wave record spanning the entire year of 2012. However, when analyzing the field data, it was noted that the water level record started to show unreliable measurements (drift) toward the end of August 2012. Therefore, the data set for analysis was truncated after August 20th, 2012. The analysis of the measured data resulted in multiple parameters such as a time series of the wave height, wave period, and water depth, which are all used to compute wave energy and wave power (see spectral analysis for detail).

Spectral Analysis

Spectral analysis is used to convert raw water pressure data into desired wave parameters since it can describe the distributions of irregular wave energies over a wide range of wave frequencies, following methods outlined in Kamphuis (2010) and Karimpour and Chen (2016). After the dynamic part of the measured pressure from the wave gauge was converted to water level, P/ρg, wave properties were calculated from a water surface elevation power spectral density, S_ηη, which is represented as:

$$ {S}_{\eta \eta}=\frac{1}{K_p^2}\times \frac{S_{pp}}{\rho^2{g}^2} $$

(1)

$$ {K}_p(f)=\frac{\cosh \left(k{d}_p\right)}{\cosh (kh)} $$

(2)

where K_p is the dynamic pressure to surface elevation conversion factor, ρ is the density of water, g is the gravitational acceleration, S_pp is the wave dynamic pressure spectrum, f is the frequency, h is the local water depth, k is the wave number, and d_p is the pressure measurement distance from the bed.

The process of using the spectral density analysis for the pressure data series requires multiple steps. First, a fast Fourier transform (FFT) is applied to transform the water surface signal from time domain to frequency domain. Next, for each frequency in the frequency domain, a wave number is calculated using the dispersion relationship in linear wave theory, and then, K_p is calculated for each frequency. The pressure conversion factors, K_p, are applied to the energy density spectra in frequency space to compensate for the wave dynamic pressure attenuation in the water column. It should be noted that K_p becomes close to zero in high frequency which leads to exaggerated energy values for those high-frequency waves, and therefore, its application must be cutoff at an empirically determined frequency to ensure that high-frequency energies and any presented noise are not unrealistically magnified. A simple high-cutoff frequency f_max can be used to remove this portion of the data set (Karimpour et al. 2016). The high-cutoff frequency in this study was set equal to f_max = 1 Hz, whereas the low-cutoff frequency, was set equal to f_min = 0.05 Hz. A smooth energy spectrum was obtained by using the Welch’s method (Welch 1967). In this study, each burst is divided into 16 overlapping segments with 50% overlap (Karimpour and Chen 2016). Each of these segments is windowed by using the Hanning window function (Oppenheim et al. 1999). After the spectrum was calculated, K_p was applied to it up to a frequency of 1.0 Hz to follow the high-cutoff frequency limit (Fig. 2).

Fig. 2

Example of wind sea-dominated (blue) and swell-dominated (orange) frequency energy spectrum for a single burst with pressure correction applied. Also shown is the calculated sea-swell partitioning frequency

Then, the sea and swell energies were seperated using the method from Hwang et al. (2012), as

$$ {m}_n(f)={\int}_f^{f_u}f{\prime}^n{S}_{\eta \eta}\left({f}^{\prime}\right){df}^{\prime } $$

(3)

$$ {I}_1(f)=\frac{m_1(f)}{\sqrt{m_{-1}(f)}} $$

(4)

$$ {f}_s=24.2084{f}_m^3-9.2021{f}_m^2+1.8906{f}_m-0.04286 $$

(5)

where m_n is the nth moment of the wave energy spectrum, f_u = 0.5 Hz, f_m is the frequency associated with the maximum value of the spectrum integration function, I₁, and f_s is the separation frequency. The left spectrum in Fig. 2 shows a typical case where the sea energy, locally generated by wind, is dominating the energy spectrum. In contrast, the right spectrum in Fig. 2 shows a case where the swell energy dominates the spectrum, typically during calm winds. The zero-moment wave height, or significant wave height, H_s, is calculated for both sections of the spectrum using the zeroth moment, m₀, of the energy spectrum as:

$$ {H}_s=4\sqrt{m_0} $$

(6)

$$ {m}_0={\int}_{f=0}^{f=\infty }{S}_{\eta \eta}(f) df $$

(7)

Moreover, the peak wave period is determined for the sea and swell energy by finding the highest point in the energy spectrum, locating the frequency associated with it, and then simply inversing this frequency. Note that this method works well when dealing with a well-defined spectrum as shown in Fig. 2. However, when the overall wave energy is low, the point defined as the peak tends to shift to an incorrect frequency because of the abnormal shaped, peaky spectral density plot. This usually occurs when small waves, typically less than 5 cm, are present. Ultimately, the observed wave parameters of the significant wave height, H_s, peak wave period, T_p, and water level are used to calculate wave power and validate the model.

Nearshore Wave Conditions

The SWAN model is a third-generation wave model used to compute waves in coastal areas from given wind, bottom, and current conditions. It is based on an Eulerian formulation of the discrete spectral balance of action density that accounts for changes in wave height, shape, and direction over arbitrary bathymetry and current fields (Booij et al. 1999). Spectral wave models are based on the wave action balance equation, which is represented as:

$$ \frac{\partial }{\partial t}N+\frac{\partial }{\partial x}{c}_xN+\frac{\partial }{\partial y}{c}_yN+\frac{\partial }{\partial \sigma }{c}_{\sigma }N+\frac{\partial }{\partial \theta }{c}_{\theta }N=\frac{S}{\sigma } $$

(8)

where the first term on the left hand side is the local rate of change of action density (N) in time, the second and third terms represent propagation of action in geographical space with c_x and c_y as the propagation velocities in x and y space. The fourth term represents the shifting of the relative frequency (σ) due to variations of depths and currents where θ is the wave direction. The fifth term represents depth- and current-induced refraction. On the right hand side, S is the source term in representing the effects of generation, dissipation, and nonlinear wave-wave interactions.

In this study, the SWAN model is set up to simulate the wave conditions mainly driven by local wind and partly by ocean swell in TB. The 2D nonstationary mode of SWAN version 40.91 is implemented in TB using the backward space backward time finite difference scheme. The physics include quadruplet interactions, triad interactions, whitecapping, breaking, bottom friction, and diffraction. The van der Westhuysen et al. (2007) wave generation is activated. This formulation is based on experimental findings that whitecapping dissipation appears to be related to the nonlinear hydrodynamics within different wave groups (SWAN Team 2016). This yields a dissipation term that depends on quantities that are local in the frequency spectrum, as opposed to ones that are distributed over the entire spectrum, as found in the default formulation of Komen et al. (1984) (SWAN Team 2016). The calculation for depth-induced wave breaking is a spectral version of the Battjes and Janssen (1978) formulation. A constant breaker parameter, γ = 0.73, is used in the breaking simulation.

Bottom friction is one of the physical mechanisms for wave energy dissipation in shallow water systems. In this study, the semi-empirical expression derived from the JONSWAP (Joint North Sea Wave Project) model of Hasselmann et al. (1973) is chosen for bottom friction dissipation. This JONSWAP formulation can be expressed in the following form:

$$ {S}_{bfr}\left(f,\theta \right)=-{C}_bE\left(\sigma, \theta \right){\sigma}^2/\left({g}^2{\sinh}^2 kd\right) $$

(9)

where C_b is the bottom friction coefficient, σ is the relative angular frequency, and θ is the propagation direction. In a study by Alkyon (2009), decreasing the bottom friction coefficient is a diagnostic modeling measure taken to help eliminate the underestimation of swell energy near the mainland coastline in the Eastern Wadden Sea. Additionally, van Vledder et al. (2010) found that for regions with typical sandy bottoms, the best wind sea hindcast results were obtained when the default bottom friction value of C_b = 0.067 m²s⁻³ was replaced by a lower value of 0.038 m²s⁻³. However, for smoother seafloors like ones in the Gulf of Mexico, a lower value of 0.019 m²s⁻³ is advised (SWAN Team 2015). Bottom friction values ranging from 0.010 − 0.067 m²s⁻³ are tested on the model. The values used in this study are 0.038 m²s⁻³ and 0.021 m²s⁻³ for wind sea and swell conditions, respectively.

In this model setup, the spectral densities cover the full directional range of 360°. The resolution in directional space is set to 10° for wind sea and 5° for swell conditions. However, the resolution in frequency space is not constant. It is defined in the SWAN manual as:

$$ \frac{\Delta f}{f}={\left(\frac{f_{high}}{f_{low}}\right)}^{1/n}-1 $$

(10)

in which f_high is the highest discrete frequency used in the calculation, f_low is the lowest discrete frequency, and n is one less than the number of frequencies. It must be noted that the DIA approximation of the quadruplet interactions in SWAN is based on a frequency resolution of ∆f/f = 0.1. For this model, f_low is set to 0.05 Hz and f_high is set to 1.05 Hz resulting in an n value of 32.

The model domain can be seen in Fig. 1. A high-resolution mesh is necessary to accurately resolve the complex shoreline geometry and nearshore bathymetry found in this region. A structured, rectangular mesh is used in this study with a 50-m resolution. The wave boundary condition is applied along the lateral and seaward boundaries, whereas the water level and wind are applied uniformly over the domain. Additionally, bathymetric data with a spatial resolution of 30 m are used from USGS EROS (2015) and interpolated onto the mesh.

The SWAN model is forced with wind and water level inputs from nearby National Oceanic and Atmospheric Administration (NOAA) monitoring stations. The Grand Isle, LA station, “GISL1”, located at 29°15′54.00″ N, 89°57′28.80″ W is used as the water level boundary condition in this study since it had the most consistent data set when compared with nearby stations for the year of 2012. For wind forcing, the “LUML1” station located at 29°15′17.99″ N, 90°39′50.42″ W inside TB is used in the SWAN model. The hourly wind speed obtained from this station is a 2-min scalar average of 1-s wind speed measurements and reported on the hour. The wind and water level forcing are applied uniformly over the entire model domain. Occasionally, the LUML1 station has gaps in the data set as a result of the passing of a severe cold front or hurricane. In order to fill these gaps to have continuous yearlong wave records, 3-hourly 32-km resolution winds from NOAA’s National Center for Environmental Protection (NCEP) North American Regional Reanalysis (NARR) were simulated through the model domain and output at the LUML1 station location. The NARR output winds were then used to fill the missing gaps in the raw LUML1 data.

Finally, an offshore wave boundary condition is necessary to incorporate the swell energy into the SWAN model. Wave data along the seaward boundary of the model domain was obtained from the United States Army Corps of Engineers (USACE) Wave Information Studies (WIS) (http://wis.usace.army.mil/). WIS has consistent, hourly, long-term wave climatology along all US coastlines and US island territories. The numerical model used to provide the results in the Gulf of Mexico is WISWAVE. WISWAVE is a discrete spectral wave model solving the energy balance equation for the time and spatial variation of a 2D wave spectrum. The framework of this code is derived from Resio (1981) and modified to include shallow-water effects by Hubertz (1992). The WIS website provides access to the database of hindcasted wave information for a series of “virtual wave gauges” in water depths ranging from 15 to 20 m. Wave properties of significant wave height, peak wave period, mean wave direction, and directional spreading are available for total, sea, and swell energy components. A JONSWAP one-dimensional spectrum with a certain imposed directional distribution is generated by SWAN from the USACE WIS wave parameters at the boundary. There are four WIS stations available at the seaward boundary of the SWAN domain (Fig. 1). It is determined that there is negligible difference between significant wave height and peak period reported by these four stations throughout the yearly time series. Therefore, using only station 73125 uniformly across the offshore boundary is justified. For the lateral boundary conditions, the sea and swell significant wave height is linearly interpolated from its original, deepwater, magnitude to zero at the land boundary. The purpose of this is to eliminate the lateral boundary effect. The amount of error that the boundary effect can induce depends on the incoming wave direction and directional spreading of the wave field. However, it must be noted that there is a slight error associated with the implementation of a linearly interpolated lateral boundary condition. In reality, waves shoal and refract as they approach the shoreline. Using the interpolated boundary conditions masks the effects of shoaling along the boundary and ignores refraction since a constant wave direction is assigned uniformly along the boundary. The effect of such a boundary error on the study area is minimized by extending the lateral boundaries away from the area of interest.

Classification of Bulk and Effective Wave Power

Clarification for the multiple wave power parameters presented in the following sections is necessary. The “combined wave power”, P_comb, is the summation of the wind sea and swell wave power values. The “bulk wave power”, P_b, is the raw wave power calculated at each time step, regardless of whether the water level exceeds the marsh surface elevation or falls below the toe of the marsh edge, or without using a screening or filtering scheme. The “effective wave power”, P_e, on the other hand, is the selected or filtered wave power value calculated using a range of water levels as defined later in this section. The bulk and effective wave powers are calculated for wind seas, swell, and combined conditions (P_b,sea, P_b,swell, P_b,comb, P_e,sea, P_e,swell, and P_e,comb).

Water levels are known to have a direct effect on the magnitude of wind waves. Through numerical and field studies, wave action is found to have the harshest impact on wetlands when the water level is just below the marsh platform (McLoughlin et al. 2015; Tonelli et al. 2010). Moreover, it has been found that wave action increases with water depth, up until the marsh becomes submerged, then rapidly decreases in magnitude (Leonardi et al. 2016b). Therefore, it is important to calculate the effective wave power because the marsh edge is considered safe from wave attack when water levels exceed a certain elevation above the marsh platform, typically during the passing of a cold front or hurricane, or when the water level is below the scarp bottom. An appropriate range of water level used for effective wave power calculation is explored herein.

Based on findings from Tonelli et al. (2010), McLoughlin et al. (2015) calculated the effective wave power by setting P_i = 0 when the instantaneous water level was above the marsh scarp top elevation. However, this definition excludes the effects of waves breaking onto the marsh edge and the effect that wave orbital velocities can have on the inundated marsh surface or mudflat (Karimpour et al. 2016). Furthermore, in order to establish these marsh elevations, a detailed topographic survey is needed at all 125 model output locations in TB (See Bulk and Effective Wave Power section for model output details). Because this was not feasible for an estuary-scale study, mean high water (MHW) and mean low water (MLW) are used as a proxy of the elevations of the marsh platform and marsh toe, respectively. The marsh platform elevation, for stable salt marshes, exists at an elevation within the intertidal zone that approximates that of the MHW (Krone 1985; Morris et al. 2002). Therefore, a new water level range for effective wave power calculation is proposed using these representative elevations.

The time series of the significant wave height, H_s, and peak wave period, T_p, output from the SWAN model, along with the water level time series, are used to calculate wave power, P, at numerous locations across TB:

$$ P=\frac{\rho g{H}_s^2}{8}{C}_g $$

(11)

where C_g is the wave group velocity and ρ is the water density. The upper water level limit is analyzed for 0, 1.0, and 1.5 ×H_s above the MHW and the lower water level limit is set as 1.0 ×H_s below the MLW. Figure 3 shows a schematic of the new water level effective range for 1.0 ×H_s above MHW. P_e is calculated when the water level is between this range using Eq. (11); otherwise, it is set to zero at that time step. This approach is applied to both the wind sea and swell SWAN results, and then added together to obtain P_e,comb in TB. Additionally, when the modeled values of significant wave height for wind sea, H_s,Sea, are below 5 cm, then all P_e,Sea values were set to zero because very small waves will not have any significant influence on the marsh edge. Moreover, SWAN overestimates the wave energy for these small waves; therefore, the results are not reliable. Note that this adjustment was not applied to P_e,Swell.

Fig. 3

A schematic of the proposed effective range for wave power computation

Five different years of wave conditions are simulated to compare the effects of hurricanes and normal conditions on the wave power magnitude across the estuary. The years 2004 and 2006 are modeled since they were relatively calm years for TB that included no hurricanes passing through the Louisiana coast. The year 2005 includes Hurricane Katrina, which made landfall east of TB, and Hurricane Rita, which passed west of the bay. Similarly, the year 2008 includes Hurricane Gustav, which passed through TB, and Hurricane Ike, which passed west of the bay. The year 2012 includes the passing of Hurricane Isaac, which made landfall just east of TB.

Results

Temporal Characteristics of Wave Power

Wave power values calculated from the SWAN model results are tested against the observed wave power at the wave gauge deployment site. A comparison of the observed and modeled daily-averaged bulk wave power for wind sea conditions (Fig. 4) and for swell conditions (Fig. 5) is shown. Very good agreement is found for daily-averaged bulk sea, P_b,sea, while the swell wave power agreement, P_b,swell, is reasonable. Note that in these figures, the vertical scale for the period of March 19th–June 6th is different from the other periods to better show the complete wave power variation of the time series. The seasonality of the wave power is clearly observed in both temporal distributions. P_b,sea is greatest in mid-March, reaching peaks of approximately 550 W/m. Furthermore, it is important to note that the modeled and observed values are well correlated for these large values of P_b,sea. For the daily averages of P_b,swell, the magnitude fluctuates consistently through first quarter of the year (Fig. 5). These pulses in high P_b,swell are a direct result of storms passing across the Gulf of Mexico. Offshore NOAA monitoring buoys “SPLL1” and “LOPL1” show high significant wave height values ranging between 1.8 and 2.4 m for March 21st, 2012, and values up to 1.7 m on June 24th, 2012. Additionally, although the general temporal trend of the P_b,swell is modeled well, for some of the larger swell events, there is a noticeable difference between the modeled and observed P_b,swell. The model slightly overpredicts swell wave power with a bias of 1.98 W/m at the measurement location in upper TB. Lastly, swell energy is less affected by local winds; however, the highest values of the P_b,swell and P_b,sea occur simultaneously.

Fig. 4

Daily averages of observed and modeled bulk wind sea wave power, P_b,sea, at the wave gauge deployment site

Fig. 5

Daily averages of observed and modeled bulk swell wave power, P_b,swell, at the wave gauge deployment site

A scatter plot of the observed and modeled daily-averaged, P_b,sea and P_b,swell, is shown in Fig. 6. It is seen that the wind sea wave power is generally higher than the swell wave power at the field observation site in upper TB. The daily averages of P_b,sea are well correlated with the observed values (R² = 0.86, RMSE = 27.1 W/m). The daily averages of P_b,swell are reasonably predicted by the model (R² = 0.54, RMSE = 6.4 W/m). The slight overprediction in P_b,swell previously discussed is reflected in Fig. 6b.

Fig. 6

A comparison of observed and modeled daily averaged bulk wave power values at wave gauge deployment site for wind sea (a) and swell (b) conditions

Bulk and Effective Wave Power

In addition to the temporal variation of wave power at the study site, it is important to understand how this parameter varies spatially across TB. Given that the model has been tested against the field data for sea and swell wave power, these values can be calculated at multiple locations inside the bay. Due to the size of the computational domain and the grid resolution, the model results were output at selected sites rather than over the whole domain. One hundred twenty-five (125) locations are chosen for the wave power analysis. These sites correspond to the location of marsh edge retreat rates established in TB by Allison et al. (2015). To resolve waves in the nearshore region, where the retreat rates are established, a phase-resolving model, such as the Boussinesq wave model (Chen et al. 2000), is required to accurately simulate the complex phenomena that occur here. Since SWAN is a phase-averaged model, the output locations are selected at a water depth of approximately 1 m, similar to the study site where SWAN was shown to run stably and produce accurate results. The output locations are extracted 100 m normal from the marsh edge; however, some of the points along this contour still showed shallow depths around 35 cm. Therefore, these points are shifted further offshore until a depth near 1 m was found. The average depth of all 125 model output locations is 1.04 m, with the furthest point approximately 200 m from the shoreline.

Five years (2004, 2005, 2006, 2008, and 2012) of wave conditions are modeled to compare the effects of hurricanes and normal conditions on the combined bulk, P_b,comb, and combined effective, P_e,comb, wave power magnitude in TB. A sensitivity analysis is performed on the three effective water level ranges of wave power discussed earlier. Annual averages of P_e,comb and percent difference between P_b,comb and P_e,comb, considering the three effective ranges, are performed at all 125 observation locations. The minimum, maximum, spatial average of all the sites (mean), and standard deviation (std) of these annual averages are shown for P_e,comb (Table 1) and the percent difference (Table 2), determined using the following relationship:

$$ \% Diff=100\ast \frac{\left|{P}_{b, comb}-{P}_{e, comb}\right|}{\frac{P_{b, comb}+{P}_{e, comb}}{2}} $$

(12)

Table 1 Values of combined effective wave power, P_e,comb (W/m), for five different years, each year considering three effective ranges

Table 2 Values of percent difference (%) between the combined bulk, P_b,comb, and combined effective, P_e,comb, wave power for five different years, each year considering three effective ranges

The upper limit of the effective range is specified in both tables and the lower limit is defined as 1.0 × H_s below MLW.

Previously, McLoughlin et al. (2015) calculated average effective wave power when the water level was at or below the marsh platform (i.e., 0.0 × H_s above MHW in this study). It is observed in Table 2 that using this water level range can filter out up to 121% P_b,comb and on average cause a 66.6% reduction in wave power. Additionally, it is important to consider the effects of wave breaking and orbital velocities on the marsh edge (Karimpour et al. 2016), which would be lost when using this range. Alternatively, using an effective range of 1.5 × H_s only shows up to 14% reduction on average (Table 2), including too much of the bulk wave energy. Tonelli et al. (2010) found that wave thrust acting on the marsh edge is significantly reduced as depth above the marsh platform increases. Considering this range will likely overestimate wave power acting on the marsh edge, including wave power values up to 175.6 W/m (Table 1). According to Tonelli et al. (2010), the harshest impact of wave thrust and energy dissipation on the marsh edge occurs when the water level is close to, or immediately above, the platform elevation. Therefore, we conclude that 1.0 × H_s above MHW is the appropriate upper boundary for effective wave power computation and is used for the remainder of this study. Using this effective range, note that P_e,comb values are nearly identical for 2004, 2005, 2006, and 2012, with values of 36.4, 35.1, 34.2, and 36.8 W/m, respectively (Table 1). Those years do not vary significantly from the 2008 year which has a P_e,comb value of 40.1 W/m (Table 1). Moreover, because it is difficult in practice to obtain the exact marsh platform elevations on a basin scale, using 1.0 × H_s above MHW as the limit for the effective wave power calculation takes into account a range of platform elevations that at some locations in TB may lie above MHW.

A spatial distribution of the annual average P_e,comb is shown in Fig. 7 for 2004, 2005, 2006, 2008, and 2012. Although results were obtained at all 125 locations, P_e,comb values are only shown for 14 locations in the figure for viewing clarity. The magnitudes of P_e,comb at these locations are representative of the other nearby locations included in the analysis. The geographic settings of these output locations vary significantly from being directly exposed to open bay waters, to being sheltered by the surrounding tracts of marsh on intertidal mudflats. The largest values of wave power are found at the western part of the bay, near the opening to the Gulf of Mexico (Fig. 7). Moreover, the output locations sheltered by other land masses have the lowest wave power values. These locations are characterized by small fetch conditions which, in turn, limit wave growth. The locations in the northern region of the bay that are directly exposed to open waters and long fetches (approximately 30 km) have intermediate wave power magnitudes. Lastly, a one-way analysis of variance (ANOVA) was performed to test for differences in annual averaged P_e,comb among the 5 years. Results from the test (p = 0.092 > α = 0.05) show that P_e,comb was not significantly different between the 5 years at all of the sites in TB. Therefore, it can be concluded that calculating this wave power value for only 1 year is sufficient. For the remainder of this study, all wave power values presented will be an average of the 5 years that were modeled.

Fig. 7

Spatial distributions of combined annual effective wave power in Terrebonne Bay, Louisiana. Fourteen representative locations are shown for viewing clarity

Spatial Distributions of Wave Power and Erosion Rates

The swell wave energy is investigated to determine its contribution to the total energy budget and to better understand how it varies spatially across TB. To quantify this relationship, the swell effective wave power, P_e,swell, is calculated from the SWAN model results, and then divided by the wind sea effective wave power, P_e,sea, to obtain the percentage, i.e., P_{e, swell}/P_{e, sea}. The values above 100% are interpreted as being dominated by swell. The spatial distribution of the annual average P_{e, swell}/P_{e, sea} at all 125 sites in TB demonstrates that the swell wave power is most dominant near the entrance of the bay, especially on the west side (Fig. 8).

Fig. 8

Spatial distribution of annual, swell effective wave power percentage (P_swell/P_sea) in Terrebonne Bay, LA. A wave rose is included to show the magnitude and direction of swell energy at the offshore boundary

One reason for the high P_e,swell at the southwest entrance of the bay is because of the dominant, southeasterly, swell direction at the boundary of the domain (Fig. 8). The magnitude of P_e,swell along the southwest entrance ranges from 5 to 250% of P_e,sea (1–94 W/m) and to the southeast 5 to 15% (1–4 W/m). By contrast, P_e,swell averages approximately 0.5% of P_e,sea (~< 0.5 W/m) in the northern portion of the bay. Swell energy is reported to undergo extensive dissipation as it propagates across complex bathymetry and extensive tracts of shallow water (Talke and Stacey 2003). This significant reduction is reflected in Fig. 8. As the swell energy propagates into these shallow waters, the interaction between the wave and the bottom friction increases, causing swell dissipation. Moreover, bottom friction has a greater effect on swell waves because of their long wave lengths and wider coverage of orbital velocities (Sorensen 2005).

The spatial distributions of effective wave power are compared to marsh edge retreat rates in the bay. Figure 9 shows the spatial distributions of the annual averaged P_e,comb and established retreat rates determined for the period of 2004–2012 by Allison et al. (2015). Similar to Fig. 7, Fig. 9 shows only 44 locations in the figure for viewing purposes. Since it was previously concluded that using the average of the 5 years to represent typical annual wave power values in TB, it is acceptable to compare this value against retreat rates, R, that were established over nearly a decade. The wave power and retreat rates are normalized by the spatial average of all the sites, i.e., ${P}_e^{\ast }={P}_{e, comb}/{\overline{P}}_{e, comb}$ and ${R}^{\ast }=R/\overline{R}$. Higher R^∗ values are matched with larger magnitudes of ${P}_e^{\ast }$ on the open coast, and the lower values of R^∗ correspond well to the smaller magnitudes of ${P}_e^{\ast }$ in the more sheltered, lagoon-type locations (Fig. 9). The linear regression does show that at 0 effective wave power, which suggests that there could be some small amount of edge erosion in sheltered areas from wind waves. For low wind wave energy conditions, there is still the possibility for erosion from tides and boat wakes. Additionally, where ${P}_e^{\ast }$ is significantly lower than R^∗, processes that are not accounted for in this study, such as marsh submergence or direct removal, could play an important role of wetland loss at these locations.

Fig. 9

Comparison of dimensionless combined effective wave power, P*, and dimensionless retreat rate, R*, in Terrebonne Bay, LA

The improvement of incorporating swell energy into wave power computation is demonstrated by plotting R^∗ against ${P}_e^{\ast }$ in TB, which is similar to the relationship R^∗ = 0.67P^∗, established by Leonardi et al. (2016a) (Fig. 10a). This relationship developed by Leonardi et al. (2016a) was obtained by averaging numerous data points from 8 different bays across the world including 6 in the USA. The best fit line, R^∗ = 0.64P^∗ + 0.36, is included to show that the general trend is similar to the relationship reported in Leonardi et al. (2016a). Note that the data are categorized by location inside of TB. The first group includes locations with high retreat rates, R^∗ > 2.1, for ${P}_e^{\ast}\sim 1.25$. The second group consists of locations influenced by swell, which is classified as the P_e,swell > 5% of P_e,sea. The rest of the data points are in the group labeled “Other”. Alternatively, without including the effects of swell, the best fit line deviates significantly from the Leonardi et al. (2016a) relationship (Fig. 10b). The ${P}_e^{\ast }$ values shown in this figure are only considering the effects of wind sea energy; the points labeled “Influenced by Swell” are just shown as a means of comparison to Fig. 10a. From these figures, it is clearly observed that the incorporation of swell energy in the wave power computation improves the general trend of the relationship between the normalized wave power and normalized marsh edge erosion.

Fig. 10

Relationship between dimensionless combined wave power (P*_e) (a), dimensionless sea wave power (P*_{e, sea}) (b), and dimensionless erosion rate (R*). The dashed line shows best fit for the data and the solid line represents the relationship, R* = 0.67P*, determined by Leonardi et al. (2016a). The mean absolute combined effective wave power, sea effective wave power, and retreat rate are 36.7 W/m, 32.5 W/m, and 4.4 m/year, respectively

The geographic locations of the point classifications used for Fig. 10 are shown in Fig. 11. A total of 39 points at the southwest and 3 points at the southeast entrance to the bay are influenced by swell energy from the Gulf of Mexico.

Fig. 11

Geographic location of data point classification shown in Fig. 10a, b

Soil properties such as shear strength and cohesiveness, marsh cover and age, root volume, and species distribution are all factors that should be accounted for to understand how the strength of the marsh varies spatially. Figure 11 shows a total of 12 locations around the bay with areas of high retreat rate, R^∗~2.1 − 6.8, but have a constant wave power, ${P}_e^{\ast}\sim 1.25$. For these locations, it can be deduced that the soil and vegetation properties are likely to be weaker here compared to other locations. The high retreat rate points correspond to areas of fragmented marsh (Fig. 11), backed by small open-water bays, indicating that the marsh islands at these sites are highly unstable. The marsh fragmentation also reduces the wetland ability to trap or retain sediment to keep pace with the high rate of relative sea level rise.

Discussion

Model Uncertainties

Results from the SWAN wave model show a slight discrepancy between the modeled and observed swell wave power. This discrepancy can be attributed to error in the magnitude of swell energy at the offshore and lateral boundary as well as the approximation of wave diffraction in SWAN. The WISWAVE model results are used as the offshore wave boundary condition to incorporate the swell energy into our SWAN model domain. Tracy and Cialone (2004) compared the WISWAVE model results with measurement sites available in the Gulf of Mexico during 1995. Tracy and Cialone (2004) reported that the average RMSE for H_s and T_p was 0.30 m and 0.67 s, respectively. This is likely to cause some inconsistency in our swell wave power calculation; however, the conclusion that swell energy contributed to marsh edge erosion remains the same. Further, it is difficult to model swell energy in a shallow estuary, such as TB. Where diffraction effects are important, e.g., within a harbor or immediately behind barrier islands, phase-resolving models are typically required (Violante-Carvalho et al. 2009). The wave diffraction physics that are associated with swell wave energy transformation in partially sheltered shallow water are only an approximation in the SWAN model. To account for these physics, a phase-decoupled refraction-diffraction approximation based on the mild-slope equation was implemented into the model by Holthuijsen et al. (2003). It is important to include diffraction physics in the model to accurately represent the transformation of waves around the barrier islands that exist in TB. Holthuijsen et al. (2003) performed extensive test showing SWAN results improved when the approximation of wave diffraction was activated. However, error may still exist given the limitation of phase-averaged models in terms of modeling wave diffraction.

The uncertainty associated with the diffraction approximation in SWAN is minimized given our site conditions. The widths of the inlets between the barrier islands in TB are quite large, which reduces the effect of diffraction that takes place. Additionally, our model output locations are not immediately behind the barrier island where the effect of diffraction would be most significant. Holthuijsen et al. (2003) stated that the diffraction approximation in SWAN could become an issue when the distance from the obstruction to the coastline is small (less than a few wave lengths).

An interpolated lateral boundary condition is used in the SWAN model to represent wave conditions along the east and west boundaries of the domain. The purpose of constructing the model in this way is to reduce computation time while maintaining a high-resolution mesh. However, using the interpolated boundary conditions masks the effects of shoaling along the boundary and ignores refraction since a constant wave direction is assigned uniformly along the boundary. To ensure that there is no significant error introduced into the model, the accuracy of the imposed lateral boundary conditions is tested against a widened domain with no lateral boundary conditions. Using the same mesh resolution, the domain is extended wide enough to where the boundary effect has no influence on wave conditions inside Terrebonne Bay. This adjustment allows SWAN to calculate the 2D, nearshore wave transformations along the east and west boundaries of the original domain as the waves propagate into the bay. The model is rerun for the representative year of the study period, 2012, using the new widened domain. The average absolute percent change in P_b,comb between the two simulations is calculated to be only 2%. Therefore, it can be concluded that any error associated with the lateral boundary is negligible.

Potential Impact of Model Results

From the spatial analysis of wave power in TB, it is clear that the swell energy is the primary driver of marsh edge retreat in the southwest part of the bay as barrier islands continue to degrade in this region. It is important to understand the impacts that barrier islands have on reducing offshore wave energy and protecting the back barrier marshland. Therefore, a test case is carried out using the existing SWAN model to simulate the effects of barrier island degradation in TB. An average of the bulk wave power values for the years of 2004, 2005, 2006, 2008, and 2012 is obtained to determine the representative year for this simulation. Adjustments to the bathymetry data include completely removing the four offshore barrier islands by replacing their existing elevation with the equilibrium offshore elevation at each barrier island. All other boundary conditions remain unchanged. This scenario is run for the full duration of representative year, 2012.

As expected, a significant increase in wave power across TB is observed, primarily at the sites near the opening of the bay. Using the same classification for swell influenced locations defined earlier (P_{e, swell}/P_{e, sea} > 0.05), the annual averaged combined effective wave power, P_e,comb, is plotted in Fig. 12 to show this increase. Without barrier islands, the wave power is 1.5 times greater on average, and at some locations, can be up to 2.6 times greater. This result can be attributed to an increase in fetch for locally generated waves, which leads to larger wave growth. Additionally, this condition allows swell energy to easily propagate into the bay without being reduced or blocked by the barrier islands. In areas currently directly exposed to open waters, there is not expected to be a significant increase in P_e,comb. These locations can be observed along the 1:1 line in Fig. 12.

Fig. 12

Comparison of the annual averaged effective combined wave power, P_e,comb, with and without barrier islands at swell influenced locations of the study sites across Terrebonne Bay, LA

With an observed increase of wave energy in TB due to the removal of the barrier islands, it can be expected that the retreat rate at these locations will increase as well. The relationship established in Fig. 10a is used to quantify this amount. Normalizing the newly calculated P_e,comb values and inserting them into the equation, new retreat rates are obtained for the condition without barrier islands. The full magnitudes of these retreat rates are compared against the existing ones established for 2004–2012 (Fig. 13). The retreat rate is 2.1 times greater on average and a max of 9.9 times greater. Additionally, there are three locations where the retreat rates did not increase. The measured retreat rate at these locations is found to be abnormally high compared to the modeled wave power. Therefore, an increase in P_e,comb due to barrier island removal would have to be significant enough to overcome this abnormality when inserted into the equation.

Fig. 13

Comparison of the marsh edge retreat rates with and without barrier islands at swell influenced locations of the study sites across Terrebonne Bay, LA

Conclusions

This study applied the third-generation spectral wave model SWAN to simulate the nearshore wave climate in Terrebonne Bay, Louisiana. Modeling both the wind sea and swell energy in TB required multiple forms of analysis to be carried out. The SWAN model results were used to calculate wave power across TB. The wind sea and swell wave powers are computed separately, then added together to get the combined wave power in the bay. Daily-averaged wave power values are computed to show the seasonal variation of both sea and swell wave power magnitudes throughout the year at the study site. Good correlation is found between the modeled and observed daily-averaged wave powers for sea conditions and reasonable agreement for swell conditions. Though swell is important in the absence of other forcing, the local, wind sea waves clearly deliver more energy to the marsh edge at the study site in northwestern TB. The wave power in TB has strong temporal variability with March being the most energetic month in the cold front season.

The bulk and effective wave powers were calculated at 125 locations to show the spatial distributions around TB. The effective wave power, P_e, has been documented to be the value that causes marsh edge erosion (McLoughlin et al. 2015; Tonelli et al. 2010). Determining P_e is important when comparing to marsh retreat rates because it provides a more accurate representation of the wave thrust that acts on the marsh edge. A new range for effective wave power computation is proposed and shows, on average, a 25% reduction from the bulk values. Additionally, this analysis revealed that the magnitude of P_e varied significantly among the sites. However, the change in magnitude between years was considered to be negligible. Therefore, it can be concluded that modeling only 1 year of wave power is sufficient when used to predict retreat rates in a shallow water estuary. Additionally, the magnitude of swell energy in TB is established at all 125 locations as a percentage by calculating the ratio of P_{e, swell}/P_{e, sea} × 100. The highest percentages of swell wave power are found at the southwest entrance of TB. This may be a result of the dominant, southeasterly, offshore wave energy direction. Further, through the quantification of swell wave power around the bay, improvements can be made to wave hindcasts in long-term wave analyses for shoreline retreat prediction. The established percentages can be applied to wind sea wave energy results from parametric models, such as Young and Verhagen (1996), which do not account for swell energy. This is especially important at locations near the entrance of the bay where the swell energy is significant and local wind wave growth is limited by short fetches.

Spatial distributions of combined effective wave power, P_e,comb, are found to be in reasonable agreement with the established retreat rates in TB. Including the effects of swell energy improved the overall trend of the effective wave power and retreat rate comparison. As the barrier islands in this estuary continue to degrade, the contribution of swell energy to the overall energy budget at the marsh edge will increase, especially at locations partially exposed to the open waters of the Gulf of Mexico. However, this is only one part of fully understanding marsh edge erosion rates in TB. Extensive field campaigns should be carried out to determine marsh scarp height and soil strength properties in the bay to gain a robust outlook on the variability of marsh erosion in TB.

The methodology of quantifying and partitioning wind sea and swell effective wave power using the SWAN wave model integrated with field observations of wind and water levels as well as ocean-scale wave hindcast can be applied to other estuaries along the northern Gulf of Mexico and beyond. The proposed method to determine the effective wave power that impacts the marsh edge stability results in consistent annual effective wave power regardless of the number of hurricanes or frontal weather system passages in a particular year. The numerical experiment on the effect of barrier island degradation highlights the importance of restoring barrier islands that shelter wetlands from the impact of swell energy. The findings of our study will benefit other marsh edge erosion studies and wetland restoration in Mississippi River Delta and beyond.