The role of climate and population change in global flood exposure and vulnerability

Flood modeling

We have developed a global flood model of flood depth inundation throughout the century for coastal ocean and lake, fluvial, and pluvial flood hazards for various years (baseline, 2020 to 2100 at 5-year increments), return periods (10 to 500 years), uncertainty levels (5th, mean, 95th), and climate scenarios (baseline, SSP1-2.6, SSP2-4.5, SSP5-8.5). The digital elevation model (DEM) used to define the Earth’s surface was based on a wide variety of publicly available sources, including MERIT19, Coastal90 dem20, USGS NED33, and other publicly available datasets in Europe, Canada, Japan, and Australia. The model has global coverage and is based on a 90-m grid scale using the EGM96 geoid. Some areas have enhanced modeling at 10-m resolution. We employed separate models for coastal, fluvial, and pluvial flooding, and then combined the results.

Coastal ocean and lake model

The coastal ocean and lake flood model used multiple climate projection datasets to estimate the effects of sea-level rise and storm surge, tides, and waves on coastal inundation, as well as storm surge and lake levels on lake shoreline inundation. The basis of the model is at coastal points globally (see Supplementary Information Fig. SI8), whereby a series of sequential operations were conducted to compute the return-period-based water levels for a range of years and climate scenarios. Data inputs for the coastal ocean and lakes are from separate sources at the coastal points.

Coastal Ocean

Beginning with a storm surge model (GTSR)21 which includes tidal effects, the extreme values statistical parameters were debiased using ocean observations (GESLA-3)34 at 14,979 coastal ocean points globally. Next, we included the effects of parameterized wave setup22 using validated data from the WaveWatch III Hindcast35. We then included relative sea level rise projections from the IPCC AR6 projections28 using nearest neighbor interpolation (which includes estimates of local land movement, glacial isostatic adjustment, thermal expansion, water storage, etc.).

Coastal lakes

Lake data was obtained from long-term remote sensing (G-REALM)36, and in-situ observations (NOAA, CO-OPS water level observation data) for 656 lakes globally. The short-term variability in historical lake levels was used to define extreme value statistical parameters. While different processes dominate in different coastal lake settings, our general approach is data-driven, extreme value analysis. The long-term lake levels follow the recent trends since 1950 in the near term but are gradually reduced to zero trends towards the end of the century. This is to account for uncertainty in forecast trends using historical data projections, which can arise from the climate but also human water use. Thus, long-term lake levels were based on historical trends without climate scenarios, which is an avenue for future development.

Combined coastal ocean and lake model

The statistical parameters used to define short-term extreme values, and the long-term water level trends for coastal oceans and lakes were then combined using a joint uncertainty analysis and a Monte Carlo sampling approach.

Wind projections were based on extreme value analysis of downscaled daily maximum wind speed from multiple CMIP6 GCMs37,38, specifically with the following models (number of ensemble members indicated in parentheses): ACCESS-ESM1-5 (3), CESM-LENS (40), CESM Med. Ens. (15), CNRM-CM6-1 (6), MIROC6 (50), MPI-ESM1-2 (10), and UKESM1-0 (6). In tropical cyclone (TC)-prone regions, wind projections were supplemented using a very large sample of synthetic TC tracks39,40 with a parametric TC intensity model adjusted for future climate conditions41,42,43. Parameterized climate effects of changing winds on short-term extreme water level values were applied to both coastal ocean and lake points. We assume wind energy is the primary driver of short-term water level fluctuations arising from a variety of mechanisms including wind-driven storm surge, wave setup, seiches, etc. This assumption may under- or over-estimate changes in extreme values locally, however, due to weak climate trends for future wind conditions, the impact on the global 100-year coastal water levels is small with no long-term trend and less than 10% variation. In tropical cyclone-prone areas there is a 4% average increase by 2100 and less than 20% variation (see Supplementary Information Fig. SI9).

Finally, datums were adjusted from mean sea level to EGM96 framework using the Mean Dynamic Topography44. These coastal water levels were then passed to the inundation model to compute depths of flooding on the grid (see below).

Fluvial Model

The fluvial flood model is a cascading sequence of hydrologic and hydraulic simulators that derive depths within a river network by routing daily runoff from a suite of 27 simulations of the CMIP6 suite of GCMs37,38 downscaled to one km2 resolution. The depths within the river network were then passed to the HAND inundation model to compute the depths of flooding on the grid (see below).

Specifically, we used the following models (number of ensemble members indicated in parentheses): ACCESS-ESM1-5 (6), BCC-CSM2-MR (1), CESM2 (1), CMCC-ESM2 (1), EC-Earth3 (1), HadGEM3-GC31-LL (1), INM-CM5-0 (1), IPSL-CM6A-LR (2), KACE-1-0-g (1), MIROC6 (3), MPI-ESM1-2-HR (1), MRI-ESM2-0 (5), NorESM2-MM (1), and UKESM1-0-LL (5). Daily runoff is defined as the excess precipitation and snowmelt available for surface and subsurface runoff after interception, infiltration, ponding, and evapotranspiration losses are subtracted. Daily precipitation runoff was spatially averaged over hydrologic response units based on sub-catchments defined by a vectorized global river network dataset extracted from high-resolution 90-m datasets. The time series of spatially averaged runoff was converted into a time series of lateral flows into the river network using a convolution with a triangular unit hydrograph adapted to local watershed characteristics. Daily time series of lateral inflows were routed through the river network using the Muskingum routing method to obtain daily stream flows. Finally, the fluvial process employed a 1D river channel model to solve for the depth, using a hydraulic model based on a quasi-steady state form of the St. Venant equations, local river geometry cross-sections at 90 m spacing, Manning’s roughness, and stream flows. The quasi-steady state form drops the local acceleration from the inertial terms but keeps the convective acceleration term. Backwater effects from boundary conditions at the ocean were considered at high tide. A maximum annual series of stream flows was used to inform the extreme value analysis to determine return periods and associated uncertainty. Channel flow depths associated with maximum stream flow at the desired return periods were used in the inundation algorithm to fill the floodplain.

Floodplain inundation for fluvial and coastal

The floodplain inundation model used for both fluvial and coastal floods was based on the Height Above Nearest Drainage (HAND) algorithm45, which applies flood levels of constant elevation along flow paths. For fluvial floods, we applied the HAND algorithm to simulate propagation away from river channels. The classic application of HAND does not conserve mass because it is typically applied using water elevations obtained from calibrated rating curves at the channel that does not adequately capture the geometry of the floodplain and its effect on energy losses. This can cause unrealistic flood extents. To resolve this limitation, we used the valley-averaged channel and floodplain cross-section to calculate the flow cross-sectional area in the 1D river channel equation. The resulting flow depth (for a given channel slope and roughness) produces an average flow cross-section that when multiplied by the reach length conserves mass (i.e. conservation of mass would be perfect if the floodplain cross-section was constant along the reach length). This causes HAND flood extents to be consistent with mass balance constraints. Discrepancies in the mass balance are introduced when the valley geometry deviates from its average shape at different points along the reach. Thus, the method conserves mass on average at a reach scale, but there is some scatter for individual cross-sections. While it contains limitations, this methodology is a major improvement in the application of the HAND method in rivers.

For coastal ocean and lake floods, HAND simulated flood propagation inland from the coastline, where we added a constant energy dissipation term to reduce storm surge and tides linearly with inland propagation. The energy dissipation term was calibrated based on a wide range of external flood inundation sources of coastal flooding including FEMA and local high-resolution studies.

Thus, the inundation model can be characterized as non-hydrodynamic, mass-conserving at a river reach scale (for fluvial flood), and non-mass-conserving energy-based (for coastal flood). These approximations may lead to differences in inundation at a local scale compared to high-resolution dynamical models, but on average, results compare well to external validation data (see Supplementary Information).

Flood defenses for rivers and coasts were included in the algorithm based on datasets covering the USA, (US National Levee Database, https://nld.sec.usace.army.mil/), the Netherlands (https://www.opendelve.eu/), portions of the UK (UK Environment AIMS 2020 Spatial Flood Defences for England), France (Bureau de Recherches Géologiques et Minières, Flood Zoning – 2020 Report) and portions of Germany (local government databases Baden-Württemberg, Bremen, Lower Saxony). Flood defense information outside these regions was not available at the time of model development but will be included in future model updates when available. This assumption implies potentially overpredicted flood inundation in regions where flood defenses exist on the ground but not in the model. Additionally, we do not consider future climate adaptation measures or improvements to existing flood defenses in the model. The defense logic was based on a return period protection level under baseline conditions (i.e. overtopping occurs if water levels exceed the baseline protection level).

Pluvial model

The pluvial flood model was forced by global extreme precipitation projections derived from CMIP6 GCM results37,38. These extreme precipitation projections were used to simulate the rain-on-grid events.

Specifically, we used the following models (number of ensemble members indicated in parentheses): ACCESS-ESM1-5 (3), CESM2 (3), CESM-LENS (40), CESM Med. Ens. (15), MIROC6 (50), MPI-ESM1-2 (10), and UKESM1-0-LL (6). These were bias-corrected, downscaled to 0.25 degrees, and summarized using a generalized extreme value (GEV) model46 and associated uncertainty. The pluvial model combined an efficient fill-spill-merge (FSM) algorithm that infills low areas47 with a modified rational method model used in combination with Manning’s equation to approximate overland flow depth48. In locations within the USA and Europe where DEM data are available at 10-m resolution, the pluvial model was run at the native higher resolution and then coarse-grained to 90 m. Elsewhere, the model was run at 90 m resolution, which may underestimate pluvial flooding in dense urban areas in the FSM method. We excluded the effects of storm drains or other stormwater control features in the pluvial model. In these cases, the rational method provides a reasonable flooding result for the overland flooding over the pixel on average.

Combined flood model

The combined flood model result was based on the maximum coastal, fluvial, and pluvial flood depth at each pixel for a given return period, year, and climate scenario. The assumption of maximum depth is based on the probability of occurrence of a flood hazard. Because we are evaluating flood risk from coastal storm surge events, small-scale high-intensity pluvial flood events, and large-scale medium-intensity fluvial events, these processes are on average independent. Coincident event types are rare, but certainly possible, leading to compound flood hazards. In general, these joint hazards can increase flood levels in areas where hazards overlap. Thus, these results may underestimate flood depth in some selected areas. While modeling compound events is possible computationally at regional scales11, it is currently not feasible at global scales, and is a topic of future development.

Flood model validation

Extensive validation of the flood results to external datasets is provided in the Supplementary Information.

Flood exposure

Based on the flood model depth results, we created a flood mask F for regions without permanent water,

where D is the flood depth and α is a flood threshold considered meaningful for impacts. The lower threshold limit is needed here where pluvial flooding is modeled because nearly all grid cells have some positive depth and are resolved to 1 cm vertical depth resolution. While there is no universal agreement on this threshold, 10 cm is used for this study and is generally assumed to represent a transition from nuisance flooding to flooding potentially causing damage49. W is a mask of typically dry areas based on the seasonal duration of surface water Sw from the Global Surface Water Explorer dataset at 30 m resolution (https://global-surface-water.appspot.com/)50,

and β is a threshold here selected as 6 months. F is thus a measure of areas within the modeled floodplain of meaningful depth that are not within permanent water bodies (ocean, lakes, river channels, etc.).

To aggregate the results in a meaningful way, we used the geoBoundaries dataset (https://www.geoboundaries.org/) selecting the average administrative level (ADM0,1,2) for each country, which is closest to 4.2×103 km2 on average on a logarithmic scale. For example, in the USA, this results in the ADM2 level (counties). The inundated area I over a given area is then given by the sum of the flooded area

where dA is the grid cell area.

Population exposure

Population used gridded baseline density Pd at 30 arc-sec resolution (approx. 1 km) (GPWv4, https://sedac.ciesin.columbia.edu/data/collection/gpw-v4) for 2020. Future population estimates PT (year, country) were based on the United Nations World Population Prospects24 medium variant projection (https://population.un.org/wpp/). For each country and year, we adjusted the summed gridded population to match the UN projection with a correction factor

The corrected gridded population density P is then given by

For statistical analysis, countries are grouped by the United Nations Geographic M49 regions (https://unstats.un.org/unsd/), Africa, Asia, Australia & Oceania, Europe, Latin America & the Caribbean, and North America.

The population exposed (PE) to flooding in a given area is then the sum of P within floodplain F,

To distinguish urban vs nonurban areas, we applied a mask U using a threshold \(\gamma\) of 1500 ppl/km2 per the World Bank51

Thus, the population exposed to flooding within urban and nonurban areas (PEu, PEnu) is

If we assume a future condition (f ) is the result of a base condition (b) plus a change,

Substituting these definitions into the expression for some future condition PEf,

where the terms from left to right are population exposed from baseline, population change, climate change, and a joint term from coincident climate and population change. This approach to modify the density assumes that the relative population spatial distribution will remain constant within a country and that society will not spatially adapt to climate change. Population exposure data is generally consistent with results from other studies (see Results).

Gross domestic product

The regional gross domestic product was taken from Kummu et al. (2020)52 gridded GDP per capita (Gd) for the year 2015 at 5 arc-min resolution (approx. 9 km). Total country GDP estimates GT(country) for 2019 are based on the Penn World Table pwt100 dataset, expenditure-side real GDP at chained PPPs (in 2017 $M USD)53. For each country, we adjust the summed gridded GDP to match the country GDP with a correction factor

The corrected GDP G for a given region is then given by

Vulnerability

While there are many definitions of vulnerability, here we define the relative vulnerability of a region as an area having high PE and low G; an area that has a relatively low ability to pay for adaptation or climate infrastructure and high exposure of the population to flooding is particularly vulnerable. Additionally, we conducted an analysis of the sensitivity of population exposure to varying flood event magnitudes following Devitt et al.27, and we fit the population exposure to an exponential form

where r is the return period of interest, the highest modeled return period (500-year) is used to normalize the function, and b is the extreme event sensitivity found using the best fit of the 10-, 20-, 100-, and 500-year events for each region. Values of b less than one are defined as populations exposed to low-frequency flood events, while values greater than one are sensitive to high return period events.