Wave thrust values at point locations along the shorelines of Massachusetts and Rhode Island

This process step and all subsequent process steps were performed by the same person, Alfredo Aretxabaleta in Matlab (ver. 2016b), unless otherwise stated.
The latest European Centre for Medium Range Weather Forecast (ECMWF) Re-Analysis (ERA-5, https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5) model solution was extracted for all areas around the US east coast shoreline. The hourly data for the east-west and north-south wind components at 10 m height were downloaded from https://doi.org/10.24381/cds.adbb2d47 on 20 November 2019. Data for the period January 2000 to December 2018 was used for each location.

The east-west and north-south wind components are combined to produce wind speed and direction (degrees from True North) using a Matlab code that transforms and rotates vectors from the vector components to magnitude and direction.
The wind speeds and directions were binned to create a set of wind roses for each coastal location. Wind direction was binned in 12 intervals using 30-degree intervals: Bin 1: -15 (345) to 15 deg N; Bin 2: 15 to 45 deg N; Bin 3: 45 to 75 deg N; Bin 4: 75 to 105 deg N; Bin 5: 105 to 135 deg N; Bin 6: 135 to 165 deg N; Bin 7: 165 to 195 deg N; Bin 8: 195 to 225 deg N; Bin 9: 225 to 255 deg N; Bin 10: 255 to 285 deg N; Bin 11: 285 to 315 deg N; Bin 12: 315 to 345 deg N
and wind speeds were binned with the following ranges: Bin 1: 0-2 m/s; Bin 2: 2-4 m/s; Bin 3: 4-6 m/s; Bin 4: 6-8 m/s; Bin 5: 8-12 m/s; Bin 6 12-38 m/s

A matrix of model simulations (72 total simulations) of the ADCIRC/SWAN (version 53.04; https://adcirc.org/home/documentation/users-manual-v53/; Dietrich et al. [2011]) modeling system were conducted with different wind directions (12) and intensities (6). A constant wind for each direction centered on the mid-value (0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330 deg N) and with speeds centered at the mid-value of each bin (1, 3, 5, 7, 10, 25 m/s) were imposed for the entire model domain.
The unstructured grid is the one from the Hurricane Surge On-demand Forecast System (HSOFS, https://www.weather.gov/sti/coastalact_surgewg; Moghimi et al. [2020]). HSOFS is an unstructured finite element grid that extends westward to the 65 W longitude and resolves the entire east coast of the United States and the Gulf of Mexico. The grid contains 1.8 million grid points resulting in horizontal resolutions as low as 130m with most coastal grid elements being 300-400m in size. The grid extends overland to approximately the 10 m elevation (NAVD88).
The model was run in coupled mode, initialized from rest, and run until steady-state was achieved. The maximum water levels, significant wave heights, and peak periods were extracted from the end of each simulation.
Many papers describe the development and usage of the ADCIRC computational model, but basic details can be found in Luettich et al. (1992) and Dietrich et al. (2011).
References:
Dietrich, J.C., Zijlema, M., Westerink, J.J., Holthuijsen, L.H., Dawson, C., Luettich Jr, R.A., Jensen, R.E., Smith, J.M., Stelling, G.S. and Stone, G.W., 2011, Modeling hurricane waves and storm surge using integrally-coupled, scalable computations. Coastal Engineering, Vol 58(1), pp.45-65.
Luettich, R.A.; Westerink, J.J.; Scheffner, N.W. 1992, ADCIRC: An Advanced Three-Dimensional Circulation Model for Shelves, Coasts, and Estuaries; Report 1: Theory and Methodology of ADCIRC-2DDI and ADCIRC-3DL; Technical Report CERC-TR-DRP-92-6; U.S. Army Corps of Engineers, U.S. Department of the Army: Washington, DC, USA.
Moghimi, S., Van der Westhuysen, A., Abdolali, A., Myers, E., Vinogradov, S., Ma, Z., Liu, F., Mehra, A. and Kurkowski, N., 2020. Development of an ESMF based flexible coupling application of ADCIRC and WAVEWATCH III for high fidelity coastal inundation studies. Journal of Marine Science and Engineering, Vol. 8(5), p.308.

In order to estimate wave thrust, the orientation of the shoreline with regard to the incident waves causing thrust needed to be determined. We calculated shoreline orientation from the unstructured grid by following these steps:
a) We extracted shoreline segments (zero meter mean sea level contour) from grid using a Matlab code.
b) We calculated the angle with respect to True North of each shoreline segment.
c) We determined the wet (ocean) side of each segment by comparing the bathymetric/topographic elevation of the segment (0 m) with that of a point perpendicular to the segment and 100m from it. If the perpendicular point is deeper, then the ocean is to that side.
Depth_100m_perpen = depth(xmid+abs((y2-y1)/100), ymid+abs((x2-x1)/100)), where xmid, ymid are the coordinates of the center point of the segment, x1, y1 are the coordinates of one of the endpoints of the segment, and x2,y2 are the coordinates of the other endpoint of the segment. Depth_100m_perpen is the depth at a point perpendicular to the segment and 100m from it.
d) Finally, we determine the angle between the shore-normal direction and the incoming wave.

The climatological wave thrust was calculated as the weighted average of the wave conditions from each direction that was exposed to the waves coming from the ocean side.
The weights for each location were extracted from the wind distribution of directions and magnitudes from the climatological wind roses. Based on the shoreline orientation, the incident directions coming from land were assigned a weight of zero. The rest of directions were adjusted so that maximum weights correspond to the direction perpendicular to shore. The wave thrust estimation was conducted following the equations:
Based on Tonelli et al. (2010), first, the above mean sea level component (accounting for the hydrostatic pressure from wind waves) is defined as: WTasl = 0.5*ρ*g*Hwave^2, where ρ is the density of water, g is the acceleration due to gravity and Hwave is the significant wave height
Second, the dynamic pressure of wind waves can be calculated as WTbsl = ρ*g*Kp*Hwave, where Kp is the pressure-response factor accounting for the dynamic component due to water particle acceleration and can be calculated as Kp = cosh(k*(h+z)) / cos(k*h), where h is the marsh elevation with respect to mean sea level, z is the water level, and k is the wave number.
As wave thrust from dynamic pressure depends on water depth, we followed the approach of Tonelli et al. (2010) with a reduction in wave thrust proportional to the water depth on the marsh platform. Next, the thrust from the two components can be added to give the total thrust due to wave attack. WTtotal = WTasl + WTbsl
Based on the direction of the waves and grid orientation, the fraction of total thrust that is normal to the marsh edge is calculated
The model solutions provided significant wave height, peak wave period, and wave direction. As the below sea level expression requires wave number, the dispersion relationship was used to calculate wave number from wave period and water depth.
Reference: Tonelli, M., Fagherazzi, S., and Petti, M., 2010, Modeling wave impact on salt marsh boundaries: Journal of Geophysical Research, Vol. 115, C09028, https://doi.org/10.1029/2009JC006026.

This processing step was performed by Zafer Defne in ArcMap (ver. 10.7.1) using tools from ArcToolbox. For complex operations, names of specific tools used are given in CAPITAL letters (any critical parameters used are given in parentheses, separated by a semicolon, immediately after the tool name). The input (GULF_ATLANTIC_ESI) and output file names are provided in [square brackets] when necessary. Units for length and area calculations are meters (m) and square meters (m2) unless otherwise stated.
In this step we generate points along the high resolution (1/24,000 scale) shoreline. We generate a 1-km spaced and a 100-m spaced dataset.
Download the 2017 National Environmental Sensitivity Index Shoreline (Contiguous U.S. shoreline) dataset from NOAA Office of Response and Restoration's Environmental Sensitivity Index (ESI) maps inventory. See https://response.restoration.noaa.gov/esi_download.
Crop the [GULF_ATLANTIC_ESI] line features in the GULF_ATLANTIC_ESI geodatabase file from the Connecticutt-Rhode Island border to the Massachusetts-New Hampshire border.
GENERATE POINTS ALONG LINES (Point Placement =Distance; Distance=0.009 Decimal degrees) along the cropped shoreline approximately with 1-km spacing to generate the point features [Massachusetts_shoreline_wave_thrust_HSOFS_dir_1km.shp].
GENERATE POINTS ALONG LINES (Point Placement =Distance; Distance=0.0009 Decimal degrees) along the cropped shoreline approximately with 100-m spacing to generate the point features [Massachusetts_shoreline_wave_thrust_HSOFS_dir_100m.shp].
The distance along the shoreline is calculated straight out from the source data. Therefore, the spacing between adjacent points in either the 100m or 1km datasets might be much smaller than 100m or 1km, as the spacing is from the original data and not between points in the extracted points.

In the final step, the wave thrust calculated on the ADCIRC grid (HSOFS) is linearly interpolated onto the shoreline positions for the States of Massachusetts and Rhode Island with either the 100-meter spacing or the 1-km spacing between shoreline points. The interpolation was conducted using the scatteredInterpolant function in Matlab with the linear interpolation option. In Matlab v2016b, the function dlmwrite.m was used to export the data to CSV format. The shapefile was created in ArcGIS from the CSV files.