Storm surge and ponding explain mangrove dieback in southwest Florida following Hurricane Irma
A combination of airborne and satellite remote sensing data were used to quantify changes in mangrove forest structure and function from Hurricane Irma (Supplementary Fig. 1). Findings based on multi-sensor airborne data were scaled to the entire study area using estimates of forest structure and vegetation phenology derived from satellite data.G-LiHT Airborne campaignDuring April 2017, NASA Goddard’s Lidar, Hyperspectral, and Thermal (G-LiHT) airborne imager conducted an extensive airborne campaign in South Florida covering >130,000 ha. The same flight lines were resurveyed with G-LiHT eight months later, during November and December of 2017, to quantify structural changes in coastal forests of South Florida and Everglades National Park (ENP) following Hurricane Irma (Fig. 1). Lidar data was collected with two VQ-280i (Riegl USA) and synced during flight using RiACQUIRE version 2.3.7. The plane flew at a nominal height of 335 m above ground level at a pulse repetition frequency of 300 kHz to collect ~12 laser pulses per square meter. The analysis of pre- and post-hurricane conditions used 1-m resolution lidar data products (Supplementary Fig. 2) and 3-cm resolution stereo aerial and ground photos to estimate changes in vegetation structure, fractional cover, and terrain heights across the study domain. G-LiHT lidar canopy height models, digital terrain models, and estimates of fractional vegetation canopy cover (FVC) were produced using standard processing methodology21. All Level 1 through 3 lidar data products and fine-resolution imagery are openly shared through the G-LiHT webpage (https://gliht.gsfc.nasa.gov/).High resolution stereo maps of canopy heightStereo imagery from high-resolution commercial satellites can be used to estimate canopy and terrain surfaces42,43. Here, we derived digital surface models (DSMs) from DigitalGlobe’s WorldView 2 Level 1B imagery. DigitalGlobe provides these data to U.S. Government agencies and non-profit organizations that support U.S. interests via the NextView license agreement44. The spatial resolution of these data depends on the degree of off-nadir pointing for each acquisition. In this study, image resolution ranged from 0.5 to 0.7 m. We selected along-track stereopairs within the study domain to identify stereo image strips (each ~17 km × 110 km) that were nominally cloud-free over the forested domain of interest for years 2012–2013, the most recent cloud-free stereo data available for the study region prior to Hurricane Irma. The DSMs were produced using the Ames Stereo Pipeline (ASP) v. 2.5.1 on the NASA Center for Climate Simulation’s Advanced Data Analytics Platform at Goddard Space Flight Center (ADAPT, https://www.nccs.nasa.gov/services/adapt). The Worldview DSMs have been shown to accurately estimate mangrove canopy height when compared to airborne lidar and radar interferometry42,43. The processing workflow was adapted from ref. 45, and was implemented semi-global matching algorithms with a 5 × 5 correlation kernel, and a 3 × 3 median-filter applied to the output point cloud prior to producing a 1 m DSM using a weighted average gridding rule46. The ASP processing yielded five DSMs at 1-m resolution that were used to capture pre-storm canopy surface elevations.Each of the five Worldview DSMs were individually calibrated using overlapping pre-storm G-LiHT lidar data to estimate mangrove canopy heights across the study region (Supplementary Fig. 1). We sampled 1000 points within the mangrove forest cover (see mangrove classification, below) to develop a bias-correction equation between G-LiHT lidar-derived canopy heights and stereo DSM elevations (Supplementary Fig. 6). The bias-corrected canopy height models from high-resolution stereo imagery were mosaicked together to generate a 1-m resolution CHM for the entire study region (Supplementary Fig. 7). A pre-storm canopy volume was calculated by summing the 1 m × 1 m WorldView CHM for the entire region of interest. Similarly, a post-storm canopy volume was derived using the canopy damage model (see the section below), the relationship between the pre-storm CHM and the max wind speed. This analysis was conducted in ArcMap 10.7.1.Landsat mangrove forest classificationLandsat 8 Operational Land Imager (OLI) imagery was used to map mangrove cover for the southern Florida study region. The imagery was preprocessed to surface reflectance47 and clouds were masked following methods outlined in ref. 48. The Surface Reflectance Tier 1 product in Google Earth Engine was used to create a cloud-free image mosaic for 2016 based on the median values of all cloud-free images for the year for all bands (Supplementary Fig. 1).Training points were hand-selected using contemporary Google Earth imagery, field photos, and expert knowledge of the region. Twenty-four polygons covering a mangrove area of 1243 ha and 17 polygons covering a non-mangrove area of 2759 ha were identified for training regions. Within each of the two classes (i.e., mangrove and non-mangrove), 100,00 points were sampled and used for the training data in a Random Forest Classification implemented in Google Earth Engine49. The Random Forest model used 20 trees and a bag fraction of 0.5. The Landsat-based mangrove map was validated using the Region 3 species land cover map developed by the National Park Service for Everglades National Park50. The National Park species map was reclassified into mangrove and non-mangrove land cover, and 500 randomly generated points were sampled within each of the two land cover classes. The resulting error matrix indicated an overall accuracy of 90.6%.Post-storm canopy coverTime series of Landsat data were used to estimate hurricane damages of mangrove forest cover through December 31, 2017. We combined data from Landsat 7 ETM+ and Landsat 8 OLI to create a dense time series of cloud-free observations. All images were pre-processed to surface reflectance and masked for clouds using the same methods as the mangrove classification. Landsat 7 and Landsat 8 data were then harmonized to account for differences in the sensor specifications following51. We calculated the Normalized Difference Vegetation Index (NDVI) and Normalized Difference Water Index (NDWI) for each image in the collection. We calculated two reference maps from the time series of Landsat imagery (Supplementary Fig. 1). A pre-storm reference was calculated as the median value for each reflectance and index band for all cloud-free imagery in the two years prior Hurricane Irma, August 31, 2015 through August 31, 2017. Similarly, a post-storm median mosaic image was made using Landsat data between October 1, 2017 and December 31, 2017.Pre- and post-storm wall-to-wall Fractional Vegetation Cover (FVC) maps were generated using a combination of lidar-based FVC metrics and Landsat imagery (Supplementary Fig. 1). First, lidar-based FVC was binned into five classes; 0–20%, 20–40%, 40–60%, 60–80%, and 80–100% (Supplementary Fig. 7). We then collected 1000 randomly generated points in each of the five FVC classes, a total of 5000 points, to be used as training data in the Landsat classification. Here, we implemented a Random Forest Classifier using 100 trees and a bag fraction of 0.5. These steps were applied to both the pre-storm and post-storm lidar-derived FVC and Landsat mosaic image metrics. Changes between the pre- and post-storm FVC were then calculated based on the five different FVC classes (Supplementary Fig. 7). For example, a pixel with pre-storm FVC of 80–100% and a post-storm FVC of 20–40%, a reduction of three FVC classes, was assigned a drop in FVC of 40–60% (Fig. 1).Recovery times and resilienceWe estimated the time to full recovery of pre-storm mangrove green canopy cover using the time series of Landsat NDVI during the first 15-months following Hurricane Irma. The pre-storm mean NDVI layer was used as a reference, as described in the previous section. Next we calculated the NDVI anomaly for each image during the post-storm period, September 17, 2017 through December 31, 2018 (Supplementary Fig. 1). We then summed the individual anomaly values from each Landsat image and normalized by the total number of valid pixels (i.e., pixels meeting quality control measures) to estimate the average change in NDVI within the 15 months after the storm. We used anomaly values to identify mangrove forests with large decreases in the 15 months after the storm using a threshold of 0.2 for the 15-month NDVI average anomaly19,52. These areas suffered large losses of canopy material and limited new growth during the post-storm period. We used the slope in NDVI values for each pixel during 2018 to estimate the time in years to full recovery to pre-storm NDVI values, excluding data from October to December 2017 to remove delayed browning of damaged vegetation and spurious NDVI values from surface water features following the storm. Areas with a negative NDVI slope were not assigned a recovery time.We used a combination of the NDVI slope, estimated time to full NDVI recovery, and the average change in NDVI between the pre- and post-storm periods to categorize mangrove forest resilience, the potential for mangroves to rebound to pre-disturbance conditions. The specific criteria for mangrove recovery rates and mangrove damage thresholds were adapted from field and remote sensing studies, respectively6,19,25. Regions of high resilience (a combination of high resistance and resilience) were identified based on rapid recovery and/or little to no immediate impact from the storm: (1) areas that were observed to recover to pre-disturbance conditions during 2018, (2) areas that were predicted to recover within 5 years regardless of the post-storm drop in NDVI6, and (3) regions with a post-storm change in NDVI 15 years or a negative NDVI slope that occurred in regions with the largest ( >0.2) post-storm drop in average NDVI25 (Supplementary Fig. 9). The resilience class map is available online for download53.Mangrove species and elevationWe used species level maps developed by the National Park Service for Everglades National Park50 to characterize the impact of Hurricane Irma on different mangrove species. For that study, dominant species were identified through photo-interpretation of stereoscopic, color-infrared aerial imagery. Grid cells of 50 m × 50 m covering an area (Region 3) of ~100,000 ha in southwest Florida were interpreted based on the majority cover type and validated using field observations. A total of 169 vegetation cover classes were identified in this region, however, only five mangrove cover classes were considered for these analyses: Avicennia germinans (Black Mangrove), Laguncularia racemosa (White Mangrove), Rhizophora mangle (Red Mangrove), Conocarpus erectus (Buttonwood), and a single mixed species mangrove class. Mangrove forest communities were defined as the dominant diagnostic species in the upper-most stratum50. The mangrove species data were reprojected to match the Landsat resolution and the resilience maps. We used the intersection of the resilience and species extent maps to estimate the proportion of each resilience class by dominant species.The USGS National Elevation Dataset (NED) was used to estimate the soil elevation across southwest Florida. The 1/9 arc second (~3 m × 3 m) products were acquired from NED, and reprojected to Landsat resolution to estimate the proportion of each resilience class by soil elevation.Additional data and analysisModeled maximum storm surge data for Hurricane Irma were acquired from Coastal Emergency Risks Assessment data portal. Storm surge is derived from the ADCIRC Prediction System that solves for time dependent, circulation, and transport in multiple dimensions54. Maximum sustained hurricane wind speed was modeled hourly at a 5 km × 5 km resolution for 2017 by NASA’s Global Modeling and Assimilation Office (GMAO)55. The storm maximum wind speed for each 5 km × 5 km grid cell was calculated and binned into six discrete classes of wind speeds at 5 m s−1 increments: 26–30, 31–35, 36–40, 41–45, 46–50, and >50.Statistical analysesCanopy height losses measured from NASA G-LiHT data were grouped by five pre-storm canopy height classes (0–5 m, 5–10 m, 10–15 m, 15–20 m, and >20 m). All valid pixels within the lidar footprint was used to calculate the mean, standard error, and area (sum of 1 m × 1 m pixels) for each class (Supplementary Table 1). These results were then tested for significant differences between canopy height losses and pre-storm canopy height classes between using a one-way ANOVA analysis with a post-hoc Tukey test in R (version 4.0.3). For testing the significance between environmental variables (i.e., pre-storm canopy height, canopy height loss, percent canopy height loss, surface elevation, and storm surge water level above ground) we employed a two-sided Kolmogorov–Smirnov test56 implemented in R (version 4.0.3). First, we created a multi-band stacked image which included each of the variable layers. Within each resilience class (i.e., Low, Intermediate, and High) with randomly selected 10,000–20,000 points using Google Earth Engine to sample from the environmental variables images. From that sample set we then randomly selected 500 samples within each of the resilience classes. Each class combination (1) Low-Intermediate, (2) Low-High, and (3) Intermediate-High were compared using the Kolmogorov–Smirnov test. We repeated this procedure using 5000 iterations in order to provide a robust estimate of the Kolmogorow–Smirnov statistic, including the mean and first and third quartiles, which were then compared to the critical value (Supplementary Table 2).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
