Role of wetlands in reducing structural loss is highly dependent on characteristics of storms and local wetland and structure conditions

Coastal communities in New Jersey (NJ), New York (NY), and Connecticut (CT) sustained huge structural loss during Sandy in 2012. We present a comprehensive science-based study to assess the role of coastal wetlands in buffering surge and wave in the tri-state by considering Sandy, a hypothetical Black Swan (BS) storm, and the 1% annual chance flood and wave event. Model simulations were conducted with and without existing coastal wetlands, using a dynamically coupled surge-wave model with two types of coastal wetlands. Simulated surge and wave for Sandy were verified with data at numerous stations. Structural loss estimated using real property data and latest damage functions agreed well with loss payout data. Results show that, on zip-code scale, the relative structural loss varies significantly with the percent wetland cover, the at-risk structural value, and the average wave crest height. Reduction in structural loss by coastal wetlands was low in Sandy, modest in the BS storm, and significant in the 1% annual chance flood and wave event. NJ wetlands helped to avoid 8%, 26%, 52% loss during Sandy, BS storm, and 1% event, respectively. This regression model can be used for wetland restoration planning to further reduce structural loss in coastal communities.

Wave. During Sandy, four wave-buoys within the domain recorded the wave data: two at the apex of NY Bight, and another two in Long Island Sound (LIS). The significant wave height H sig and peak wave period T p simulated by the Simulating Waves Nearshore (SWAN) model in the CH3D-SSMS were compared with measured data and Wave Watch III (WW3) operational run results (SI Fig. S5) 25,26 . SWAN more accurately captured the evolution of H sig in NY Bight and LIS. Maximum significant wave height over land was as high as 2.17 m at NY Bight and 2.60 m along NJ coast.
Impact of wetland on surge and wave during super storm sandy. To estimate the value of the coastal wetlands in reducing flood, we calculated the following four metrics for inundation: the Average Inundation Height ( AIH ), the Maximum Inundation Height ( MIH ), the Total Inundation Area ( TIA ), and the Total Inundation Volume ( TIV ) with and without wetlands. The TIA and TIV are defined as where H max x, y and H 0 x, y are the maximum water level and the land elevation at land cells x, y , respectively 11,12 . The wave analysis was carried out by calculating the average wave height ( AWH ), the maximum wave height ( MWH ), and the total wave energy ( TWE ) which is defined as where H rmsmax is the maximum root-mean-square wave height over the flooded land. The wave energy, instead of wave height, is more directly related to the wave-induced structure loss. We define the relative inundation reduction ( RIR ) as the difference in TIV (value without wetland minus value with wetland), divided by the TIV with wetland. The relative wave reduction ( RWR ) is defined accordingly using the TWE . The inundation and wave analysis were carried out at the regional level (SI Table S2) and zip-code resolution (Figs. 1, S6, S7).
If the wetlands were absent during Sandy, the RIR of the entire model domain would have been 4% and the RWR would have been 19%. Breaking the entire model domain into 6 regions (Table S2): South NJ (SNJ), Central NJ (CNJ), North NJ (NNJ), CT, Long Island (LI), and Mainland NY (MNY), it was found that the RIR during sandy was low (< 10%) for all regions, while the RWR was moderate for CNJ and SNJ (10-25%) where the wetland cover was high, and low for the other regions. Although CT ranked 3rd in wetland cover, its RWR ranked 4th during Sandy due to Sandy's landfall location and the sheltering of waves provided by LI which ranked 3rd.
Model simulations show that, during Sandy, if the wetland in the entire region were completely made up of woody wetland, TIV and TWE would decrease by 9.1% and 12.4%, respectively, while increase by 1.2% and 2.2%, respectively, if the wetland were composed of marsh only. RIR and RWR for the all-marsh wetland are 6.5% and 19.1%, respectively, while those for the all-woody-wetland are 17.6% and 36.8%, respectively. These results demonstrate that the woody wetland is much more effective than the marsh in buffering flood and wave.
Benefit of wetland on surge and wave during a Black Swan (BS) storm. We consider a rare "Black Swan" (BS) storm (with a 0.0034 annual frequency vs. 0.0014 for Sandy 27 ) which made landfall in Staten Island with a storm surge that surpassed the maximum surge during Sandy with maximum coastal values of 7.17 m at NY Bight and 4.51 m at NJ coast of (Fig. S8). On the other hand, maximum significant wave heights over land were not greater than those of Sandy: 1.99 m at the NY Bight and 2.12 m on the NJ coast, but average wave height was higher than during Sandy (Fig. S8). During the Black Swan storm, wetlands created the biggest RIR and RWR . MNY and NNJ experienced significant RIR , while the remaining regions experienced low RIR . In contrast to Sandy, in which CT, MNY, and NNJ had less RWR , for the Black Swan storm these regions experienced high (> 50%) RWR , while LI. CNJ and SNJ had significant (~ 25-50%) RWR . The regions with the highest RIR and RWR were closest to the storm landfall location.
Benefit of wetland on surge and wave for the 1% annual chance flood. The maximum 1% annual chance maximum flood elevation (Fig. S9) was 5.15 m for NY Bight while 5.34 m for the NJ coast. The 1% annual chance maximum wave height (Fig. S9) overland is 2.24 m at NY Bight and 1.90 m at NJ coast. All the regions had moderate RIR except for SNJ which had a significant RIR due to highest (25.4%) wetland cover (see Table S2) and NNJ experienced significant RWR while the remaining regions had moderate values. CT ranked 2nd in RIR although it has less wetland cover than CNJ because the mostly woody wetlands in CT are more effective in buffering storm surge than the marshes in NJ and NY. On the other hand, CT ranked 3rd in RWR due to the blocking of offshore wave energy by LI. RIR and RWR are found to be functions of storm characteristics, wetland type, and cover, and local conditions. Fig. S10 shows the percent wetland cover, RIR (relative TIV reduction), and RWR (relative wave energy reduction) in six regions (New York, North New Jersey, Long Island, Connecticut, Central www.nature.com/scientificreports/ New Jersey, and South New Jersey) during 1% annual chance events. As the wetland cover increases from less than 5% (NY, NNJ, and LI) to more than 10% in CNJ and SNJ, RIR and RWR generally increase, showing the increasing role of wetland in reducing inundation and wave. Relative reduction in inundation and wave energy are modest, between 10 and 30%. NNJ has properties behind the relatively sparse marsh, followed by woody wetland which protects properties behind them. Connecticut has less wetland than Central Jersey, but the mostly woody wetland is more effective in reducing flood and wave.
Benefit of wetlands on reducing residential structure loss. The monetary loss of residential structures was estimated using the simulated inundation and wave results while employing damage functions from the United States Army Corps of Engineer (USACE) North Atlantic Comprehensive Coastal Study (NACCS) and was validated using the NFIP building loss payouts aggregated by zip-code 21,24 . Direct simulation of waveinduced damage requires understanding and calculation of wave loads on structures using a depth-and phaseresolving model, which is beyond the capability of the models used in this study. Therefore, in this paper, we did not directly simulate wave-induced damage, but are accounting for wave-induced change in total water level which results in increase in estimated damage based on depth-damage functions. Overall, 96 coastal zip-codes in the state of NJ were used to validate the estimated loss. The model showed a correlation coefficient (CORR = 0.69) between simulated structure losses and NFIP payouts (Fig. 2). In NJ, as of 2019, the total NFIP payout was $3.9 billion USD, in comparison to the estimated total structure loss of $3.6 billion USD (SI Table S3), with an absolute error of 7.7%. This good agreement, plus the good agreement between the simulated and observed surge and wave reported earlier, confirms the validity of our "dynamics-based" loss assessment. We define the structural loss reduction ( SLR ) as the structural loss without wetland minus the structural loss with wetland, and the relative structural loss reduction ( RSLR ) as the ratio between SLR and the loss with wetland. A state-level analysis of structure loss in NJ showed a RSLR of 8.5%, 26.0%, and 52.3% for Sandy, BS storm, and 1% chance flood/wave, respectively (Table S3). Analysis of losses due to flood and wave indicated that for Sandy and the 1% event, most of the loss came from flood, while most of the loss in the BS storm came from waves. Avoided wave-induced loss was comparable to the avoided flood-induced loss during the BS storm and the 1% event, but much higher during Sandy, suggesting that NJ wetlands are more effective in reducing wave-induced loss vs. flood-induced loss. Results from the zip-code scale analysis during Sandy (Fig. 3) showed www.nature.com/scientificreports/ less variability: for most north NJ, RSLR ranged from 10 to 100% except for those along the Hudson River that had small increased loss (negative avoided loss). Most zip-codes in SNJ had RSLR between 5% and more than 100%. RSLR during the BS storm ( Fig. S11) was more notable: ~ 50-100% for north NJ, 0-100% for central NJ.
The above results showed that the value of coastal wetlands for flood/wave protection varies significantly with the storm. While wetlands may be more effective in reducing wave loss in some storms but flood loss in other storms, they may be ineffective in extreme storms. The 1% annual chance flood and wave event, which resulted from an ensemble of many less extreme but more frequent storms 28 , provides a more reasonable integrated scenario for the loss analysis. This is similar to the preferred use of the 1% flood map, instead of the flood map associated with a single design storm, for assessing the flood risk in any coastal region.
As shown in Fig. 5, in zip-codes with larger wetland coverage area in SNJ, wetlands could only prevent < $25,000 annual loss (loss in the 1% event divided by 100), while wetlands in zip-codes with medium wetland coverage area in central NJ could prevent > $60,000 annual loss. On the other hand, wetlands in north NJ zip-codes with smaller wetland coverage area would increase the loss. The annual loss in most of the zip-codes is < $10,000 except in one zip-code where the annual loss is ~ $50,000.
To explain the significant spatial variation in SLR by wetlands, we developed a generalized linear regression model using the normal distribution with the logit link function, which was fitted with data for all coastal zipcodes during three flood events and with and without wetlands. Percent wetland cover, percent at-risk structural value, and average wave crest height in the zip-code were used as predictors to estimate the percent structural loss (PSL = structural loss value divided by at-risk structural value) in the zip-code (Table 1). For Sandy, the regression model had a significant correlation (R 2 = 0.75, p < 0.001). For the BS storm and the 1% annual chance flood, there was a strong correlation with R 2 = 0.87 and R 2 = 0.84, respectively, both with p < 0.001. Using the three storm scenarios, each with and without wetland, the regression fit (Fig. 6) had a significant correlation with R 2 = 0.75 and p < 0.001.

Discussion
In this paper, we presented the benefits of wetlands for reducing storm surge, waves, and structural losses during Hurricane Sandy. Recent papers on this topic, such as NAR17 and LAT19, focused on the effects of wetlands on reducing structure loss due to storm surge but did not include the effect of wetlands on reducing wave-induced structural loss 13,15 . Estimated loss in NAR17 was not compared to NFIP payouts. LAT19 pointed out that vegetation has a dampening effect on wave energy, but waves were not included in their empirical NFIP pay-out model due to lack of data 15 . NAR17 did not include the wave in their hydrodynamic model explicitly 13 , although it is known that wave setup could contribute up to 20% of the peak surge 29 and wave can introduce additional loss. Here, we showed how valuable these wetlands are in terms of reducing wave energy which could damage structures and even cause deaths. Wetland effect on attenuating storm surge has been widely studied and the attenuation effect increases with the wetland thickness, density, and height 11 . Our results showed that regions with thicker, denser, and taller wetland cover reduced more TIV and TWE during Superstorm Sandy. Regions www.nature.com/scientificreports/ with highest wetland cover, hence more TIV and TWE reduction, were found in SNJ, CT, LI, CNJ, NNJ and MNY, in descending order. LAT19 observed that the extent of wetland on Barnegat-Great Bay did not have enough spatial extent and height to dissipate the Superstorm Sandy surge 15 . This is consistent with our finding that wetland dissipated TIV by − 0.5% to 0.5%. Like NAR17, we found that landward zip-codes benefited from the surge attenuation at coastal zip-codes 13 . Moreover, we found that wave attenuation at coastal zip-codes helped to reduce wave energy at landward zip-codes. Our study showed notable TWE dissipation in SNJ where the landward zip-codes experienced a cumulative reduction in surge and wave which translated into less structural loss (Fig. 3). For example, in zip-code 08330 (Hamilton Township), wetland helped to avoid 41% of structural loss, far less than the 139% reported in NAR17 13 . In contrast, we found that zip-code 08226 (Ocean City) had only a ~ 3% loss avoidance. The annual avoided losses for zip-code 08330 and 08226 were ~ $11,000 and ~ $3.6 million, respectively. This showed that percent avoided loss is not a good index for characterizing the wetland value for protection against storm surge and waves. Zip-code with low total structural value and low avoided loss, such as 08330, was estimated by NAR17 to have a misleadingly high value of percent avoided losses.
Our estimated loss compared reasonably well with the NFIP loss payouts. The squared correlation coefficient of 0.47 is much higher than the value of 0.184 of LAT19 who found no evidence that increasing marsh buffer distance would reduce the amount of NFIP payouts for the nine south NJ townships 15 . LAT19's suggestion that NFIP payouts will increase by $4.05 per foot of marsh width was not supported by our dynamics-based economic analysis 15 .
Our economic analysis slightly under-estimated the total NFIP payouts, probably because of some conservative assumptions. Our hydrodynamic model did not consider major rivers in NJ and did not include precipitation. For simplicity, we used the 2D version of CH3D-SSMS with the Manning's n coefficient representing the friction by vegetation, land, and structures, which likely overestimated the friction effect and reduced the loss. During Sandy, marsh stem height-usually lower than 30 cm in south NJ according to the 2012 land cover database of www.nature.com/scientificreports/ the NJ Department of Environmental Protection (NJDEP)-was much lower than the storm surge, therefore the actual marsh benefits were probably somewhat lower than that represented by the 2D model 30 . The effect of wetland structure, as well as marsh distribution, on flood and wave dissipation, could be better represented by using a 3D vegetation resolving model, e.g., Sheng et al. 11 , if detailed wetland data were available. Our structural loss model did not include unavailable basement data which could have increased the losses. The Manning's n used in our simulation was 0.045 which is smaller than the value of 0.068 for rigid vegetation 31  During the BS storm, NJ coast experienced a slightly higher surge but much higher wave, both significantly buffered by the wetlands. Hence, the removal of wetlands would significantly increase the flood, wave, and loss along the NJ coast than during Sandy. This showed that wetlands' buffering capacity depends significantly on the wetland structure/distribution, the characteristics of the hurricane, and the specific location of interest. To assess the overall buffering capacity of tidal wetlands over a large region, an ensemble of storms should be considered instead of 1-2 specific storms. To generate the 1% flood and wave event, 300+ "optimal storms" were generated with the JPM-OS method and the CH3D-SWAN model based on the ensemble of storms generated by a statistical hurricane model 26 . Flood and wave associated with these optimal storms are being used to develop a Rapid Forecast and Mapping System (RFMS) 33 for NJ/NY/CT region, which can be used to predict real-time surge and wave including wetland effects during any hypothetical storm.
The regression analysis shows that wetlands are more effective in reducing losses when the average wave crest height is lower, which confirms that wetlands' buffering capacity is more effective for low-intensity and more frequent storms 16 . Loss in a zip-code increased with the at-risk structural value and the average wave crest height but decreased with percent wetland cover. Zip-codes with higher structural value benefited more from wetland protection. Without any structure behind a wetland, the ecosystem service value of the wetland for flood protection is zero. The regression model developed in this study can be used to estimate the impact of changing wetland coverage on loss, e.g., doubling the wetland coverage in every NJ zip-code would decrease the average loss in Sandy from 4 to 3%.
Although coastal wetlands were found to reduce structural loss during Sandy, their percent loss reduction was found to be very modest (− 10% to + 50% for most except one zip-code) compared to 22-139%, reported      This study showed that the percent structural loss (PSL) varies significantly across a coastal region depending on percent wetland cover, at-risk structural value, and storm characteristics 11 (represented by average wave crest height). As shown in Fig. 7, the PSL for NJ is between low (< 10%) and moderate (< 25%). These maps, together with the regression model can be used to develop wetland restoration plans (e.g., increasing wetland coverage area) to reduce structural loss in the future. Adding the regression model to the RFMS will enable forecasting of structural loss during a hurricane.
Strong winds during storms can directly damage coastal wetlands, resulting in reduced buffering capacity of wetlands for coastal protection. Further studies are needed to fully understand the complex interactions among wetlands (with different type, structure, density, and coverage area, etc.), storm surge and wave, erosion, wind, and sea-level rise (SLR). With 1-2 m of SLR projected by 2100, tristate tidal marshes (with a height of 0.3-1 m) can be easily overwhelmed and lose their buffering capacity 34 . The approach presented in this study can be expanded to develop and assess wetland restoration plans to mitigate future flood risk in the twenty-first century.

Methods
The role of coastal wetland in reducing coastal flood is analyzed by determining (a) the surge peak and extent, (b) the significant wave height and peak wave energy, and (c) the resulting structural losses at all flooded locations during three flooding events (Sandy, "Black Swan" storm, and the 1% annual chance flood and wave event). The surge and significant wave height were estimated using a coupled surge-wave modeling system CH3D-SSMS and the structural loss was estimated using the USACE depth damage functions 19 .
Hydrodynamic modeling. The coupled surge-wave model is based on CH3D-SSMS 29,35,36 . CH3D-SSMS couples 2D/3D storm surge and wave models of coastal and basin scales using a non-orthogonal curvilinear grid for the horizontal directions and a terrain-following sigma coordinate in the vertical direction for an accurate representation of estuaries and coastal waters. In the coastal region, the surge model CH3D is dynamically coupled to the wave model SWAN 25,37,38 . The governing modeling equations and surge-wave coupling mechanisms of CH3D-SSMS were well documented in Sheng et al. 29 .
The study domain is represented by a high-resolution curvilinear grid which consist of 1461 by 2182 cells with 40-50 m resolution in NYC and about 20 m in the low-lying land area, such as lower-Manhattan, to resolve the local complex geographic features. The grid domain covers the coasts of New Jersey, New York, and Connecticut. It also includes the Hudson River from the New York Harbor to the Federal Dam at Troy. It extends from the continental shelf to elevations that are higher than the possible extent of flooding by storm surge. The grid in CT and NY reach elevations that are greater than 200 m and in NJ it is expanded 40+ km inland or reaching elevations that are greater than 10 m. The model's bathymetry and topography were generated using NOAA's National Centers for Environmental Information (NCEI) bathymetric and topographic telescoping digital elevation models (DEMs) of the U.S. Atlantic Coast impacted by Hurricane Sandy. This DEM resolution varies from 1/9 arc-sec (3.4 m) in the coastal zone to 3 arc-sec (90 m) in the offshore. For the locations where NCEI's DEM was not available, the 1/9 arc-sec DEMs from the U.S. Geological Survey (USGS) National Elevation Dataset (NED) was used for land elevations, while NCEI's Coastal Relief Model (CRM) was used for the water cells. The Estuary Bathymetry data from the New York State Department of Environmental Conservation (NYSDEC) was used in the Hudson River between Newburgh and Troy. These elevations were interpolated to the grid cells.
For this study, a 2D version of CH3D-SSMS was implemented with the use of Manning's friction coefficient (n) to account for the friction produced by vegetation and other land cover types. Higher values of Manning's n imply greater resistance to flow due to higher horizontal shear stress throughout the water column. Land Cover types were obtained from the USGS 2011 National Land Cover Database (NLCD) which have a resolution of 30 m 39 41 . The surge boundary conditions for CH3D-SSMS were provided by an ADCIRC large-scale model to incorporate the remote effects of storm surges 9,10 . Offshore wave conditions came from NOAA's WW3 and river discharge data from the USGS. Two major rivers were included: the Hudson River at Green Island, NY (USGS 01358000), and Connecticut River at Thompsonville, CT (USGS 01184000). Two high-resolution wind fields, computed by the Regional Atmospheric Modeling System (RAMS), were applied as the meteorological forcing. A wind field that covers a smaller spatial area with a spatial resolution of 6 km and a temporal resolution of 15 min from Weatherflow, Inc. was merged with a larger spatial area wind field with a spatial resolution of 4 km and a temporal resolution of 1 h. These wind fields were also used by Wang et al. 42 successfully to simulate storm surge during Sandy.
A "Black Swan" storm. The Black Swan (BS) storm was designed as an extreme scenario of a Cat-3+ hurricane such that winds of Category 1 hurricane intensity would be felt at the Village of Piermont, approximately 40 km upstream of NYC along the Hudson River, after the landfall at Staten Island (Fig. S13). This storm would create an extreme event in NJ because it moved northward very near the NJ shoreline. The wind field was generated using the Holland wind model 43  www.nature.com/scientificreports/ 0.0034 with a return period of 294 years. Along the NJ coast, extreme Black Swan winds would generate much stronger surge and wave, and result in much higher structural loss and different wetland effects. At the Village of Piermont, the peak wind speed during the Black Swan storm would be stronger than the peak wind speed during Sandy, and the wind direction would be from the southeast instead of the east. This Black Swan wind in the Piermont region was used in a separate study to investigate the role of the Piermont marsh in protecting residential structures during storms.
Storm ensemble. Here, we consider a storm ensemble generated by Hall using the North Atlantic Stochastic Hurricane Model (NASHM) 44,45 . Then, we used JPM-OS to generate a set of optimal storms to represent all the possible storms described by Hall's original storm ensemble 28,46 . Finally, we simulated the 1% annual chance flood and significant wave height with and without the wetlands and produced the corresponding maps according to the next equations: Here, is the mean annual rate for all storms on the site, f x (x) is the joint probability density function of the storm characteristics, and the conditional probability that a storm with certain characteristics x will cause a water level height or significant wave height above η . This integral is evaluated for all possible combinations of storm characteristics. The integral in Eq. (4) is not easily determined and is usually approximated as a weighted summation of discrete storm parameters value 28,46,47 . Economic analysis. Following the procedure described in Loerzel et al. 17 , the economic analysis focused on the state of NJ for simplicity. To quantify the value of the coastal wetland for flood and wave protection, two simulations were performed, one with wetland and another without wetland. The without-wetland scenario was done by replacing Manning's n values of the wetlands with 0.02, which represents the open water 13 . For the economic analysis, we identified the residential parcels that were flooded, as well as the inundation levels and significant wave heights of each parcel. The structure value and the number of stories at every parcel are needed to use the depth damage functions from USACE 21 .
The structural data used in this study included the residential parcels in the 'Parcels and MOD-IV Composite of New Jersey' shapefile from the New Jersey Office of Information Technology, Office of GIS (NJOGIS) 48 . Basements were unaccounted for due to lack of information. The value of the structure was assumed to be uniform on any single parcel. Following Kousky and Walls 49 , a partially inundated parcel was included in the economic analysis if and only if the parcel's centroid is inside the inundation layer. The monetary value used in the analysis for each parcel was the "improvement value" reported in 2018 USD. The loss percent was calculated using the interpolated flood depth at each parcel centroid. Residential parcels without "improvement value" or building description in terms of the number of stories were removed from the analysis. The numbers of stories were rounded to the closest integer.
We combined flood and wave for economic analysis. First, the depth-limited controlling wave height is calculated as H c = min 0.78d, 1.6H sig where d is the flood, and H sig is the significant wave height. The next step is to determine if the structure is in A zone, Coastal A zone, or V zone, according to FEMA. A zone is defined as where H c is lower than 1.5 ft, Coastal A zone is where waves are between 1.5 ft and 3.0 ft, and V zone is where H c is higher than 3.0 ft 50 . Next, the total water depth above the first-floor elevation is calculated as D = d + 0.7H c − FFE for Coastal A and V zones and D = d − FFE for A zone, where FFE is the first-floor elevation. Then we used this newly updated flood elevation which accounted for the wave crest for the USACE NACCS damage functions. The wave damage functions were applied to structures in the V zone, and the flood damage functions were applied to structures in Coastal A and A zones. The structural loss of each parcel ( PL ) was calculated by multiplying the structure damaged by the parcel "improvement value". Finally, for validations purpose only, we incorporated the maximum NFIP coverage for small structures by setting the transformed structure loss PL T = min(PL, $250, 000) . The total residential structural loss was then computed by adding together the structural losses of each parcel ( PL ) in the same zip-code. Avoided losses for every zip-code were calculated by subtracting the structural losses of the wetland-absent case from that of the wetland-present case. We normalized the NFIP payouts to 2018 USD using the county housing units [51][52][53] and data reported in Weinkle et al. 53 .
It should be noted that the surge and wave were calculated by process-based dynamically coupled surgewave models and hence the simulated surge and wave contain nonlinear surge-wave interactions. On the other hand, the damage assessment was based on empirical engineering understanding of flood and wave impact on structures and used empirical formulations including empirically determined damage functions. In this study, we followed the approach of USACE and FEMA without reinventing their empirical methods and formulas, by using the 1% flood elevation and 1% wave height simulated by the coupled CH3D-SWAN model to better represent the wave effect on total flood depth and damage. While this may not be the perfect approach for damage assessment, our analysis provided a conservative damage estimation which probably slightly over-estimated the total damage by using the sum of maximum flood and maximum wave as the total maximum flood depth during each storm. Nevertheless, the good agreement between our assessed damage during Sandy and the NFIP payouts confirmed the robustness of our method. Further improvement of the empirical damage assessment may require additional data, model simulation, and parameter tuning which could be the topic of a future study.

Data availability
The topography, bathymetry, and land-use datasets used in the regional and local studies are available via the sources described above. The data on New Jersey structure value is available from https ://njogi s-newje rsey. opend ata.arcgi s.com/datas ets/406cf 68603 90467 d9f32 8ed19 daa35 9d.The NFIP payout data is available from https ://www.fema.gov/media -libra ry/asset s/docum ents/18037 4. Damage functions are available from the U.S. Army Corps of Engineers (Physical Depth Damage Function Summary Report, North Atlantic Comprehensive Coastal Study: Resilient Adaptation to Increasing Risk). All derived data such as differences in losses and flood and wave heights between the two scenarios for the regional study, are available from the corresponding author on reasonable request and may be subject to a suitable Non-Disclosure Agreement.

Code availability
CH3D and SWAN have both been described in detail in the literature cited. Any inquiry about these models should be directed to the corresponding author and TU-Delft, respectively. Wetland coverage data are available from the corresponding author upon reasonable request.