The global flood protection savings provided by coral reefs

Coral reefs can provide significant coastal protection benefits to people and property. Here we show that the annual expected damages from flooding would double, and costs from frequent storms would triple without reefs. For 100-year storm events, flood damages would increase by 91% to $US 272 billion without reefs. The countries with the most to gain from reef management are Indonesia, Philippines, Malaysia, Mexico, and Cuba; annual expected flood savings exceed $400 M for each of these nations. Sea-level rise will increase flood risk, but substantial impacts could happen from reef loss alone without better near-term management. We provide a global, process-based valuation of an ecosystem service across an entire marine biome at (sub)national levels. These spatially explicit benefits inform critical risk and environmental management decisions, and the expected benefits can be directly considered by governments (e.g., national accounts, recovery plans) and businesses (e.g., insurance).

ranging from 0.1 to a coarse 25 km resolution. We chose different coastline resolutions on a case-by-case basis for the 39 sub-regions, depending on coastal sinuosity, representativeness of the coastal transects and computational efficiency. The intermediate resolution (1 km) was used in most of the sub-regions.
We divided the coastline in to 2-km cross-shore profiles globally. These profiles were defined orthogonal from the shore to deep water (e.g., see Supplementary Figure 5). As the initial points of all profiles had to be in deep water, the length of the profiles varied. We aggregated the 2-km coastal profiles into larger "study units" that were approximately 20 km wide (along the shoreline) for analysis and representation purposes. The study units were identified by: (1) a coastline segment 20 km wide, (2) a landward buffer parallel to the coast at variable distances (ranging from 1 to 20 km) and (3) a seaward buffer (parallel to the coast) at variable distances (from 30 to 100 km). Irregularities in the coastline such as estuaries and bays resulted in study units that were sometimes wider than 20 km. Supplementary Figure 6 shows examples of these study units in the Caribbean.
Offshore hydrodynamics. Hydrodynamic conditions at offshore locations (Supplementary Figure 3) are obtained by combining global data on wave parameters, mean sea level (MSL), astronomic tide (AT), and storm surge (SS) in to 280,320 1-hourly sea states covering the period 1979-2010. Each single dataset has been used previously in global climatology studies and validated with buoy measurements, satellite altimetry and tidal gauges 2,3 . To propagate the offshore data to the nearshore and calculate flooding, the 32-year long time series of offshore wave and sea level data are classified into a set of 500 sea state conditions using prior published, data-mining based methods 13,20 . Nearshore hydrodynamics data. We used the global ETOPO bathymetry 26 and combined it with the SeaWiFS (Sea-Viewing Wide Field-of-View Sensor) bathymetry for coral reefs 11,12 . SeaWiFS reef bathymetry data has a 1 km spatial resolution and was collected over the period 1997-2002. The current SeaWiFS-derived depth map, updated on 6 November 2002, represents the contribution of 28,937 Global Area Coverage (GAC) resolution files. Each SeaWiFS map pixel (0.01 by 0.01 degrees) averages data from no fewer than four input pixels, and many had more than 200 input pixels.
Nearshore hydrodynamics reef wave model. The reef model assumes a longshore uniform, two-dimensional vertical (2DV) coordinate system, with waves propagating shoreward. Wave propagation over the reef is calculated from linear wave theory. Wave propagation is modeled at shore-perpendicular one-dimensional transects, therefore along shore processes such as longshore currents are neglected. The evolution of a wavefield of root-mean square (rms) wave height H with weak mean currents is computed by solving the wave energy balance equation: where Ew is the wave energy density and Cg the group velocity. The dissipation of wave energy flux is caused by wave breaking (Db), bottom friction (Df), and the presence of vegetation in the water column (Dv ), which is not considered in this study. Equation (1) is widely applied in coastal studies to assess wave propagation (e.g. SWAN) 27 and previously applied to reef environments 28 . As applied here, Db and Df are expressed following Thornton and Guza 29 : where  is water density, g the constant of gravity, k the wave number,  the angular wave frequency and f p the peak frequency. As the incident waves propagate into shallow water on the reef flat, they also are increasingly affected by frictional effects. On some barrier reefs, friction has been shown to be the dominant dissipative process 31 . Wave setup depends strongly on the profile geometry and reef architecture. To define the bottom friction, we build from specific work by Sheppard and others 32 and Nunes and Pawlak 33 on coral reefs and wave damping. Shepard and others 32 have identified friction factors (f w ) based on live coral cover as summarized in Table S3. We have used these recommendations to identify different friction factors by oceanic region (Supplementary Figure 4) based on overall coral condition 34 . We assume the greatest live coral reef cover in the Pacific Islands and Micronesia (f w =0.2); lower live coral reef cover in the Indian and west Pacific Oceans (f w =0.16) and lowest overall live coral reef cover in the Caribbean and Latin America (f w =0.14).
Nearshore hydrodynamics total water level model. The total water level (i.e., flood height) along shorelines is a function of mean sea level, astronomical tide, storm surge, and the run-up of waves 35 .The run-up represents the wave-induced motion of the water's edge across the shoreline and is built of two contributions, namely the wave setup at the shoreline and the swash representing oscillations about the setup. The run-up calculation requires obtaining the local wave conditions at the shoreline using the reef wave model above.
Nearshore hydrodynamics computation of wave setup. Wave set-up is a significant contributor to water levels, particularly in coral reefs environments, and hence is an important factor in determining flooding. The wave-setup is obtained from the conservation of mass and the momentum equations 36 . In our one-dimensional setting, the computation of the wave-induce setup is based on the vertically integrated momentum balance equation 37 . Similar implementations have been used in previous work to evaluate the effect of vegetation on waveinduced setup 7 and in coral reef environments 38 . Comparison of natural forms of reef profiles and idealized horizontal reefs show that the relative setup on natural reef profiles is found to be less than the idealized ones 39,40 , partly because they do not consider longshore effects.
where H 0 is the offshore significant wave height, L 0 represents the deep-water wave length, T p the peak period, and m the bathymetry slope in the foreshore beach slope. This equation expresses runup as a function of empirical estimates of incident wave setup at the shoreline (first term of the equation) and the swash incident and infragravity band frequency band components (second term of the equation).
The wave setup is calculated from the propagation of waves in the reef model, as described above. For the swash and infragravity band frequency components, we apply the second term of the equation to find their contribution to the total water level as follows:  Using the wave propagation model over the reef, we calculate the breaking point position (x b ), breaking depth (h b ) and breaking height (H b ). The breaking point depends on the incident wave conditions, the sea level, the reef geometry and friction and is calculated from D b (equation 2).  Following Stockdon et al. recommendations we deshoal the wave height to deep water to obtain H 0 L 0 is calculated for the corresponding peak period (T p ).  The foreshore slope is obtained for each of the different profiles from the DIVA-GIS dataset (http://www.diva-gis.org/).
In our approach, we assume that Stockdon et al. 41 can be applied to coral reefs as the model was developed to include barred beaches, which resemble coral reef protected beaches. Modifications of the same formula have been applied previously to estimate the effect of vegetated ecosystems on runup 7 . In reef environments, wave runup depends strongly on the reef complexity and architecture 42,43 . An alternative to this approach to calculate the runup swash component would be the application of nonlinear process based models to coral reef environments [43][44][45] . However these nonlinear models are computationally demanding and require high resolution bathymetry. Therefore, they are only feasible for small spatial domains and a very limited set of sea states. For validation purposes, we have compared the runup predicted by our reef model versus XBeach 42 for an idealized fringing reef profile. XBeach is a state-of-art hydrodynamic model that has been applied in reef environments (see section below Validation of flooding model).
Extreme water levels and flood height reconstruction. From the propagations of waves and the calculation of total water levels onshore, the reconstruction of the flood height time series at the most onshore points is based on multi-dimensional interpolation techniques 13 . We apply a Peak Over Threshold method to select extreme flood heights and fit a General Extreme Value Distribution 15 to obtain the flood heights associated with the 10-, 25-, 50-and 100-year return periods. The methodology has been tested in case studies 46,47 and validated with observations showing a good reproduction of time series structure and statistical parameters in shallow water. The methodology is not location dependent and is applicable globally.

Estimating reef benefits.
To examine the current value of reefs for coastal protection, we compared flooding under current conditions, "with reef", to the flooding in a scenario "without reefs". In our "without reefs" scenario, we do not assume the loss of the entire reef habitat; we assume only the loss of the top 1m in height across the reef bathymetric profile. Many ecosystem service assessments assume the entire loss of a habitat for estimating benefits. For example, the replacement cost method, which is the most commonly used method for estimating the benefits from mangrove and reef habitats 48 , identifies the flood reduction benefits from habitats by estimating the cost of replacing them with seawalls or breakwaters. Many problems have been identified with the replacement method and it provides estimates of values 10 times higher than the recommended Expected Damage Function approach that we follow 48,49 .
The "without reefs' scenario is not meant to be a prediction of site-specific trajectories for reefs, but nonetheless this level of loss is already observed to be happening 6,50,51 and is conservative relative to future predictions of reef loss [52][53][54] . In addition to the widely observed declines in coral cover, growth and condition, all of which affect reef height 55,56 , new measures of seafloor elevation show that bioerosion and carbonate dissolution are degrading height across all reef habitats including on reef flats 51 . Damage from storm events can also create losses in reef height of 1-3 m 57,58 and can devastate whole shallow reef frameworks 59 . Past storms have removed many branching and massive corals at the shallowest depths 57,60 . Shallow corals have evolved with intermittent storms and can recover from them, but this is more difficult when reefs are exposed to multiple stressors 24 .
Calculating people and assets flooded. We assessed flood heights along each coastal profile and then identified the area flooded within each coastal study unit. We extended the flood heights inland by ensuring hydraulic connectivity between points at a 90m resolution; a significant advance over more common bathtub approaches in earlier global flooding models. From the flooding levels and flooding extent, we calculated the total area of land affected and damages at each study unit. Flooding maps were also intersected with population data 17 after resampling from the original 1 km resolution to the 90 m of the digital elevation model. Existing artificial defenses such as seawalls were not assessed, because data on defenses only exist for a very few areas globally; these built defenses are also less common in tropical, developing nations.
We expanded on earlier approaches to infer built capital from population data 1,18 by identifying the ratio between built capital per capita and the Gross Domestic Product (GDP) per capita for each country 1 in 2011 US$ using information from the World Bank 61 . We filled data gaps for several countries by using the average from countries with similar income levels and affiliation to the Organization for Economic Cooperation and Development (OECD). The overall global mean ratio that we obtained (2.67) is similar to that obtained by Hallegatte and others (2.8) 18 . However, we did identify significant differences in the ratios across some countries and regions (e.g., Cuba -4.53, Vietnam -3.22, Australia -3.17, Philippines -2.68, United Arab Emirates -1.98, Micronesia -1.38).
Assessing damages and estimating annual benefits. We followed existing approaches for assessing the damages to built capital as a function of the flooding level 62 . We calculated the percentage of built capital that has been damaged (D) for a given flooding level and a certain coefficient that must be calibrated as D(h) = h/(h +k). This curve indicates that as flooding level increases, the percent of damages to built capital also increases. While there is debate about the right k to choose, we have followed others in using k = 1.0, which means that the built capital flooded at 1m of depth loses 50% of its value 62 . We follow standard terminology where the total built capital flooded is the exposure of assets and the value lost is the damages. The economic benefits of flood protection are the avoided damages.
In addition to assessing risk and damages for particular events (e.g., 100-year storm event), we also examined average annual expected loss 63 . To estimate annual risk, we integrated the values under the curve that compares built capital damaged by storm return period, i.e., the integration of the expected damage by the probability of the storm events 62 .
Sea level rise. We used regional sea level rise projections and associated uncertainty estimates produced Slangen, et al. 19 on a 1 o x1 o grid. Projections account for changing ocean circulation, increased heat uptake and atmospheric pressure in 21 CMIP5 climate models 64 , combined with model-and observation-based regional contributions of land ice, groundwater depletion and glacial isostatic adjustment, including gravitational effects due to mass redistribution. In our analysis we select the estimates for the Representative Concentration Pathway 8.5, which corresponds to a global mean surface temperature of 2.7-5.4 o C and a global mean sea level rise of 0.71  0.28 m (from 1986-2005 to 2081-2100). Regional variations deviate up to 20% higher from the global mean in the subtropical and equatorial regions. We calculate the regional sea level rise at the end of the century, at each profile, by storm return period and add these values to the estimates of the total water level offshore. The wave propagations over the reefs are recalculated and the flood height statistics reconstructed as described above.
Sensitivity analyses. Hinkel, et al. 62 carried out a global analysis of coastal flood damage under sea level rise. Our two approaches are in general very similar with some differences in the global data sets selected and the fact that we use a process-based flooding model. They analyzed the sensitivity of the flooding models to bathymetry and topography and found that corrections needed to reduce uncertainties are feasible at local scales but cannot currently be applied at global scale for technical constraints. This conclusion is applicable to our analysis.
To identify the sensitivity of our flood height estimates to the factors affecting wave run-up, we examined two major reef types: barrier and fringing reefs. From our global set, we selected 1000 profiles for each reef type that were representative of their geographic variations (2000 profiles total). We examined the effects of reef friction (C f_reef ), wave breaking (γ coral ), and water depth (h coral ) on wave run-up across the 2000 profiles. We varied the reef friction parameter in increasing steps of 0.01, from 0.08-0.20, as taken from Sheppard and others 32 , which cover the range of friction factors from sand to 75%-100% live coral. We varied the wave breaking parameter proposed by Yao et al. (2013) for γ coral = 0.2 to 1.2, in increasing steps of 0.1. We varied water depth above the reef in increasing steps of 0.1 m, from 0.1 m to 1 m.
Supplementary Table 4 summarizes the results of the sensitivity analysis. For each of the reef types, we show the maximum, minimum and mean run-up variation observed across all the incremental changes in the parameters (e.g., for the water depth analysis, the table, first column, summarizes the results from 1,000 profiles at each of 10 water depths for a total of 10,000 runs for the barrier reef type).
For the reef friction factor, the results indicate that the overall effects of changes in friction on run-up are very low, with a 13% difference across the entire range of factors examined (a ~1% change in run-up for each increment in the friction factor that we tested).
For the wave breaking coefficient, results show that run-up varies on average from 1.71%-2.53% per increment. Based on the formulation by Yao et al. (2013), and considering a variation in the breaking coefficient from 0.2 and 0.6, our results indicate a total variation in results from 6% to 10% (Supplementary Table 4).
Changes in water depth resulted in slightly larger variations in run up in our sensitivity analysis. On average, and depending on the type of reef, every 10 cm increase in water depth results in an increase of 1.68%-2.75% in run-up. This is mainly due to the effect that depth has on breaking or shoaling conditions, and clearly indicates the importance of an accurate bathymetry (Supplementary Table 4).
In sum for the sensitivity analysis, we conclude that the range of friction coefficients recommended by Sheppard and others 32 produces a variation of 13% in wave run-up, whereas variations in the wave breaking coefficient recommended by Yao et al. (2013) introduces variations up to 10%. Changes in water depth may change run-up from 16% to 27% over a 1 m range. Hence it can be concluded that the most relevant aspect to reduce uncertainties in the results is using high quality bathymetry with a good definition of coral reef height. The importance of bathymetry is relevant for all coastal flooding models.
Validation of the wave and flooding model. Our physical modeling approach relies on linear wave theory to calculate the effect of reefs on waves and wave setup. We use a modification of an empirical formula to estimate the swash component of run-up (see above) that was originally calculated for beach environments. Models based on the same equation (e.g., SWAN) have been used to study the effect of reefs on flooding at regional scales 28,65,66 . More complex but still 1D models have been used for the study of the hydrodynamics of coral reefs in more detail at specific sites and they also use transect-based approaches 42 .
We compare and validate our modeling approach against a more complex model for reef environments, XBeach 43,67 in a theoretical fringing reef profile. The comparison was performed in a frictionless setting to focus on comparisons of wave breaking and run-up. XBeach was originally derived from sandy beaches and has hitherto been successfully applied to predict erosion and overwash under hurricane forcing 68 and more recently, successfully applied in reef environments 42,43 .
We considered three different reef depths or sea levels (0, 1.5 and 3m) and three different offshore significant wave heights (1, 3 and 6m) and estimated flood height at the shore with each model. The comparison shows that the flood height values of our approach are comparable to those from a higher resolution model (Supplementary Figure 7, r 2 =0.978) given the simple geometries that are characteristic of the global bathymetric data.