Uncertain seas: probabilistic modeling of future coastal flood zones

ABSTRACT Future sea-level rise will likely expand the inland extent of storm surge inundation and, in turn, increase the vulnerability of the people, properties and economies of coastal communities. Modeling future storm surge inundation enhanced by sea-level rise uses numerous data sources with inherent uncertainties. There is uncertainty in (1) hydrodynamic storm surge models, (2) future sea-level rise projections, and (3) topographic digital elevation models representing the height of the coastal land surface. This study implemented a Monte Carlo approach to incorporate the uncertainties of these data sources and model the future 1% flood zone extent in the Tottenville neighborhood of New York City (NYC) in a probabilistic, geographical information science (GIS) framework. Generated spatiotemporal statistical products indicate a range of possible future flood zone extents that results from the uncertainties of the data sources and from the terrain itself. Small changes in the modeled land and water heights within the estimated uncertainties of the data sources results in larger uncertainty in the future flood zone extent in low-lying areas with smaller terrain slope. An interactive web map, UncertainSeas.com, visualizes these statistical products and can inform coastal management policies to reduce the vulnerability of Tottenville, NYC to future coastal inundation.


Introduction
The rate of global sea-level rise will likely increase due to increases in ocean temperature and increases in land-ice melt (Parris et al. 2012, Church 2013. Coastal properties are already prone to flooding and coastal flooding is expected to be exacerbated by future sea-level rise (Wong et al. 2014, Sweet et al. 2017, Fleming 2018. In 2012, Superstorm Sandy resulted in approximately $19 billion in damage and 48 lives lost in New York City (NYC, Blake et al. 2013). There were two deaths in the Tottenville neighborhood of NYC during Superstorm Sandy, and this neighborhood is the case study for modeling future flood zones from storm surge inundation enhanced by sea-level rise.
Sea-level rise in NYC has already increased the number and magnitude of coastal flood events (Sweet et al. 2013, Talke et al. 2014, and coastal flood zones will likely CONTACT Christopher J. Amante Christopher.Amante@colorado.edu; Christopher.Amante@noaa.gov expand in NYC in the future (Horton et al. 2015). Sea-level rise is a relatively slow process, however, it will likely continue to increase the frequency, magnitude, and duration of storm surge inundation (Parris et al. 2012). Future storm surge inundation enhanced by sea-level rise will, in turn, likely increase the vulnerability of the people, properties and economy of NYC. Modeling future storm surge inundation enhanced by sea-level rise uses numerous data sources with inherent uncertainties. There is uncertainty in (1) hydrodynamic storm surge models, (2) future sea-level rise projections, and (3) topographic digital elevation models (DEMs) representing the height of the coastal land surface above established vertical datums. This study implements a Monte Carlo approach to incorporate the uncertainties of these data sources and model the future 1% flood zone extent in the Tottenville neighborhood of NYC in a probabilistic, geographical information science (GIS) framework.

Literature review
The following sections provide background information on storm surge, sea-level rise, and DEMs, and their inherent uncertainties, as well as previous research on modeling future storm surge inundation enhanced by sea-level rise.

Storm surge
Storm surge is the build-up of water onto coastal land from wind shear stress associated with intense low-pressure weather systems, such as tropical and extratropical cyclones (Murty et al. 1986). Present-day coastal flood zones from storm surge inundation are typically determined by hydrodynamic models including the Sea, Lake, and Overland Surges from Hurricanes (SLOSH, Jelesnianski et al. 1992) and the ADvanced CIRCulation (ADCIRC, Luettich et al. 1992) models.
Many studies use hydrodynamic models, such as SLOSH and ADCIRC, directly to delineate potentially inundated areas (e.g., McInnes et al. 2003, Kleinosky et al. 2007, Frazier et al. 2010, Shepard et al. 2012, Atkinson et al. 2013, Ding et al. 2013, Maloney and Preston 2014. Other coastal inundation studies use products derived from storm surge model outputs, such as the Federal Emergency Management Agency (FEMA) Flood Insurance Rate Maps (FIRMs, e.g., Patrick et al. 2015). FEMA FIRMs are typically derived from outputs from ADCIRC and coupled wave models, such as the Simulating WAves Nearshore (SWAN) model (Algeo and Mahoney 2011, FEMA 2014, Kress et al. 2016. FIRMs determine flood insurance rates based on the 1% annual chance of inundation and identify where flood insurance is required as a condition of a federally-backed mortgage (Burby 2001).
Sources of uncertainty in storm surge models and, consequently, in derived products such as FEMA FIRMs, originate from the input parameters, including wind speed and direction, bathymetry, topography, friction coefficients, and boundary conditions (Atkinson et al. 2013). There is also uncertainty in the historical meteorological data on extreme events and the choice of statistical functions that represent their occurrence (McInnes et al. 2003).

Sea-level rise
Sea-level rise provides an elevated water base for storm surge to build on, and future sea-level rise will likely continue to increase the frequency, magnitude, and duration of storm surge inundation (Parris et al. 2012). Coastal flooding occurs at the land-water interface, and, therefore, local information on the relative vertical movement between the land and water surface is required to model future flood zones. Local sea-level projections can diverge from the global mean due to differences in ocean temperature, salinity, currents, and land elevation changes due to glacial isostatic adjustment, human extraction of ground water, and tectonic processes (Nerem and Mitchum 2001, Cazenave and Nerem 2004, Lombard et al. 2005, Milne et al. 2009, Church 2013. Relative sea-level rise in NYC is expected to exceed the global average primarily due to local land subsidence and to increases in regional sea-level, due in part to projected weakening of the Gulf Stream current (Yin et al. 2009, Kenigson and Han 2014, Horton et al. 2015).
Other regional factors that can affect future sea-level in NYC include climate modes, such as the North Atlantic Oscillation (NAO, Barnston and Livezey 1987, Hurrell 1995, Kopp 2013, Han et al. 2017) and the Atlantic Multidecadal Oscillation (AMO, Trenberth and Shea 2006, Wang and Zhang 2013, Han et al. 2017. Han et al. (2017) identified limitations of understanding regional sea-level change, including representation of these climate modes, interactions among climate modes, effects of anthropogenic forcing on the modes, effects of ocean internal variability, and limited observational records. These limitations must be addressed to fully understand the climate modes' impact on regional sea-level, and to reduce uncertainty in future sea-level rise projections (Han et al. 2017).

Digital elevation models (DEMs)
The topographic elevation above sea-level is a primary factor that determines the inland extent of coastal flooding. Accordingly, one of the most common methods for modeling coastal inundation is 'bathtub modeling.' Bathtub modeling delineates flood zones as areas where coastal land elevation values are below a modeled water height that represents a sea-level rise projection and/or storm surge level (e.g., Poulter and Halpin 2008, Gesch 2009, Leon et al. 2014, Schmid et al. 2014, West et al. 2018. Different representations of the topographic elevations, consequently, can result in large differences in the modeled flood zones (Coveney and Fotheringham 2011, Zhang 2011, West et al. 2018, especially for low-lying coastal areas (Van de Sande et al. 2012). Estimations of potential DEM errors, i.e., DEM uncertainty, should, therefore, be incorporated in flood models to produce a range of possible flood zones.
DEMs can be generated from various data sources and methods (Nelson et al. 2009, Wilson 2012, Maune and Nayegandhi 2018, including traditional ground-surveys with theodolites (Nelson et al. 2009), satellite-derived stereo photogrammetry (Shean et al. 2016, Almeida et al. 2019, structure-from-motion (Mancini et al. 2013), radio detection and ranging (RADAR) interferometry (Rabus et al. 2003), and light detection and ranging (LIDAR, Maune and Nayegandhi 2018). These various data sources and methods result in DEMs with different spatial resolutions and vertical accuracies (Nelson et al. 2009). Furthermore, coastal DEMs are often generated by integrating numerous elevation data sources of disparate age, quality and measurement density (Eakins and Taylor 2010, Eakins and Grothe 2014, Amante 2018b, resulting in spatially-varying vertical uncertainty (Amante 2018b).
Modeling future storm surge inundation enhanced by sea-level rise There are two general methodologies to model future storm surge inundation enhanced by sea-level rise: dynamic (also known as hydrodynamic or numerical simulation) and static (also known as bathtub) modeling. Dynamic modeling evaluates sea-level rise projections, and changes storm surge and wave model input variables, such as ocean depths and bottom friction coefficients, before modeling storm surge (Atkinson et al. 2013, Orton et al. 2015.
Conversely, static modeling evaluates sea-level rise projections after obtaining the output from present-day storm surge and wave models, or indirectly, such as with the present-day FEMA FIRM 1% flood zone, in a GIS framework , Leon et al. 2014, Patrick et al. 2015. The static, linear addition by expansion method adds a sealevel rise increment to the present-day modeled storm surge heights. The future storm surge inundation area is then delineated as areas where the cumulative water levels are greater than the DEM values (McInnes et al. 2013. The linear addition by expansion method does not incorporate physical forces on the water movement considered in dynamic modeling, such as differences in bottom friction due to changing water depths and landscape characteristics. However, the linear addition by expansion method more closely mimics the results of dynamic methods compared to other static methods . Inland local depressions that have elevations below a projected water level and are not connected to the ocean are incorrectly inundated in GIS-based, static methods that do not implement water connectivity algorithms (e.g., Titus and Richman 2001). In these cases, terrain barriers exist between the ocean and the low-lying areas that would prevent inundation (Li et al. 2009). In studies that implement water connectivity algorithms, areas are considered inundated only where their elevations are less than the modeled water levels and are also adjacent to the current or modeled future ocean area (e.g., Poulter and Halpin 2008, Gesch 2009, Li et al. 2009. Adjacency is typically defined by either 4 or 8 neighbors in a 3 × 3 kernel. The 4-neighbor case only considers the 4 cardinal directions (i.e., N, E, S, W), and the 8-neighbor case considers all adjacent cells (i.e., N, NE, E, SE, S, SW, W, NW).
Dynamic methods of modeling future storm surge inundation enhanced by sea-level rise are more complex than GIS-based, static methods, which results in larger computational expense. Consequently, dynamic methods typically cannot simulate numerous sea-level rise projections, storm surge scenarios, and DEM realizations to incorporate the uncertainties of these data sources (Orton et al. 2015). Zhang et al. (2013) found that dynamic methods require approximately 30 to 60 times more computation time than the static linear addition by expansion method to model future storm surge inundation enhanced by a given sea-level rise projection. Therefore, computational expense is greatly reduced if static methods can approximate dynamic methods of modeling future coastal flood zones ).
GIS-based, static methods assume storm surge inundation dynamics are the same with current and future sea-levels, and the enhanced storm surge heights are simply the linear addition of the present-day storm surge heights and the future sea-level rise. However, the relationship between storm surge heights and sea-level rise can be nonlinear in shallow water due to changes in bottom friction and shoreline configuration (Lowe et al. 2001, Lowe and Gregory 2005, Loder et al. 2009, Ding et al. 2013).
There are a few, notable studies that directly compare the results of dynamic and GISbased, static modeling of future coastal inundation (e.g., Zhang et al. 2013, Orton et al. 2015, Patrick et al. 2015, Seenath et al. 2016. Zhang et al. (2013) determined that inundation areas and peak storm surge heights generated by the static, linear addition by expansion method differed from the dynamic method by only 7% and 4% on average, respectively. Seenath et al. (2016) also determined that a GIS-based, static method produced an inundation area that was within 5% of dynamic methods.
Two studies modeled future storm surge inundation enhanced by sea-level rise in NYC to compare GIS-based, static methods (Patrick et al. 2015) and dynamic methods (Orton et al. 2015). The static and dynamic methods resulted in similar storm surge heights for most locations (usually within approximately ± 0.15-m) and resulted in similar future flood zone boundaries (Orton et al. 2015). The main purpose of these NYC studies was to compare static and dynamic methods, however, both studies acknowledged limitations due to the lack of incorporation of input data uncertainties into their analyses. Specifically, Patrick et al. (2015) stated estimates of uncertainty associated with the elevation, sea-level rise, and FEMA flood heights should be used to determine the degree of confidence in the flood model outputs. Lin et al. (2012) also determined that static methods were an excellent approximation of dynamic methods in NYC. The bottom friction is expected to remain relatively constant under future sea-level rise in urban areas, where the existing land cover is not dominated by vegetation, and, consequently, the effect of changes to bottom friction is expected to be minor compared to the effect of the water level increase from sea-level rise (Atkinson et al. 2013). These previous studies indicate that GIS-based, static methods can approximate dynamic modeling of future flood zones in NYC. Importantly, static methods enable the incorporation of input data uncertainties to determine the likelihood of the future flood zone extent in a probabilistic framework, which would be computationally impractical with dynamic methods.

Probabilistic future flood models
Deterministic flood models that do not consider input data uncertainties can result in inaccurate and misleading future coastal flood zones, and potentially inappropriate coastal management policies (West et al. 2018). Probabilistic flood models that consider these data uncertainties can reduce the conveyed over-confidence of a single, deterministic flood model, demonstrate the effect of data uncertainties on the flood model, and produce a range of possible alternative flood model outcomes (Cowell et al. 2006, Hare et al. 2011, Wallentin and Car 2013, Winter et al. 2018. Flood models should, therefore, incorporate the uncertainties in storm surge models, sea-level rise projections, and DEMs in a probabilistic framework to determine future coastal flood zones. Probabilistic flood models can then inform coastal management policies by indicating the likelihood of the future flood zone extent.
Several previous studies modeled future flood zones by incorporating the uncertainty in one or two of the major sources of input data, i.e., present-day storm surge, sea-level rise projections, and DEM (e.g., Gesch 2009, Li et al. 2009, Neumann et al. 2010, Zhang 2011, Strauss et al. 2012, Albert et al. 2013, Cooper and Chen 2013, Leon et al. 2014, Schmid et al. 2014, Kane et al. 2015, West et al. 2018. There is limited research on modeling future storm surge inundation enhanced by sea-level rise that incorporates these three major sources of uncertainty in a probabilistic framework. Furthermore, almost all previous studies focused on the spatial uncertainty of inundation extent that results from the uncertainties of these data sources. One notable exception, Kane et al. (2015), determined the time frame when sea-level rise could result in a rapid increase in areas prone to flooding. Leon et al. (2014) incorporated DEM uncertainty in modeling the combined effect of a uniform 1% storm surge height of 2.9-m and a 1-m sea-level rise projection. Leon et al. (2014) improved the incorporation of DEM uncertainty in coastal flood models by spatially distributing DEM errors based on land cover and terrain variables. Specifically, Leon et al. (2014) quantified vertical errors as differences between the DEM and more accurate real-time kinematic (RTK) GPS elevation measurements. The errors were then correlated with land cover and terrain variables and 1,000 vertical error surfaces were generated using sequential Gaussian simulations and regression-kriging in a Monte Carlo, probabilistic approach. A generated error surface was added to the DEM in each simulation and the probability of future inundation was calculated as the proportion of times a DEM cell was inundated from 1,000 simulations. Leon et al. (2014), however, only modeled one sea-level rise projection, and, therefore, did not incorporate the uncertainty of future sea-level rise. Furthermore, Leon et al. (2014) modeled a uniform storm surge height of 2.9-m for the entire study area. Storm surge heights typically vary along a coastline due to the offshore bathymetry and nearshore topography (Jelesnianski et al. 1992, McInnes et al. 2003. The methods in Leon et al. (2014), and their limitations, informed the GIS framework for probabilistic modeling of future coastal flood zones in this study. Specifically, this study implemented a similar probabilistic approach as Leon et al. (2014) that incorporated DEM uncertainty, and also incorporated the uncertainties in the present-day 1% flood zone from storm surge and future sea-level rise projections.

Methodology
This study implemented a Monte Carlo approach to model the future 1% flood zone extent in Tottenville, NYC in a probabilistic, GIS framework. The future flood model incorporated the uncertainties in a topographic DEM, present-day 1% flood zone from storm surge, and future sea-level rise projections. A 500-member, future flood model ensemble was generated from random combinations of input data realizations within estimated uncertainty bounds. Open source software, including MB-System (Version 5.4.2220, Caress and Chayes 1996) and Python Version 2.7, was used to develop the future flood model and to derive spatiotemporal statistical products from the Monte Carlo ensemble. The derived statistical products can inform coastal management policies by indicating the likelihood of the future 1% flood zone extent in Tottenville, NYC.

Study area
Tottenville is the southernmost neighborhood of Staten Island, NYC ( Figure 1). The landuse in Tottenville is primarily residential, and the neighborhood is surrounded by water in three cardinal directions. The Arthur Kill is located to the west, and the Raritan Bay is located to the south and to the east of the neighborhood. Tottenville has a maximum elevation of~30-m above the North American Vertical Datum 1988 (NAVD 88). Most of the neighborhood is currently protected from coastal flooding at these higher elevations, however, additional areas will likely become prone to future storm surge inundation enhanced by sea-level rise and be in the future 1% flood zone.

Digital elevation model (DEM)
Woolpert, Inc. collected topographic LIDAR between March 21-21 April 2014 for the U.S. Geological Survey (USGS) Coastal Marine and Geology Program (CMGP) Post-Sandy project (Woolpert 2014). The data set was collected to achieve a maximum nominal post spacing of LIDAR returns of 0.7-m, and 0.42-m horizontal accuracy at 95% confidence (Woolpert 2014). The overall vertical accuracy for bare-earth LIDAR returns was 5.3-cm root mean square error (RMSE z , Woolpert 2014).
A bare-earth, topographic DEM ( Figure 2, Left Panel) was generated from the LIDAR bare-earth returns with spline interpolation from the open-source software tool, MB-System's 'mbgrid.' The generated DEM was referenced horizontally to the World Geodetic System 1984 (WGS 84) and vertically to NAVD 88. The DEM spatial resolution, 1/9 th arc-second (~3-m), followed the highest-resolution of the framework collaboratively developed by NOAA National Centers for Environmental Information (NCEI) and USGS (Eakins et al. 2015). Refer to Eakins and Taylor (2010) and Eakins and Grothe (2014) for additional details on NOAA NCEI coastal DEM generation and associated challenges.
An uncertainty surface that estimated potential DEM cell-level vertical errors at one standard deviation (Figure 2, Right Panel) was generated with the methods from Amante (2018b). The uncertainty surface represented contributions from the (1) measurement uncertainty provided in the LIDAR metadata (i.e., 5.3-cm RMSE z ) with an assumed nonbiased, normal distribution, (2) number of measurements per DEM cell, (3) subcell measurement variance, and (4) interpolation uncertainty in cells unconstrained by measurements. In DEM cells with LIDAR bare-earth returns, the cell-level vertical uncertainty was calculated from a pooled standard deviation that incorporated the measurement uncertainty of 5.3-cm RMSE z, and subcell measurement variance from the mean elevation of the DEM cell. The pooled standard deviation was divided by the square root of the number of measurements in the DEM cell to calculate the standard deviation of the mean, i.e., the standard error, to represent the DEM cell-level vertical uncertainty. There was additional uncertainty in interpolated DEM cells with no LIDAR bare-earth returns (e.g., water bodies and within building footprints), with larger distances from elevation measurements resulting in larger interpolation uncertainty Eakins 2016, Amante 2018b). Refer to Amante (2018b) for additional details on the methods to estimate coastal DEM uncertainty at the individual cell-level.
Amante (2018b) created a DEM uncertainty surface that was spatially correlated with land cover and terrain slope. To confirm these results from Amante (2018b), this study . The largest uncertainties, and consequently, the largest differences between the DEM realizations in the future flood model, are located where there are large interpolation distances due to no LIDAR bare-earth returns, including inland bodies of water and within building footprints. calculated the average cell-level uncertainty for land cover classes from the 2010 NYC Department of Parks and Recreation land cover raster (NYC Department of Parks and Recreation 2010). Furthermore, this study quantified the effect of both the land cover and terrain slope derived from the DEM (Figure 3) on the magnitude of the cell-level uncertainty with Spearman's rank correlation coefficients (Spearman 1904).
The two land cover classes with the largest, average cell-level uncertainty were the 'water' and 'building' classes due to few LIDAR bare-earth returns and interpolated DEM values in cells unconstrained by elevation measurements (Table 1). The two land cover classes with the smallest, average cell-level uncertainty were the 'road/railroad' and 'bare Figure 3. The terrain slope derived from the DEM. The slope varies within the study area, with larger terrain slope along the western coastline and smaller terrain slope along the southern coastline of the neighborhood. soil' classes due to effective laser pulse penetration in open terrain, and consequently, many LIDAR bare-earth returns per DEM cell (Table 1). The Spearman's rank correlation coefficients indicated a negative correlation between the magnitude of the DEM uncertainty and the number of measurements for these four land cover classes (Table 1), i.e., more measurements per DEM cell resulted in smaller cell-level uncertainty, as expected per Amante (2018b). With the exception of the 'water' class, the terrain slope was positively correlated with the cell-level uncertainty (Table 1), as larger terrain slope resulted in larger subcell measurement variance, as expected per Amante (2018b). Topographic LIDAR systems typically do not produce bare-earth returns within water bodies. In such areas, DEM values were interpolated based on the elevations of bare-earth returns on the perimeters of these water bodies and produced relatively flat elevation surfaces. DEM uncertainty, therefore, was negatively correlated with terrain slope for the 'water' class because areas with the largest elevation uncertainty were in the center of water bodies due to large interpolation distances from the LIDAR bare-earth measurements, and these areas have small slope values. The DEM uncertainty was most positively correlated with terrain slope for the 'bare soil' class (Table 1), due to larger subcell measurement variance in steep, near-shore sand dunes.
This study generated DEM realizations for the future flood model using the DEM and derived cell-level uncertainty surface at one standard deviation. A normal, non-biased distribution was assumed, and the 'best-case' and 'worst-case' error realizations were created by multiplying the entire uncertainty surface at one standard deviation by 1.96 and −1.96 to represent the 95% confidence interval, respectively. The resulting error realizations were each added to the original DEM to create the 'best-case' and 'worstcase' DEM realizations, respectively. 498 additional, intermediate DEM realizations were then randomly generated using the derived uncertainty surface at one standard deviation and a mean of zero with the Numerical Python (NumPy) package function 'numpy. random.normal'. A 3 × 3 low-pass filter was applied to reduce the introduced noise and maintain the spatial autocorrelation of the terrain (Wechsler andKroll 2006, Wechsler 2007), and each resulting error realization was added to the original DEM. The 500 generated DEM realizations, with a normal distribution, represented the DEM uncertainty incorporated in the future flood model.

Present-day storm surge
Two versions of the present-day FEMA 1% flood zones, the current, accepted 2007 FIRM, and the 2013 Preliminary FIRM, represented the 'best-case' and 'worst-case' scenarios of the present-day storm surge and wave heights in the future flood model, respectively. On 17 October 2016, FEMA announced the administration of the mayor of NYC, Bill de Blasio, won its appeal of the 2013 Preliminary FIRM, and FEMA agreed to revise the NYC FIRM (FEMA 2016). The appeal cited two primary sources of bias in the storm surge and offshore wave models that resulted in more inland areas within the 1% flood zone than warranted: (1) insufficient extratropical storm model validation and (2) misrepresentation of tidal effects for extratropical storms (Zarrilli 2015).
498 additional, intermediate present-day 1% flood zone realizations were randomly generated between the 2007 and 2013 FIRM 1% flood zone extents. Intermediate, present-day total flood areas were randomly sampled between the total areas in the 1% flood zone for the 2007 and 2013 FIRMs with the python function 'numpy.random. uniform.' Next, initial water levels were estimated as the elevation values from the original DEM along the inland extent of the 2007 FIRM 1% flood zone. The initial water levels were iteratively raised by 0.01-m and the linear addition by expansion method was used to flood additional areas where the cumulative water levels were greater than the DEM values. Additional areas were flooded beyond the 2007 FIRM, but within the 2013 FIRM 1% flood zone, until the randomly sampled total flood area criteria was met. Generating intermediate present-day 1% flood zone realizations based on the DEM was more realistic than realizations based on arbitrary distances between the 2007 and 2013 FIRM 1% flood zone extents. Figure 4 illustrates the intermediate presentday 1% flood zone realizations between the 2007 and 2013 FIRMs.
Storm surge is typically a larger factor in coastal flooding than wave action for locations with wide and gently sloping continental shelves (Walsh et al. 2012), such as off the U.S. East Coast. Accordingly, this study refers to the initial water level heights in the future flood model derived from the 1% flood zone realizations simply as storm surge heights. The 500 generated present-day 1% flood zone realizations, with a uniform distribution, represented the present-day storm surge uncertainty incorporated in the future flood model. Intermediate, present-day 1% flood zone realizations were randomly generated between these two FIRMs based on the DEM; darker colors between the two FIRMs represent areas more likely to be in the present-day 1% flood zone realizations.
This study incorporated local sea-level rise projections from the NYC Panel on Climate Change (NPCC) 2015 report in the future flood model (Horton et al. 2015). The NPCC 2015 report aggregated individual components of sea-level rise to generate low (i.e., 10 th percentile) and high (i.e., 90 th percentile) sea-level rise estimates (Horton et al. 2015). The individual components included global thermal expansion, local changes in ocean height, loss of ice from Greenland and Antarctic ice sheets, loss of ice from glaciers and ice caps, gravitation, rotational, and elastic 'fingerprints' of ice loss, vertical land movements/glacial isostatic adjustments, and land-water storage (Horton et al. 2015). Refer to Table 2.2 in Horton et al. (2015) for the NPCC specific estimates of sea-level rise for the 2020s, 2050s, 2080s, and 2100.
The NPCC sea-level rise projections were referenced to a 2000-2004 baseline. This study modified the NPCC sea-level rise projections to be relative to 2014, which is the year of the LIDAR collection that generated the DEM. The relative sea-level trend from a nearby NOAA Tides and Current station (Sandy Hook, New Jersey, NOAA 2018) indicated a linear rate of 0.406-cm sea-level rise per year. Accordingly, the year 2002 was considered the middle of the baseline, and the NPCC low and high estimates were reduced by 4.872-cm to approximate the future sea-level rise projections relative to 2014. Two future sea-level rise projections were derived for each decade through 2100 from the NPCC low and high estimates with 2 nd degree polynomials using the python function 'numpy.polyfit', to represent the 'best-case' and 'worst-case' scenarios, respectively. 498 additional, intermediate sea-level rise projections were then derived randomly between the low and high estimates with the python functions 'numpy. random.uniform' and 'numpy.polyfit' (Figure 5). The 500 generated sea-level rise projections, with a uniform distribution, represented the sea-level rise uncertainty incorporated in the future flood model.

Future flood model
The future flood model implemented the GIS-based, static linear addition by expansion method with Python Version 2.7, and, primarily, the NumPy and Scientific Python (SciPy) packages. First, initial storm surge heights were estimated as the elevation values from the original DEM along the inland extent of the randomly generated present-day 1% flood zone realization. A randomly generated sea-level rise projection was then added to these initial water heights and additional inundated areas were determined iteratively with array convolution and arithmetic functions from the SciPy and NumPy packages, respectively ( Figure 6).
Starting at the inland extent of the randomly generated 1% flood zone realization, the water heights of the dry cells in the randomly generated DEM realization were calculated as the average water height of adjacent inundated cells. The dry cells then became inundated if their elevations were less than these calculated average water heights. This process was repeated until no new dry cells were inundated. The flood model recorded the year at which the calculated average water heights from the randomly generated present-day storm surge heights and future sea-level rise projection became greater than each DEM realization cell. An 8-neighbor water connectivity algorithm defined the cell adjacency and was implemented to avoid incorrectly inundating inland local depressions not connected to the present-day or modeled future ocean. The 8-neighbor algorithm was implemented to model the 'worst-case' scenario of water connectivity, as the 4-neighbor algorithm is a conservative estimate of water connectivity and can potentially result in less inundation due to terrain barriers of water flow.

Monte Carlo simulations
This study implemented a Monte Carlo approach to incorporate random combinations of input data realizations and create a 500-member, future flood model ensemble (e.g., Heuvelink et al. 1999, Cooper and Chen 2013, Leon et al. 2014, Schmid et al. 2014, West et al. 2018. Each of the 498 intermediate simulations consisted of randomly generated input data realizations within the minimum and maximum DEM, sea-level rise, and storm surge values in Table 2. The year inundated for each DEM cell was recorded in each simulation, resulting in a distribution of years inundated for each cell based on the random combinations of the input data realizations. Spatiotemporal statistical products were then generated from the cell-level distributions of the year inundated to indicate the likelihood of the future 1% flood zone extent in Tottenville, NYC. The effect of the input data realizations on the modeled total area in the future 1% flood zone was also quantified with Spearman's rank correlation coefficients (Spearman 1904).

Interactive web map to host statistical products
An interactive web map, UncertainSeas.com, visualizes the spatiotemporal statistical products that depict the future 1% flood zone in Tottenville, NYC. One statistical product depicts the annual probability of the future 1% flood zone for each decade through 2100. The annual probability was calculated for each DEM cell as the number of times the cell was flooded by the specified year out of the 500-member ensemble, and larger probabilities indicate areas more likely to be located within the future 1% flood zone.
Another statistical product, alternatively, depicts the year in which each DEM cell exceeds a given annual probability of being in the future 1% flood zone. This probability exceedance year statistical product can be especially useful for local coastal management policies because communities can have much different risk tolerances to coastal flooding. This statistical product can indicate when a community should implement policies to reduce the vulnerability of an area to flooding, depending on their risk tolerance. For example, 'Community Z' will implement policies to reduce the vulnerability of an area to coastal flooding only when there is a 95% chance of the area being in the future 1% flood zone (i.e., they have a high-risk tolerance). Conversely, 'Community Y' is risk-averse and will implement policies to reduce the vulnerability of the same area when there is a 5% chance of the area being in the future 1% flood zone. These different risk tolerances would result in 'Community Y' implementing policies to reduce the vulnerability of an area to coastal flooding earlier than 'Community Z.' Table 2. Input data in the Monte Carlo simulations, and their 'best-case' and 'worst-case' values that delineate the future minimum and maximum flood zone extents, respectively. 498 additional, random combinations of each input data source realization within the minimum and maximum flood extent values generated the 500-member ensemble. The DEM realizations were randomly sampled from a normal distribution, whereas the sea-level rise and storm surge realizations were randomly sampled from a uniform distribution between the minimum and maximum values in this

Results
The future flood model for Tottenville, NYC incorporated multiple input data uncertainties in a GIS, probabilistic framework. Figure 7 indicates the Tottenville land area in the future 1% flood zone through the year 2100 for each of the 500 members of the Monte Carlo ensemble. The area in the future 1% flood zone varies between approximately 1.1 to 1.8 km 2 by the year 2100 and corresponds to the input data values that produce the minimum and maximum flood extents in Table 2, respectively. Figure 8 illustrates the effect of the input data realizations on the total area in the future 1% flood zone for every decade through 2100. The storm surge realizations have the largest correlation coefficients with the total area in the future flood zone between 2020-2050. After 2050, the sea-level rise realizations have larger correlation coefficients with the total area in the future 1% flood zone due to larger uncertainty in sea-level rise projections in these more distant decades. The correlation coefficients between the storm surge and sea-level rise realizations and the total area in the future 1% flood zone all have p-values less than 0.001 for each decade through 2100, except for the sealevel rise realizations for the year 2020. The DEM realizations have much smaller correlation coefficients with the total area in the future 1% flood zone and the p-values are greater than 0.37 for each decade through 2100.
Spatiotemporal statistical products derived from the Monte Carlo ensemble hosted on the interactive web map, UncertainSeas.com, address limitations of the non-spatial metric of the total land area in the future 1% flood zone. These statistical products indicate areas with high elevations in the center of the neighborhood should remain protected from coastal inundation into the distant future (i.e., 2100), even with the 'worst-case' DEM, sealevel rise, and storm surge realizations that produced the maximum flood extent in Table 2. Conversely, low-lying areas, such as along the northern and southern coasts of Tottenville, are already prone to coastal inundation, and the likelihood of adjacent areas being in the 1% flood zone increases in future decades, even with the 'best-case' realizations of the input data in Table 2.
Probability of future flood zone Figure 9 visualizes the annual probability of coastal locations in Tottenville being in the future 1% flood zone for the years 2020, 2060, and 2100, respectively. The inland extent of the future 1% flood zone remains relatively constant over time due to high elevations in the center of Tottenville that approach~30-m above NAVD 88, however, the likelihood of many coastal locations being in the future 1% flood zone increases over time due to projected sea-level rise.
The area in the box shown in Figure 9 has larger uncertainty in the future 1% flood zone extent due to low-lying elevations with smaller terrain slope. The probability of the future 1% flood zone in the year 2100 in this highlighted area is shown in the left side of Figure 10. The horizontal band of probabilities that range from 0 to 100% in the left side of Figure 10 illustrates the uncertainty of the future flood zone extent. The right side of Figure 10 indicates the annual probability of being in the future 1% flood zone for three locations (A, B, and C) through 2100. The annual probability of being in the future flood zone increases over time for these three locations due to projected increases in sea-level rise.  Probability exceedance year Figure 11 illustrates the year in which each DEM cell in this same highlighted area exceeds the 95%, 50%, and 5% annual probability of being in the future 1% flood zone. These different highlighted probabilities correspond to potential differences in a community's risk tolerance to coastal flooding. For example, the highest risk tolerance to flooding, the 95% probability, results in most areas not being in the future flood 1% zone until after 2100. Conversely, the lowest risk tolerance, the 5% probability, results in many of these same areas being in the 1% flood zone in the near future. Temporal information on the expansion of future flood zones is critical for coastal management policies, and the individual community's risk tolerance determines when to implement policies to reduce the vulnerability of an area to coastal flooding.

Uncertainseas.com
The interactive web map, UncertainSeas.com, visualizes the spatiotemporal statistical products and provides a suite of decision-making tools to inform coastal management policies in Tottenville (Figure 12). The web map hosts the statistical products that indicate the DEM cell-level annual probability of being in the future 1% flood zone for every decade from 2020 through 2100. The web map also hosts the novel statistical products that spatially depict the year at which various annual probabilities of being in the future 1% flood zone are exceeded in Tottenville (i.e., exceeds the 5%, 25%, 50%, 75%, and 95% probability).

Discussion
Studies that model future coastal inundation typically do not incorporate all major sources of input data uncertainties, i.e., the storm surge, sea-level rise, and DEM uncertainties. The probabilistic framework in this study incorporates these input data Figure 10. The annual probability of areas being in the future 1% flood zone in the year 2100 is highlighted for a portion of the study area (Left Panel). The band of probabilities that range from 0 to 100% illustrates the horizontal uncertainty of the future flood zone extent in this portion of the study area. The line graph (Right Panel) indicates the annual probability of being in the future flood zone increases over time for all three locations (A, B, and C), due to projected sea-level rise. Figure 11. Example of the annual probability exceedance year data product for the same portion of the study area shown in Figure 10. The maps depict the year at which each DEM cell exceeds a given annual probability of being in the future 1% flood zone. Three probabilities, 95%, 50%, and 5%, are shown to represent potential differences in a community's risk tolerance to future coastal flooding. uncertainties and accounts for possible non-linear interactions in modeling the future 1% flood zone extent in Tottenville, NYC.
The spatiotemporal statistical products hosted on UncertainSeas.com have many advantages over standard reports, paper maps, and the figures in this manuscript in depicting the future 1% flood zone and its uncertainty in Tottenville. The web map can facilitate interactive planning for a specific location, house, business, or infrastructure over time as coastal flood zones in Tottenville expand in the future due to projected sealevel rise. However, web map-users must remain aware of the limitations of the generated spatiotemporal statistical products due to inherent uncertainty in geospatial data, including the estimated storm surge, sea-level rise, and DEM uncertainties, and their incorporation in the future flood model.
The probability exceedance year statistical product hosted on UncertainSeas.com can be especially useful for individualized coastal management policies. This statistical product enables a uniform probability threshold to be established based on the individual community's risk tolerance, and the year at which areas are prone to flooding based on that probability threshold are depicted. The temporal information is important for coastal management policies, and the high temporal resolution of the statistical products, i.e., every decade, also advances previous studies, which often model only one year in the distant future, such as the year 2100.
Future flood zone in Tottenville, NYC Spatiotemporal statistical products generated from the Monte Carlo ensemble indicate the future 1% flood zone extent and its uncertainty in Tottenville, NYC. One statistical product illustrates the annual probability of being in the future 1% flood zone for each decade through 2100. This annual probability is calculated as the proportion of times DEM cells are flooded by the given year from the 500-member ensemble and should not be confused with the traditional, cumulative probability of flooding. For example, the traditional, cumulative probability of flooding indicates a building located in the current 1% flood zone has a greater than 26% chance of being flooded by at least one 1% magnitude flood over the course of the standard 30-year mortgage (Pielke, Jr. 1999, Burby 2001. This larger, cumulative probability of flooding further emphasizes the need to protect areas of Tottenville immediately, such as along the northern and southern coasts of the neighborhood. Figure 8 indicates that the storm surge realizations have the largest correlation coefficients with the total area in the future flood zone in the coming decades (i.e., 2020-2050). This is not surprising, given the accepted legal challenge to the 2013 Preliminary FIRM, and the large, present-day storm surge uncertainty illustrated in Figure 4. Furthermore, there is little consensus on how storms in NYC will change in future climates (Orton et al. 2015, Sweet et al. 2017. There is medium confidence in a projected increase in the intensity of hurricanes in the North Atlantic (Sweet et al. 2017). However, there is low confidence in the projected increase in frequency of intense Atlantic hurricanes, and the associated amplification of flooding could be offset or amplified by changes in overall storm frequency or tracks (Sweet et al. 2017). Future sea-level rise can also reduce the rate at which low-lying areas drain, increasing the likelihood of flooding from rainfall (Titus et al. 1987), which was not incorporated in this study. Potential changes in future storm climatology and additional flooding from rainfall provide justification for incorporating the 'overly pessimistic' 2013 Preliminary FIRM in the future flood model.

Effect of input data uncertainties
The 2007 and 2013 FIRMs have differences in the flooded elevations of greater than 1-m in the same, general vicinity along the southern coastline of Tottenville. The differences between sea-level rise projections in the near-future and the DEM realizations are much smaller than this 1-m difference in the storm surge height representations. For example, in the year 2025, the difference between the low and the high sea-level rise estimates from the NPCC 2015 report is~0.2-m (Horton et al. 2015). Differences in DEM realizations in Tottenville are also relatively small, as the average cell-level uncertainty at one standard deviation is only~0.05-m due to modern, LIDAR technology (Amante 2018a). Incorporating DEM uncertainty in flood models is more important for coastal communities with topographic elevations mapped with older, less accurate technologies, such as with Shuttle Radar Topography Mission (SRTM) interferometric radar data (Rabus et al. 2003). For example, SRTM elevation products have absolute height errors of 5.6-m at the 90% confidence level for the continent of Africa (Rodriguez et al. 2005). Differences in DEM realizations would likely have a larger effect on the modeled future flood zone for coastal communities mapped with these older, less accurate technologies.
The effect of the input data realizations on the modeled total area in the future 1% flood zone illustrated in Figure 8 also changes over time. The sea-level rise uncertainty increases over time due to increasing uncertainty in the magnitude of ice sheet melt (Horton et al. 2015), and, consequently, the sea-level rise realizations have the largest correlation coefficients with the total area in the future 1% flood zone after 2050. It should be noted sea-level rise uncertainty and DEM uncertainty are intertwined, as the sea-level rise estimates in the NPCC 2015 report consider future, local land subsidence and its uncertainty (Horton et al. 2015).
The generated spatiotemporal statistical products illustrate the uncertainty in the future flood zone extent that results from the uncertainties in the input data sources and from the terrain itself. Figure 9 indicates there are portions of Tottenville with larger horizontal bands of uncertainty in the future flood zone extent. This occurs in low-lying areas with smaller terrain slope, such as the area highlighted along the southern coast of Tottenville in Figures 3, 9, 10, and 11. In such areas, small changes in the modeled land and water heights within the estimated uncertainties of the data sources results in larger uncertainty in the future flood zone extent. Conversely, areas with larger terrain slope, such as along the western coast of Tottenville highlighted in Figure 3, have a much narrower band of horizontal uncertainty in the future flood zone extent. These results support previous findings on the impact of terrain slope on the horizontal uncertainty of the future flood zone extent (Gesch 2013, West et al. 2018. Incorporating input data uncertainties in coastal flood models is especially important for communities with heterogeneous terrain, such as Tottenville, which can result in low-lying areas with smaller terrain slope having greater uncertainty in the future flood zone extent.

Limitations and future research
The estimation and incorporation of the input data uncertainties in the coastal flood model can be improved in future research. The DEM cell-level uncertainty estimation is limited by the incorporated measurement uncertainty from the LIDAR dataset's metadata. This uniform measurement uncertainty of 5.3-cm RMSE z should be refined in future research using accurate ground control measurements because LIDAR measurement uncertainty is typically correlated with land cover and terrain (Su and Bork 2006, Bater and Coops 2009, Coveney and Fotheringham 2011, Spaete et al. 2011, Leon et al. 2014, Goulden et al. 2016, West et al. 2018). This study did, however, implement methods from Amante (2018b) that partially incorporate land cover and terrain effects on the magnitude of DEM cell-level uncertainty. These methods consider the number of measurements per DEM cell and the subcell measurement variance, resulting in a spatially-varying estimate of DEM cell-level uncertainty that is correlated with land cover and terrain slope, respectively (Table 1). Accurate ground control elevation measurements are also needed to identify any systematic, vertical errors in the DEM, which were not rigorously quantified in this study, but could cause differences in the modeled future coastal flood zone. Future morphologic changes, such as coastal erosion or accretion, were also not incorporated in the estimated DEM uncertainty. Additional components of DEM uncertainty are difficult, if not impossible to incorporate, including future drainage modifications (e.g., canals, ditches, culverts), and engineered barriers (e.g., levees, seawalls, flood gates, Gesch 2013).
The effect of the DEM spatial resolution on the modeled future flood zone was also not evaluated. The DEM spatial resolution must resolve important terrain features, such as channels or levees, which can enhance or impede the flow of water, respectively, by considering the Nyquist-Shannon sampling theorem. The Nyquist-Shannon sampling theorem states that terrain features must have dimensions at least twice the spatial resolution to be resolved by the DEM (Nyquist 1928, Shannon 1949, McBratney et al. 2003, Hengl 2006. Previous research indicates the modeled flood area generally increases with coarser-resolution DEMs (e.g., Saksena andMerwade 2015, Hsu et al. 2016), which is a manifestation of the scale effect of the modifiable areal unit problem (Openshaw 1977, 1983, Fotheringham and Wong 1991. DEMs with finer spatial resolutions than the 1/9 th arc-second (~3-m) DEM in this study that are supported by the LIDAR dataset's post spacing of 0.7-m could be incorporated in future research to investigate the sensitivity of the modeled future flood zone to the DEM spatial resolution.
Another avenue of future research is to compare DEM realizations created using the different methods described in Wechsler (2007), including various spatial moving averages, pixel swapping, spatial autoregressive models, and sequential Gaussian simulation, and determine the sensitivity of the future flood model to these different DEM realizations. These future research endeavors to improve the estimation and incorporation of DEM uncertainty in coastal flood models should provide greater benefits to communities with topographic elevations mapped with older, less accurate technologies than topographic LIDAR.
GIS-based, static methods of modeling future storm surge inundation enhanced by sea-level rise typically use only a topographic DEM. Dynamic methods directly use a bathymetric DEM in the hydrodynamic storm surge models, and Amante (2018b) indicates the vertical uncertainty of offshore bathymetry can be much larger than the uncertainty of nearshore topography. The depth and shape of the seafloor can influence storm surge heights, and relatively wide and shallow continental shelves can amplify storm surge (McInnes et al. 2003, Walsh et al. 2012. The bathymetric uncertainty is often not incorporated in hydrodynamic storm surge models due to computational expense, and, therefore, this unincorporated bathymetric uncertainty provides further justification for using the 2013 Preliminary FIRM to represent additional storm surge uncertainty. Future research should model future flood zones with dynamic methods that use bathymetric DEMs directly to validate the GIS-based, static methodology implemented in this study. Future research should also investigate optimal mesh node locations for storm surge models based on the vertical uncertainty and horizontal precision of the underlying DEM (Amante 2018a), and the implications of storm surge mesh node locations in comparisons between dynamic and static methods of modeling future coastal flood zones.
This study focused on modeling the future 1% flood zone in Tottenville. Future research should also identify areas most vulnerable to both the physical hazard of future coastal flooding, as well as the social vulnerability to coastal flooding (e.g., the SOcial Vulnerability Index (SOVI), Cutter et al. 2003). Important social factors to consider include the community's experience with coastal flooding, and the community's ability to respond to, cope with, recover from, and adapt to coastal flooding, which is influenced by the community's economic, demographic, and housing characteristics (Cutter et al. 2003). The social data uncertainties should also be incorporated in a similar probabilistic framework to model the areas and people of Tottenville that are most vulnerable to future coastal flooding and to quantify its uncertainty.
One large source of uncertainty that is difficult to quantify is a community's potential response to expanding flood zones, such as by engineering flood mitigation structures or with planned retreats from coastlines (Fleming 2018, Lempert 2018). An important factor in a community's potential response is the perceived loss of culture and identity associated with unique cultural heritage sites and a sense of place (Lempert 2018). Collaborations between physical and social scientists are needed to incorporate all these sources of uncertainty in both physical and social models and to determine the overall vulnerability of communities to future coastal flooding.

Conclusions
This study implemented a Monte Carlo approach to model the future 1% flood zone extent in Tottenville, NYC in a probabilistic, GIS framework. The future flood model incorporated the uncertainties in a topographic DEM, present-day 1% flood zone from storm surge, and future sea-level rise projections. Generated spatiotemporal statistical products from the 500-member Monte Carlo ensemble indicate the future 1% flood zone extent and its uncertainty in Tottenville.
The spatiotemporal statistical products indicate areas with high elevations in the center of the neighborhood should remain protected from coastal inundation into the distant future (i.e., 2100). Conversely, low-lying areas, such as along the northern and southern coasts of Tottenville, are already prone to coastal inundation, and the likelihood of being in the 1% flood zone increases in adjacent areas in future decades. The spatiotemporal statistical products also indicate the uncertainty in the future flood zone extent that results from the input data uncertainties and from the terrain itself. Small changes in the modeled land and water heights within the estimated uncertainties of the data sources result in larger uncertainty in the future flood zone extent in low-lying areas with smaller terrain slope.
An interactive web map, UncertainSeas.com, visualizes the likelihood of the future expansion of the 1% flood zone in Tottenville. The spatiotemporal statistical products hosted on the web map, combined with additional information on social vulnerability, can inform coastal management policies to reduce the overall vulnerability of the people, property and economy of Tottenville, NYC to future coastal inundation.