Modeling passive marine and groundwater flooding at high tide

Passive flood modeling is described in the main text, under the section: Approach and Models, Passive flooding at high tide. For validation, modeled flooded areas were compared against known elevations during present day conditions.

In passive modeling, we assume that the amount of SLR is known. Thus, we estimate the error in the baseline scenario \({\sigma }_{base}\) as the increase in exposed land area that corresponds to ~10 cm of SLR (10 cm is the vertical error in the DEM described below):

$${\sigma }_{base}=Are{a}_{SLR=0.18m}-Are{a}_{base,(SLR=0.08m)}$$ (2)

For higher SLR values, we assume that the error increases proportionally to projected exposure area; and, the constant of proportionality is the fractional uncertainty (error in area divided by estimated area) of the baseline area. This helps account for larger errors that correspond to larger increases in exposure area when, at higher SLR scenarios, ground slopes are generally milder than the baseline scenario. At milder slopes, the vertical error in the DEM leads to larger errors in exposed area projections. The resulting error for projected exposure under future SLR is:

$${\sigma }_{future}=Are{a}_{future}\cdot \frac{{\sigma }_{base}}{Are{a}_{base}}$$ (3)

Digital Elevation Model for Passive Flood Modeling

The DEM for passive flood modeling is a digital surface that represents the local terrain. It consists mainly of a 1 m resolution bare-earth DEM that was derived from light detection and ranging (LiDAR) data collected in 2013 by the U.S. Army Corps of Engineers (USACE) (http://shoals.sam.usace.army.mil/). This DEM extends seaward to depths of approximately 35 m, while the landward coverage extent is typically between 400 and 1000 m inland from the shore. The 2013 DEM was augmented with the 3 m resolution NOAA SLR Viewer DEM36 to fill in low elevations that extend far inland. Only the terrain portion of the DEM (not the bathymetric portion) was used for passive flood modeling.

Due to the collection and processing methods used on the 2013 USACE LiDAR data, many streams and wetlands are poorly classified or characterized. Where possible, satellite imagery of a similar date was used in conjunction with a hillshade rendering of the DEM to identify misclassified or misinterpreted areas and apply a corrected elevation to the depicted stream or wetland area that approximates nearby valid elevations. This “hydro-flattening” is intended to better depict low-lying wetland areas and streams in the modeling, and to ensure waterway connectivity to the ocean in the presence of bridges and similar structures.

Modeling wave inundation

The non-hydrostatic mode of XBeach was used in this study to capture short wave contributions to water levels, which are often required to accurately model wave overtopping56. For validation in Hawai‘i, we modeled wave heights and water levels using XBeach, and compared them against measurements from an array of bottom-mounted pressure sensors located on O‘ahu.

Projections of water depths and velocities were estimated, for each SLR scenario, as follows. We run a 1-hour (model time) XBeach simulation using maximum annually recurring wave conditions along shore-normal profiles spaced 20 m apart. Then, we define the extreme water level as the average of the five highest water elevations at each grid point along each profile, and identify the corresponding velocity. Water depth is calculated as the difference between the annual extreme water level surface and the ground elevation. Depths and velocities are then interpolated to create 2D surfaces using ArcGIS, and resampled at 5 m grid spacing. Modeled parameters (depth, velocity) at grid locations with velocity <0.1 m/s or depth <0.1 m are omitted to ignore static, inland water that does not result from waves, and to remove thin veneers of water excursion. The error in land area exposed to wave inundation is determined using the same method described in modeling passive flooding, except the \({\sigma }_{direct}\) value is calculated using a 17 cm vertical error, rather than 10 cm, to account for additional uncertainty in the wave model.

Data and boundary conditions

We augmented the 1 m resolution 2013 USACE DEM with lower resolution (10–50 m) bathymetric data collected by the Hawai‘i Mapping Research Group for water depths that exceed the maximum LiDAR depth range (roughly >35 m depth). However, the lower resolution data is rarely used since most modeling is conducted over shallower water depths. To fill bathymetry gaps within nearshore areas, we use a 10 m DEM that was created from all 1999–2003 USACE bathymetric and topographic surveys available from NOAA’s Digital Coast36. Data gaps in shallower water (<10 m) are left unfilled because we find that these gaps only occur in small patches whose spatial footprints are much smaller than the resolution of the 1999–2003 10 m DEM; and using the lower resolution DEM values to fill these gaps results in obvious elevation errors near the shoreline. Instead, we linearly interpolate across any small gaps in elevation after we have extracted the 1D profiles from the DEM.

Mean water levels are determined by increasing the water level above MHHW by the desired amount of SLR. To obtain wave information at the local level, we use the wave forecast available from PacIOOS57, which has been providing a 5-day hourly forecast that is calibrated using local wave buoys, for Hawaiian regions since 2009. This forecast is produced by using the WaveWatchIII and SWAN wave models to dynamically downscale atmospheric forcing from the NCEP Global Forecast System weather model57. From this forecast, we use the hourly time series of significant wave height, peak wave period, and dominant direction, for the years 2009–2015.

We then follow the method presented by Vitousek and Fletcher39 and extract peak wave heights from the time series of significant wave height, and fit those to a generalized extreme value (GEV) distribution, from which the maximum annually occurring significant wave height and the associated peak period are found. To determine the dominant direction, we perform the GEV fit on subsets of wave series that are limited to incoming wave directions within 30-degree bins, and let the bin ranges shift by 15 degree increments (a maximum of 24 bins; there will be overlap between bins). The direction from which waves will produce the highest run-up is typically shore-normal, and has the largest wave height and period. On the rare occasion that these three conditions do not align, we use expert knowledge and generally select the direction corresponding to longer period waves. See Supplementary Tables S2–S4 (excel spreadsheets) for wave parameters used.

We repeat this procedure on wave data at locations spaced 1.5–2.5 km apart, in roughly 25–35 m depth around each island. The annually-recurring maximum wave parameters are then interpolated along the offshore boundary to match profile locations. Since projections of future wave conditions due to climate change do not indicate significant changes in wave characteristics near Hawai‘i, we use the same wave parameters for all future SLR scenarios58. From the wave parameters, we use XBeach routines to generate the 1-hour time series of incoming waves.

Friction due to bottom roughness affects wave flow velocity. To determine bottom roughness, we use two mapping products (GIS polygon layers) that classify the geology and structural type of land and sea topography. For land areas, we use the USGS geology classification map for Hawai‘i59; for submerged areas, we use benthic structural classification maps (e.g. sand, spur and groove reef, pavement reef) produced as part of the Mapping of Benthic Habitats for the Main Eight Hawaiian Islands project for NOAA’s National Ocean Service Center for Coastal Monitoring and Assessment60. We then assign friction values to each geologic type and structural classification (see Supplemental Tables S5–S6 (excel spreadsheets) for friction values). Then, via a lookup table, we assign these values to each polygon. The polygon layers are merged and gridded to match the DEM grid size. Because rough reef platforms can episodically become exposed and/or covered by sand, and because the USGS geology maps are lower-resolution than our DEMs, we edited some areas of the USGS and NOAA maps to match conditions at the time of the DEM data collection. When needed, we hand-digitized replacement polygons that more accurately represent structural classification by looking at the existing DEM and most-recent aerial photo, then merging this new layer with the existing one.

Modeling chronic coastal erosion

Following the probabilistic method described in Anderson et al.17 (and introduced in equation (1) above) we produce a probability density function (pdf) for the projected vegetation line y veg (t) by numerically combining the pdfs for each term on the right hand side of equation (1). We then determine the erosion hazard line as the 80th percentile,\({y}_{veg,80 \% }={{F}_{veg}}^{-1}(80/100)\) where \({F}_{veg}={\int }_{-\infty }^{80/100}pd{f}_{veg}\) is the cumulative distribution function of the pdf for the projected vegetation line. The 80th percentile is used to align the erosion hazard estimate with the other models (passive flood and wave) that both use the upper limit (83rd percentile) of the IPCC probability range for future sea level.

Calculations are performed at shore-normal transects, spaced 20 m apart along the shore. The erosion hazard line at each transect (the 80th percentile of the pdf for the projected vegetation line) is then projected back to map coordinates and exported as a GIS line shapefile. These calculations are done using Matlab. Some edits are manually performed in ArcGIS to smooth areas where transects overlap along coasts with high curvature. Erosion hazard areas are created in ArcGIS by connecting the erosion hazard line with the existing shoreline (MSL elevation contour).

The error in new land area exposed to erosion for future sea level S f and time t is determined as \({\sigma }_{area({S}_{f},t)}=\sum _{i}{\rm{\Delta }}x\cdot {\sigma }_{i,({S}_{f},t)}\), where \({\sigma }_{i,({S}_{f},t)}\) is the standard deviation for the ith projected shoreline position, and Δx is the distance between transects (20 m).

Data used in erosion modeling

We use historical shorelines and vegetation lines that were digitized as part of the USGS National Assessment of Shoreline Change48 and augmented with newer digitized shorelines for the north Maui and west Maui regions49. The positional error in each digitized shoreline was estimated as the root mean square of up to seven sources of positional error. The reader is referred to Fletcher et al.48 for a detailed description of the historical shorelines.

At some beaches, past human activity, such as sand mining in the early 20th century, can influence the long-term shoreline change rate and violate the assumption that shoreline change remains fairly constant over the roughly 80 years of available data. To capture ongoing behavior of the beach, we calculate shoreline change rates only over the period of time following any significant beach alternation, such as sand mining. In the case where a beach is completely lost to erosion and a seawall has been built, we calculate the shoreline change rate up to and including the first shoreline where the wall was erected. Since we cannot predict if a wall will be built or removed, we use the preexisting rate (before the wall was built) to provide an estimate of exposure to erosion, in the beach’s natural state.

A straight line is fit to the historical shoreline data using weighted least squares regression61,62. Rates are then smoothed in the alongshore direction using a weighted [1 3 5 3 1] moving average to reduce spurious variations between adjacent beach locations.

To determine beach slope, we use the collection of biannual cross-shore beach surveys (1994–1999) on beaches of O‘ahu and Maui63, as well as additional surveys at 35 locations on O‘ahu and 27 locations on Kaua‘i over the period 2006–200817. Following Anderson et al.17, we extract, for each repeated beach survey, the slope between two morphologic features: the beach toe64, which is located at base of the foreshore, and; either the point at which the repeated surveys converge (depth of closure) or the intersection of the sandy profile with the reef platform65. The final beach slope, tanβ, is determined as the average of the individual slopes. We use the standard deviation of the slopes as our estimate for the error associated with tanβ. If more than one profile location exists along a study area, slopes are interpolated using cubic splines.