Methods at a glance

This work measures the flood protection service of mangroves all over the world for two climatic conditions: (1) Cyclonic- i.e., the conditions high-intensity extreme waves and storm surge induced by tropical cyclones and (2) Non-cyclonic, i.e., the “regular” waves generated by low-intensity local storms. We followed the Expected Damage Function (EDF) approach50, recommended by the World Bank36, previously applied in coral reefs ecosystems34 and commonly used in engineering and insurance sectors51,52. We examine the role of mangroves in reducing flood risks by measuring the impacts of flooding on people and property under two different scenarios: with and without mangroves. The “without mangroves” scenario assumes the complete loss of mangrove habitat and the consequent erosion of the intertidal area with a smoothened bottom roughness. We use a regression model globally, to calculate coastal flooding by analyzing more than 7,000 historical cyclones53 and 32 years of regular waves and sea level (storm surge, astronomical tide and mean sea level). Flood impacts (i.e. land flooded) for the with and without mangrove scenarios are combined with global distributions of people and property11, and with vulnerability based on global “Flood Depth-Damage Functions”54 to assess baseline flood damages and flood damages after mangrove loss for multiple storm events and on an annual basis.

To identify the mangroves that influence a given coastal region and evaluating nearshore hydrodynamics and flood height we define cross-shore coastal profiles (Supplementary Fig. 6d). Then, we follow a multi-step framework whose key aspects are described here and in the Supplementary Material: (1) Estimate offshore dynamics produced from both, tropical cyclones and regular climate conditions. (2) Estimate nearshore dynamics by downscaling offshore waves and storm surge until shallow water, just before mangrove habitats. (3) Propagate waves and storm surge through mangroves and obtain the flood height behind the mangroves at the shoreward end of each profile. (4) Estimate the land flooded (impact) due to extreme water levels along the shore by intersecting the flood height at the shoreline with inland topography (5) Calculate land, people and property located in the flooded area and, finally, apply the corresponding damage functions to obtain flood damages with and without mangroves.

Study domain description

This global study covers 700,000 km of coastline that includes more than 141,000 km2 of mangroves, spread over 4 continents and more than 9,500 islands. To reduce the vast computational requirements such a large domain requires, the global domain is divided at three levels, from global to regional and local (Supplementary Fig. 6). The first level is a global division into six macro-regions corresponding to the following ocean basins of tropical cyclone generation53: East Pacific, North Atlantic, North Indian, South Indian, West Pacific and South Pacific (Supplementary Fig. 6, panel “a”). The second level divides the 700,000 km of global mangrove coastline into 68 sub-regions considering coastline transects of similar coastal typology (e.g. islands or continental coasts) and similar ecosystem characteristics (Supplementary Fig. 6, panel “b”). The third level of disaggregation is done at a local scale, defining units with 20-km length of coastline and extending up to 30-km inland and 10 km seaward (Supplementary Fig. 6, panel “c”). Within these units cross-shore profiles perpendicular to the mangrove habitats are created for each kilometer of mangrove coastline, totaling 700,000 profiles (Supplementary Fig. 6, panel “d”).

Building the global model based on the Philippines results

Global reanalysis of ocean and coastal waves55,56 and storm surge57 exist though not for tropical cyclones during the course of this work. We develop the global model based on an extensive re-analysis of tropical cyclone climates developed for the Philippines29. We chose the Philippines as the baseline case to develop our own tropical cyclone reanalysis at high resolution and estimate flood damages in presence and absence of mangroves. There are three main reasons that make the Philippines an excellent pilot case for valuation of the coastal protection ecosystem service provided by mangroves: (i) Almost 10% (548 events) of the global tropical cyclone records from IBTrACS database affected the Philippines53. The worldwide distribution of tropical cyclone parameters (velocity track, wind speed) closely resemble these events in the Philippines (Supplementary Fig. 7) (ii) The islands of the Philippines present high climatic variability and it is at particularly risk from natural hazards like typhoons and regular storms, which are the cause of 80% of the total losses from disasters [average loss totaling nearly $US 3 billion, 29% of this damage is due to coastal flooding58]. (iii) The Philippines ranks in the top 15 most mangrove habitat-rich countries, with 2,630 \(k{m}^{2}\) in 2010, representing 2% of the world total59. These mangrove habitats show extensive variation in both cross-shore length and average depth. Mangroves in the Philippines rage between 0.1 km and 8 km wide and between 0 m and 10 m depth (Supplementary Figs. 8 and 9). We valued flood protection service of mangroves in the Philippines by using the numerical model Delft3D considering both historical tropical cyclones and regular climate conditions with and without mangroves. We use these results to build two global statistical models. The first global model is developed to obtain offshore and nearshore ocean dynamics produced by tropical cyclones (wave height, peak period, storm surge and storm duration), and the second global model to estimate how the presence and profile of mangrove habitats influence the total water level on the shoreline. Further details of the two models are developed below:

Model 1: Offshore and nearshore dynamics generated by tropical cyclones

Offshore waves and storm surge generated by tropical cyclones (IBTrACS database) were numerically simulated in the Philippines by using Delft3D modules “Flow”60 and “Wave”61. Both modules were run simultaneously in a 2-dimensional grid of 5 km cell-size with a time step of 30 s, forced with hourly wind data and sea level pressure fields obtained from parametric model, in which the non-linear interaction processes of tide, wind setup, inverse barometers and wave setup are considered. The model was validated by comparing the storm surge generated by typhoon Rammasun, in Legaspi and Subic Bay. We use tidal gauges registers from the Global Sea Level Observing System (GLOSS, http://www.gloss-sealevel.org) for validation (Supplementary Fig. 11). Using the results of the numerical simulations carried out with the Delft3D model in the Philippines we look for statistical relationships between cyclone parameters and oceanographic variables to create a new predictive model, where oceanographic variables (wave height, period, weather tide and duration of the storm peak) are predicted based on cyclone parameters (distance, wind speed, track velocity, wind angle of incidence). In the Philippines, 548 events were simulated on a two-dimensional grid of 5 km cell size, finally creating a database of 58 million results. We randomly select 90% of the generated results to build our predictive model, and use the other 10% for validating the predictive models. We estimate the correlation between physical tropical cyclone parameters (distance from the trace to the profile D [km], wind speed W [km/h], cyclone travel speed V [km/h], wind direction from north θ WN , [in degrees] and the angle between the wind direction-profile θ WP [in degrees]) and the oceanographic variables at the target point (maximum significant wave height produced during the event at the target point H smax [in m], peak period Tp [in s], maximum storm surge SS max [in m] and maximum storm surge duration, T SSmax ). We increase the accuracy of the analysis by dividing the data into two groups: Coastal areas directly exposed to tropical cyclones and areas protected from the direct impact of tropical cyclones (Supplementary Fig. 13). For each combination (5 cyclone variables x 3-time instants x 4 oceanographic variables = 60 cases), we estimate the Pearson coefficient (P xy ), which statistically quantifies the degree of correlation between the cyclone variables “X” and the oceanographic variables “Y” (equation S.1). We then adjust ocean climate variables (Yi) to our parametric model [equation S.2]. We test this adjustment for one, two, three and four independent variables (Xi), so that we can cover all the alternatives and, based on the correlation coefficient of each one, choose the best regression model.

$${{\rm{Y}}}_{{\rm{s}}}={{\rm{a}}}_{0}+{{\rm{a}}}_{1}\cdot {{{\rm{X}}}_{1}}^{{\rm{\alpha }}1}+{{\rm{a}}}_{2}\cdot {{{\rm{X}}}_{2}}^{{\rm{\alpha }}2}+\ldots +{{\rm{a}}}_{{\rm{n}}}\cdot {{{\rm{X}}}_{{\rm{n}}}}^{{\rm{\alpha }}{\rm{n}}}$$ (1)

Where Y could be either the maximum wave height (Hs max ), the peak period (Tp), the maximum meteorological tide (SS max ) and the duration of the meteorological peak produced by the cyclone (Tss max ). Meanwhile, X could be any of the following predictor variables (see equations S.3 to S.10): minimum distance between the storm track and the target point (D min ), the wind speed when the tropical cyclone is at the closest location to the target point (W dist_min ), the average wind speed during along the storm length (W mean ), the mean direction of wind respect to the North (θ WN_mean ), the wind direction respect to the North at the minimum distance point (θ WN_dist_min ), the average track velocity (V mean ), the track velocity at the minimum distance point (V dist_min ) and is the track velocity at the moment of maximum wind speed (V wind_max ).

Model 2: The role of coastal habitats in nearshore dynamics (flood height)

Coastal vegetation provides resistance to the energy and flow of waves and water as they come onshore which is modeled by using a friction factor. Mangroves are then modeled in terms of an equivalent roughness [e.g. Sheppard friction for coral reefs62] based on Manning coefficient. In the Philippines we classify surface types into three groups: sandy soil (n = 0.02)63, mangroves (n = 0.14)63 and coral reefs (n = 0.05)64. 1-dimensional numerical propagations are carried out using the Delft3D model to obtain flood heights along the coast (Supplementary Figs. 8 and 9). We use these numerical results to create two interpolation tables (for both, regular climate and tropical cyclones) that correlate the climatic information at seaward side of the profile (Hs, Tp, SS and Tss, being this last only specific of tropical cyclones) and the characteristics of the mangrove profiles (width and average depth) with the flood height (i.e. total water level along the coast). These tables contain 37,500 tropical cyclone simulations (50 cyclones × 750 profiles) and 90,000 regular climate simulations (120 sea states × 750 profiles).

These two models described above are integrated in the multi-step framework applied globally for valuing the flood protection role of mangroves:

Step 1: Offshore dynamics

The offshore hydrodynamic conditions (wave height, wave period, storm surge and astronomical tide) were subdivided in two groups: (1) those produced by local extreme events (tropical cyclones) and (2) those produced by less intense local climate or extreme climate generated far away from the study area (regular climate). Regular climate is defined by different datasets within the period 1979–2010: a global wave reanalysis56, a global storm surge reanalysis57, astronomical tide65,66 and mean sea level compiled from historical numerical reconstruction and satellite altimetry67. Waves and sea level conditions due to tropical cyclones are excluded and studied separately to avoid double counting. Tropical cyclones were considered separately from regular climate only if two conditions are satisfied: (1) they are generated within the same ocean basin than the study area and the cyclone passes closer than 500 km from the coastline, and (2) 10-minute sustained wind speeds (W 10m ) exceed 118 km/h. Tropical depressions (W 10m ≤ 62 km/h) and tropical storms (63 km/h ≤ W 10m ≤ 118 km/h) are studied together with regular climate. For historical tropical cyclones, we used IBTrACS database53, which provides 6-hourly data of wind speed, atmospheric pressure and position. Since global reanalysis of tropical cyclones that include waves do not exist, we use a statistical model (Equations S.3 to S.10) created from the Philippines results to calculate offshore wave height, peak period, storm surge and storm surge duration just in the limit between deep and shallow water.

Step 2: Nearshore dynamics

Once we resolve offshore dynamics, we obtain waves and storm surge in the seaward side of each cross-shore profile. Waves interact with the bottom and other obstacles (e.g. islands) as they approach the coast and modify height and direction through shoaling, refraction, diffraction and breaking processes. Regular climate is propagated following a hybrid downscaling. The 32-year long series, from 1979 to 2010, include 280,000 sea states (1 sea state is 1-hour register of wave height, peak period and total water level). To each profile, we allocate the closest point of the offshore databases. Considering all, the 700,000 coastal profiles and the 280,000 sea states results in an unmanageable number of cases. We reduce the number of sea-state propagations by, firstly, considering only the 3,787 non-repeated combinations of wave height, peak period and total water level (SS + AT + MSL) and, then, applying The Maximum Dissimilarity Algorithm (MDA68,69) to finally obtain 120 sea states to be propagated with Snell law and shoaling equation (Eqs. 1 and 2). Meanwhile, tropical cyclone nearshore hydrodynamics are obtained by means of the previously derived regression model (equations S.3 to S.10). We apply regression models in each profile, and we obtain the same parameters as for regular climate, in addition to the time duration of the meteorological tide (T ss ).

$${H}_{s\_perfil}={H}_{s\_off}\cdot \sqrt{\frac{{C}_{off}}{{C}_{perfil}}}\cdot \sqrt{\frac{\cos \,{\theta }_{off}}{\cos \,{\theta }_{perfil}}}$$ (2)

$${\rm{Snell}}:\frac{{C}_{off}}{\cos \,{\theta }_{off}}=\frac{{C}_{perfil}}{\cos \,{\theta }_{perfil}}$$ (3)

Step 3: Modeling the role of coastal habitats in nearshore dynamics, flood height

The next step consists on propagating ocean hydrodynamics over mangroves forest which dissipate wave and surge energy, and, consequently, reduce flood height. Flood height is a function of mean sea level, astronomical tide and run-up of waves. Mangroves dissipation takes place by means of breaking and friction processes. Given the large scale of this global analysis, we follow a simplified approach for vegetation modeling. We use the interpolation table from the Philippines to infer the resulting flood height given mangrove length and depth, significant wave height, peak period and total water level at the head of each cross-shore profile. Then, we apply the statistical reconstruction technique RBF (Radial Basis Functions)70 to calculate in each profile the complete historical flood height time series. Next, we carry out an extreme value analysis. First, we select maximum values on a variable threshold (minimum, 1-in-5-year event). We adjust these selected values to a Generalized Pareto-Poisson distribution, and we obtain the flood height vs return period curves for both scenarios: with and without mangroves. We observe a high spatial variability of flood height produced by tropical cyclones along worldwide coastlines, which highlights the importance of addressing global flood risk analysis at high resolutions to consider local topographic and bathymetric variations (e.g. 1-in-100-year flood in Vietnam, Supplementary Fig. 10). We assume that countries with less than 100 ha of mangroves were excluded from the analyses as there were too few mangroves to reliably estimate benefits using a global model. This excluded 15 countries in total, including Bahrain and Benin, which had some of the most over-estimated values of benefit/ha; as well as eight Caribbean Small Island Developing States (Supplementary Table 9).

Step 4: Calculating impacts: flooding maps

Uncoupling wave and surge propagation from the flooding process allows us to freely choose the most accurate flooding method. Other alternative strategies exist, but they are unapproachable at global scale due to its computational cost and high-quality data required, usually unavailable at global scale (e.g. using coupled phase resolving models or phase averaged models). To obtain flood maps by means of uncoupling propagation and flooding processes require: (i) the flood height along the coast with high enough resolution to avoid significant longshore gradients, (ii) a good DTM (Digital Terrain Model) and (iii) a flooding method (flood models). Separating flooding process from waves and sea level propagation gives us more flexibility to adapt the flooding approach to the elevation data. Local scale analyses (<100 km of coastline) with high resolution DTM (<10 m) could be addressed by using process-based flood models, like RFSM-EDA (Rapid Flood Spreading Method – Explicit Diffusion wave with Acceleration term)71,72. However, larger scale´s (>100 km) with coarser DTM (>10 m), require fast and less precise techniques, such as “bathtub” method, based on hydraulic connectivity which consists of merge points below the flood height. We use this last strategy to address global flooding in presence and absence of mangroves. The flood extent is estimated globally by using a 30-meter SRTM-DTM (Shuttle Radar Topography Mission)73.

Step 5: Assessing global flood consequences in mangroves protected areas

Mangrove benefits are assessed in terms of avoided damages to people and property. Property value is directly obtained as the sum of industrial and residential stock from GAR15, at 5 km resolution worldwide11. The suitability of this database to be used in global assessments of coastal flooding exposure and damage lies on the fact that it integrates homogeneous global population and country‐specific building typology, use and value data74. The consistency of the methodological approach used in the development of GAR15, as well as the choice of the best data currently available for its implementation, have produced a product fully adapted to the needs of the global model of the evaluation of probabilistic risk74. Consequently, GAR15 is the most appropriate source of data available for global scale analysis, looking for an order of magnitude of the value of adequate protection for mangroves, usually used by critical stakeholders such as the World Bank75. People distribution comes from the freely available 1 km resolution database GPW (Gridded Population of the World), from SEDAC (Socioeconomic Data & Applications Center). To be consistent with flood layer grid resolution (30-m) it is necessary to redistribute people and property over a finer mesh. We apply a downscaling method, in which the re-distributed values are calibrated with the other existing data of people distribution (WorldPop) by imposing as boundary condition that the total sum of the assets re-scaled to 30 meters in each region was equal to the sum of those same assets without re-scaling in the same control zone. The sensitivity of people and stock to different levels of flooding is obtained through different damage functions. Damage functions provide information of the number of people affected by coastal flooding and the stock losses, according to the water depth. We use different damage functions for population and for stock. Population damage is based on the hypothesis that water depths below 0.5 meters do not affect people, while water depths above 0.5 meters affect 100% of people hit by flooding. It is a common practice in the scientific literature not to use damage functions to calculate the population affected by floods3. This option overestimates the results obtained; therefore, it is recommended to opt for a certain threshold below which the effects of flooding are not considered54. This threshold is set at 0.5 meters because it is a common value used by emergency services (Japan, Netherlands, USA) in determining whether or not it is necessary to evacuate people from an area under threat. In case of stock, we adapted the global flood depth-damage functions from Huizinga/JRC (Joint Research Centre) broken down by continent (Africa, Asia, Oceania, North America, South America and Central America) and by asset type: residential and industrial54. Average values of damage at different water depths are provided in Supplementary Table 5. Finally, flood risk (magnitude and probability) is obtained by combining damage curves with people and property exposure distribution. Then, we integrate the return period curves to obtain the Expected Annual Damages and Benefits at each 20-km study unit. We can thus show global information on annual flood damage anywhere and with a spatial resolution high enough to be incorporated into coastal planning and ecosystem conservation policies.

Modeling assumptions

To provide a more nuance discussion of the strengths and weaknesses of the modelling approach, we have included a table with all the assumptions stablished at each step of the methodology, as well as the corresponding reference to the existing literature where this assumption is applied and validated (Supplementary Table 9). This table summarizes the assumptions considered in this work and may help the reader to assess how strong the assumptions are and potentially identify areas for future work.