Current and future rainfall-driven flood risk from hurricanes in Puerto Rico under 1.5 and 2 ∘C climate change

Abstract

Flooding associated with Hurricane Maria in 2017 had devastating consequences for lives and livelihoods in Puerto Rico. Yet, an understanding of current and future flood risk on small islands like Puerto Rico is limited. Thus, efforts to build resilience to flooding associated with hurricanes remain constrained. Here, we take an event set of hurricane rainfall estimates from a synthetic hurricane rainfall simulator as the input to an event-based rainfall-driven flood inundation model using the hydrodynamic code LISFLOOD-FP. Validation of our model against high-water-mark data for Hurricane Maria demonstrates the suitability of this model for estimating flood hazard in Puerto Rico. We produce event-based flood hazard and population exposure estimates for the present day and the future under the 1.5 and 2 ${}^{\circ }\mathrm{C}$ Paris Agreement goals. Population exposure to flooding from hurricane rainfall in Puerto Rico for the present-day climate is approximately 8 %–10 % of the current population for a 5-year return period, with an increase in population exposure to flooding by 2 %–15 % and 1 %–20 % under 1.5 and 2 ${}^{\circ }\mathrm{C}$ futures (5-year return period). This research demonstrates the significance of the 1.5 ${}^{\circ }\mathrm{C}$ Paris Agreement goal for Small Island Developing States, providing the first event-based estimates of flooding from hurricane rainfall under climate change for a small island.

1 Introduction

Climate change is amplifying the probability of high-intensity tropical cyclone events globally (Patricola and Wehner, 2018; Kossin et al., 2020; Mei and Xie, 2016; Knutson et al., 2020), compounding the rising social and economic costs associated with disasters due to increasing population and asset exposure (Jiménez Cisneros et al., 2014). The adoption of the Paris Agreement in 2015 aimed to limit global warming to well below 2 ${}^{\circ }\mathrm{C}$ above pre-industrial levels and if possible to 1.5 ${}^{\circ }\mathrm{C}$ (United Nations Framework Convention on Climate Change, 2015). Following this, numerous studies have investigated how these global temperature changes could impact societies, ecosystems, and places (IPCC, 2018; Mitchell et al., 2016). Under the Paris Agreement goal of 2 ${}^{\circ }\mathrm{C}$, there will likely be a higher proportion of tropical cyclones that become the most intense storms (i.e. Categories 4 and 5 hurricanes), with an increase in precipitation intensity (Knutson et al., 2020). Whilst flooding accounts for the largest proportion of loss of life and economic damages from tropical cyclones (Rappaport, 2014; Czajkowski et al., 2017), there is a lack of literature exploring how flooding might be affected by changes in tropical cyclone characteristics under climate change. This is particularly pertinent for Small Island Developing States where the difference between the 1.5 and 2 ${}^{\circ }\mathrm{C}$ temperature goals may be critically important (Hoegh-Guldberg et al., 2018).

Small Island Developing States (SIDS) are a group of small island nations and territories with an acute risk of disasters and the impacts of climate change and were an instrumental force in the implementation of the 1.5 ${}^{\circ }\mathrm{C}$ goal in the Paris Agreement (Ourbak and Magnan, 2018). Considering risk as a function of hazard, exposure, and vulnerability (United Nations Office for Disaster Risk Reduction, 2024), high hazard frequency, high exposure in relation to size, and underlying vulnerabilities drive the risk of hydrometeorological disasters and climate change in SIDS (Nurse et al., 2014; Mycoo et al., 2022). Climate change is likely to exacerbate current flood risk in SIDS (Joyette et al., 2014; Thomas et al., 2017) based on projected changes in tropical cyclone precipitation (Vosper et al., 2020), increased coastal storm surge heights (Knutson et al., 2020; Monioudi et al., 2018), and sea level rise (Storlazzi et al., 2018; Nicholls et al., 2018; Rasmussen et al., 2018). Yet, very little island-scale quantitative assessment of flood risk has been conducted in SIDS. This is largely due to the inadequacy of existing methods as well as insufficient data resolution and quality suitable for the scale of small island modelling (typically $<$ 10 000 ${\mathrm{km}}^{\mathrm{2}}$) (Thomas et al., 2019).

Recent work by Vosper et al. (2020) demonstrates that total rainfall associated with tropical cyclones (also known as hurricanes) in the Caribbean will increase under both the 1.5 and the 2 ${}^{\circ }\mathrm{C}$ Paris Agreement goals in comparison to the present-day climate. They also estimate that a 100-year return period event similar to Hurricane Maria in Puerto Rico would be twice as likely to occur under the 2 ${}^{\circ }\mathrm{C}$ scenario than under the 1.5 ${}^{\circ }\mathrm{C}$ scenario (Vosper et al., 2020). Puerto Rico is an unincorporated territory of the US, located in the Greater Antilles islands of the Caribbean (see Fig. 1). The urgent need to understand both current and future flood risk was recently reinforced following Hurricane Maria in 2017, which made landfall as a high-end Category 4 hurricane, causing catastrophic wind and flood damage (Pasch et al., 2018). Hurricane Maria was the strongest hurricane to hit Puerto Rico since Hurricane San Felipe II in 1928, resulting in at least 2975 deaths (Audi et al., 2018). The estimated economic loss of USD 90 billion made it the third costliest disaster in US history (Pasch et al., 2018). Despite the underlying structural failures and inadequate emergency response that also contributed to the scale of the disaster in Puerto Rico (Towe et al., 2020; Rivera, 2020; Caban, 2019; Willison et al., 2019), the volume and intensity of the rainfall associated with Hurricane Maria were unprecedented and exacerbated the scale of the impact on communities on the island (Keellings and Hernández Ayala, 2019; Ramos-Scharrón and Arima, 2019). Historically, hurricane rainfall has been the key cause of flooding in Puerto Rico (Hernández Ayala et al., 2017; Smith et al., 2005). Consequently, it is pertinent that estimates of current and future rainfall-driven flood risk associated with these hurricane rainfall events are developed to assist disaster risk management in Puerto Rico. Yet, there are currently no complete estimates of flooding associated with Hurricane Maria or indeed for any other events in Puerto Rico. Dated Federal Emergency Management Agency (FEMA) flood zone maps do exist for larger river systems in Puerto Rico, but these do not include pluvial flooding, which is a key focus of this paper. They are therefore likely to provide a considerable underestimate of risk (Wing et al., 2017).

Figure 1

Map showing the island of Puerto Rico within the Caribbean region.

[Figure omitted. See PDF]

Tropical cyclones can generate pluvial, fluvial, and coastal floods, all of which interact. Of these, pluvial flooding is a comparatively lesser modelled phenomenon (Blanc et al., 2012; Rözer et al., 2019; Tanaka et al., 2020). Pluvial flooding is defined here as “flooding resulting from rainfall-generated overland flow and ponding before the runoff enters any watercourse or drainage system, or cannot enter it because the network is full to capacity” (Falconer et al., 2009, p. 199). There has been a historical split between the modelling and assessment of pluvial and fluvial – or river – flooding. However, in reality both of these inland flood types are in a continuum, and both are driven by rainfall. Thus, the distinction between the two is unhelpful in many cases. This is particularly true for small islands where much of the inland flooding is primarily driven by heavy rainfall (Jetten, 2016; Burgess et al., 2015). Pluvial flooding is also a contested term, with some defining it as including small river channels (Wing et al., 2018) and others defining it as completely independent of rivers (Rosenzweig et al., 2018; Hankin et al., 2008). The rain-on-grid approach documented here therefore overcomes this pluvial–fluvial distinction by explicitly modelling both flood types and their interactions. Here we define the flooding modelled in this approach as “rainfall-driven flooding”.

Rainfall-driven flood events can often occur with a high frequency but low magnitude. This can lead to a significant cumulative impact on a community's resilience over time, which can undermine efforts to reach the UN's Sustainable Development Goals (Moftakhari et al., 2017; Hamdan, 2015; United Nations Office for Disaster Risk Reduction, 2019). However, most studies investigating flooding under climate change focus on changes in the 100-year flood extent because this is often used as a design standard (Hirabayashi et al., 2013; Arnell and Gosling, 2016; Lehner et al., 2006). This means a gap remains with respect to understanding how smaller, more frequent events might vary under climate change, which is crucially important for improving resilience-building and climate change adaptation needed in local communities (Moftakhari et al., 2017). This paper aims to address this gap by investigating how changing hurricane rainfall characteristics influence rainfall-driven flood risk estimates in the Small Island Developing State of Puerto Rico, with an emphasis on understanding changes in lower-magnitude, higher-frequency events ($<$ 30-year return period).

Figure 2

Diagram outlining the modelling steps involved in simulating the synthetic hurricane rainfall event set and its application in the event-based rainfall-driven flood model.

[Figure omitted. See PDF]

Currently, the predominant method for understanding changes in flooding on small islands under climate change uses changes in precipitation as a proxy for changes in flood hazard, leading to uncertainty in flood hazard changes under climate change (Seneviratne et al., 2021; Ranasinghe et al., 2021). Examples of pluvial hydraulic flood modelling on small islands have previously relied on spatially uniform rainfall estimates derived from historical data for a set of design return period events (World Bank, 2015; Pratomo et al., 2016; Lumbroso et al., 2011). This approach takes a set of rainfall intensity estimates for a given duration and return period, often derived from an intensity–duration–frequency (IDF) curve, using historical rainfall data. Rainfall is typically applied uniformly across a model domain to produce design event flood extents (World Bank, 2015). Yet, this approach does not necessarily represent flooding at a particular return period, as a flood is a signature of the rainfall, the topography, and the topology of a catchment (Guerreiro et al., 2017; Skougaard Kaspersen et al., 2017). More recently, studies have highlighted the importance of representing rainfall spatially and temporally for a more realistic representation of flooding (Aldridge et al., 2020; Bernet et al., 2019; Guerreiro et al., 2017; Schaller et al., 2020). One way of incorporating these features is through an “event set approach”, which involves utilizing an event set of synthetic rainfall events (Nuswantoro et al., 2016; Tanaka et al., 2020). Nonetheless, data such as these are still limited or non-existent – particularly for small islands – and thus the aforementioned traditional approach has, until now, been the only way to represent flood hazards for small islands. Climate change is often assessed by applying an uplift factor to account for changes in rainfall associated with climate change projections (Sayers et al., 2020). However, this approach also fails to account for non-stationary effects of climate change on flooding, including changes to the different spatial and temporal characteristics of rainfall that are important for flood generation (Rosenzweig et al., 2018).

This paper details the first example of an event-based assessment of flood hazard on a small island under current and future climate change. We utilize a synthetic hurricane rainfall dataset (Vosper et al., 2020) as the input to an event-based rainfall-driven hydrodynamic flood model of Puerto Rico. We model rainfall-driven flood hazard and population exposure at the island scale in Puerto Rico (9100 ${\mathrm{km}}^{\mathrm{2}}$), at 20 $\mathrm{m}$ resolution under present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change. As part of this work, we also include novel methodological developments, including the representation of rainfall and river channels in the model. The model is validated against flood hazard simulations using two estimates of Hurricane Maria-observed rainfall (Integrated Multi-satellitE Retrievals for GPM (IMERG) and NCEP Stage IV) and high-water-mark data collected from the event. To our knowledge, these are the first published estimates of rainfall-driven flooding from Hurricane Maria. This work thus demonstrates a step change in the capacity to estimate flood hazard on a small island, superseding the information available using the traditional approaches. Within this, two key questions are investigated:

• What is the current rainfall-driven flood hazard and population exposure associated with hurricanes in Puerto Rico?

• How does population exposure to flooding change from the present day to 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change scenarios?

2 Methods

To address these questions, we first describe the application of the hurricane rainfall event set in Sect. 2.1. We explain how the event-based model was set up (Sect. 2.2), including the novel methodological applications of spatially varying rainfall in the hydrodynamic model (Sect. 2.2.1) and the parameterization of river channel bathymetry using the input rainfall event set climatology (Sect. 2.2.2). In Sect. 2.3, we describe the combination of population estimates with the flood hazard data to derive population exposure estimates under present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change scenarios. The method for validating the model is described in Sect. 2.4.

2.1 Hurricane rainfall data

The synthetic hurricane rainfall event set was developed to estimate hurricane rainfall in the Caribbean under present-day (2005–2016) as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ equilibrated climate change, using a physics-based tropical cyclone rainfall model (Vosper et al., 2020) – see Fig. 2. The model produces spatial (10 $\mathrm{km}$ resolution) and temporal (2-hourly) rainfall estimates along a synthetic hurricane track, considering four key rainfall-generating mechanisms: wind shear, topography, vortex stretching, and surface frictional convergence. Inputs to the tropical cyclone rainfall model were atmospheric temperature, specific humidity, sea surface temperature, and wind vectors, which are typically taken from global climate models (GCMs) or reanalysis products. This model has been validated against gauge-based and radar observations in several studies in the US – including in Puerto Rico – showing good agreement (Feldmann et al., 2019; Lu et al., 2018; Zhu et al., 2013).

Table 1

Table outlining the resolution of the global climate models used to drive the synthetic hurricane rainfall model from the HAPPI climate ensemble.

HAPPI climate model	Horizontal resolution	Number of simulated years of climate model data			Reference
		Present day	1.5 ${}^{\circ }\mathrm{C}$	2 ${}^{\circ }\mathrm{C}$
CanAM4	2.81${}^{\circ }$ $\times$ 2.81${}^{\circ }$	332	346	332	Wehner et al. (2014)
CAM5-1-2-025degree	0.31${}^{\circ }$ $\times$ 0.23${}^{\circ }$	409	365	396	Von Salzen et al. (2013)
ECHAM6-3-LR	1.88${}^{\circ }$ $\times$ 1.88${}^{\circ }$	427	378	383	Stevens et al. (2013)
NorESM1-HAPPI	1.25${}^{\circ }$ $\times$ 0.94${}^{\circ }$	423	382	351	Bentsen et al. (2013)Iversen et al. (2013)Kirkevåg et al. (2013)

To provide driving climate model data to the synthetic hurricane rainfall events under current as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change, four climate models from the Half a degree Additional warming, Prognosis and Projected Impacts (HAPPI) ensemble were utilized (CanAM4, CAM5-1-2-025degree, NorESM1-HAPPI, ECHAM6-3-LR; Mitchell et al., 2017) – see Table 1. Representative Concentration Pathway (RCP) 2.6 was used for model boundary conditions at 1.5 ${}^{\circ }\mathrm{C}$, using a weighted combination of RCP2.6 and RCP4.5 at 2 ${}^{\circ }\mathrm{C}$. These were selected based on the availability of variables at the required atmospheric levels with at least daily temporal resolution for input into the hurricane rainfall model, which are described in Fig. 2. HAPPI was developed to document climate change impacts under 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change above pre-industrial levels and has been a key source of climate data for such studies, including the IPCC Special Report on 1.5 ${}^{\circ }\mathrm{C}$ (IPCC, 2018).

The hurricane rainfall event set consists of 59 000 events, with each climate model scenario equivalent to between 332–427 simulated years of data depending on the climate model (Vosper et al., 2020). A total of 59 000 events were generated, corresponding to approximately 5000 events per climate model and climate scenario. For each climate model, the number of simulated years was calculated as the sum of the number of simulated events per year divided by the simulated annual frequency of events in the climate model data (see Table 1). The simulated time period for the present day is 2005–2016, representing a global average temperature of around 0.9 ${}^{\circ }\mathrm{C}$ higher than a pre-industrial climate. The 1.5 and 2 ${}^{\circ }\mathrm{C}$ time periods are for 2106–2115. This future time period was selected in the HAPPI climate ensemble as the future time slice, chosen to represent a 1.5 and 2 ${}^{\circ }\mathrm{C}$ world at around 2100 (which was the generally accepted time period for these temperature scenarios in the IPCC Special Report on 1.5 ${}^{\circ }\mathrm{C}$; IPCC, 2018), whilst also providing 100 years of simulated GCM data following the present-day time slice (2006–2015) (Mitchell et al., 2017). Each synthetic hurricane rainfall event was simulated at a 2 $\mathrm{h}$ time step and 10 $\mathrm{km}$ spatial resolution before being employed as the input to the event-based rainfall-driven flood model.

2.2 Event-based rainfall-driven flood model

LISFLOOD-FP is the hydraulic engine used to simulate channel and floodplain flow in two dimensions in our rainfall-driven hydrodynamic model (Bates et al., 2010; LISFLOOD-FP Developers, 2020). Rainfall is the key input to the model, and water flow is routed in one of two ways. Firstly, very shallow ($<$ 1 $\mathrm{cm}$) overland flows are routed using a constant-velocity rain-on-grid routing scheme (Sampson et al., 2013). Rain falls directly onto the cells and is routed through the model using a slope-dependent fixed-velocity algorithm. Secondly, flow above 1 $\mathrm{cm}$ deep (i.e. the majority) is routed hydraulically using the inertial form of the shallow-water equations (Bates et al., 2010), with river and drainage channels represented using a subgrid approach (Neal et al., 2012). Typical channel (0.035) and floodplain (0.040) Manning's coefficient friction values were applied. As Puerto Rico is an island, all downstream boundaries are the ocean. The downstream boundary conditions in the model are set to sea level, and this could be used in future work to simulate sea level rise and storm surge.

As digital elevation data are the most important input to a hydrodynamic model (Hawker et al., 2018), lidar data were used as the digital elevation model (DEM). Lidar coverage for Puerto Rico is almost complete ($>$ 99 %) (United States Geological Survey, 2017) and was resampled from its native 1 $\mathrm{m}$ resolution to 20 $\mathrm{m}$, reprojected to WGS84 and hydrologically conditioned using the priority-flood method (Zhou et al., 2016). The $\sim$ 55 ${\mathrm{km}}^{\mathrm{2}}$ of Puerto Rico not covered by lidar was patched with the globally available MERIT DEM (Yamazaki et al., 2017). This area is mountainous and sparsely populated, meaning the use of MERIT here does not affect the exposure results.

Whilst high-resolution DEMs are important for simulating floods, halving the model grid resolution leads to an increase in simulation time by an order of magnitude (Savage et al., 2016). For example, run on a 2 $\times$ 2.6 $\mathrm{GHz}$ 8-core Intel E5-2670, one example model in this event set for the 9100 ${\mathrm{km}}^{\mathrm{2}}$ domain, covering the entire island of Puerto Rico, takes 3 $\min$ to run at 90 $\mathrm{m}$, 77 $\min$ at 20 $\mathrm{m}$, approximately 770 $\min$ (12.8 $\mathrm{h}$) at 10 $\mathrm{m}$, and 7700 $\min$ (5.3 $\mathrm{d}$) at 1 $\mathrm{m}$ resolution. As a result, and given the fact that we had thousands of events to simulate, the event set was run at 20 $\mathrm{m}$. This resolution balances the need for high-resolution flood hazard outputs with the computational costs associated with employing a high-resolution event-based model at the island scale and also reflects state-of-the-art model resolutions used in other locations, such as the UK (Bates et al., 2023). Our study is the first known study to employ an event set approach at such a high hydrodynamic model resolution over such a large domain.

Infiltration was not included in this model approach for several reasons. As hurricanes take place during the hurricane season (North Atlantic: June–November), soils in Puerto Rico are often saturated, meaning infiltration is low (Smith et al., 2005). Many pluvial modelling studies do not include infiltration, as the appropriate parameter values are highly uncertain and vary widely across space and time (Bernet et al., 2018; Guerreiro et al., 2017; Hall, 2015). Although antecedent conditions are expected to vary, the infiltration is likely to be of lower importance relative to other factors, since infiltration will be minimal under extreme rainfall events – such as those associated with hurricanes (Wehner and Sampson, 2021).

To improve the representation of islands and hurricane rainfall in the model, two novel model developments were incorporated into the model set-up.

2.2.1 Spatially varying rainfall

Spatiotemporal representation of rainfall is important for accurate simulation of pluvial flood events (Blanc et al., 2012). Previous pluvial models using LISFLOOD-FP covered only small domains and relied on time-varying but spatially constant rainfall input (Sampson et al., 2013, 2015; Wing et al., 2019). This study demonstrates the first use of spatially and time-varying rainfall in a LISFLOOD-FP rainfall-driven hydrodynamic model, using a new routine to read spatiotemporal rainfall in NetCDF format. For each hurricane, a grid of rainfall at $\sim$ 10 $\mathrm{km}$ resolution across the island was input to the model domain at each time step (2-hourly), although the hydrodynamic model calculations are simulated with much shorter time steps (order of seconds). To model all 59 000 hurricane rainfall events would be computationally intractable and was not necessary considering many of the hurricane rainfall events produced no or very little rainfall. Thus, to select events to simulate in the model, all hurricane rainfall events above a threshold of 3.75 $\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}$ peak rainfall intensity were simulated – a total of 4909 events (8.3 % of the total). Within this, 1464 events were in the present day, 1801 events were at 1.5 ${}^{\circ }\mathrm{C}$, and 1644 events were at 2 ${}^{\circ }\mathrm{C}$. This threshold was selected as the minimum number of events necessary to calculate a robust estimate of the 2-year return period flood hazard, which is used as the lowest modelled return period event in the event set. Events below this threshold were not considered significant enough in terms of rainfall to run. An additional 8 $\mathrm{h}$ of simulation time was added to the end of each simulation based on our inspection of the time it took for the rainfall to move through the model and reach either the ocean or the lowest points of the DEM. These decisions were based on trial and error and inspection of the rainfall and resulting flood hazard events.

2.2.2 River channels

Including river channels in flood models is necessary to produce accurate estimates of flood hazard (Hall, 2015; Neal et al., 2021), but most pluvial flood models do not explicitly include river channels or drainage networks (Blanc et al., 2012). Here, a subgrid approach was used to represent river channels and drainage networks in the rainfall-driven modelling framework (Neal et al., 2012). Rivers and drainage channels were represented using the US National Hydrography Dataset v2.1 (Simley and Carswell Jr, 2010). River widths in Puerto Rico are inadequately represented in global hydrographic datasets such as MERIT Hydro (Yamazaki et al., 2019), as most channels are smaller than the resolution of the DEM data used to create such products (e.g. MERIT at 90 $\mathrm{m}$ in the case of MERIT Hydro). As a result, width was estimated using a power law regression with upstream-accumulated area (Leopold and Maddock, 1953). Widths used here were sampled using satellite imagery along the 13 main rivers across the island. The upstream-accumulated area was calculated using the lidar DEM at 20 $\mathrm{m}$ resolution by first generating a flow direction map and then using the RichDEM algorithm outlined in Barnes (2017).

River depth estimates are also unavailable for Puerto Rico, as is typical for most locations globally (Sampson et al., 2015). To parametrize the river channel depths, the present-day synthetic hurricane rainfall events for each climate model (total: 1464) were first simulated through a model with arbitrarily deep river channels ($-$10 $\mathrm{m}$) to get estimates of the channel water depth for each event. Using these, the water depth at a given return period was calculated empirically. Information on flood defences was also not available, so in this study we parametrize bankfull river depth by calculating the bed elevation to ensure that each channel conveyed the present-day 1-in-2-year discharge (Pickup and Warner, 1976; Williams, 1978; Wolman and Miller, 1960) generated by the present-day hurricane ensemble and subtracted from the bank height derived from the DEM to get a calibrated estimate of the channel depth value. Inevitably this means that in locations where rivers do have defences, the model is likely to overpredict flood hazard. If defence standard information was to become available, it would be a simple matter to retrospectively apply these to the output flood hazard layers.

2.3 Population exposure estimates

Population exposure was calculated for each flood event as the total number of people exposed to flood depths above 10 $\mathrm{cm}$. The WorldPop 90 $\mathrm{m}$ top-down constrained population dataset (2020) was used to estimate the number of people per 90 $\mathrm{m}$ grid cell (Tatem, 2017; Bondarenko et al., 2020). WorldPop was chosen because total population estimates are adjusted to 2020 UN population estimates, meaning out-migration trends following Hurricane Maria in 2017 are accounted for. The WorldPop data were downscaled from 90 $\mathrm{m}$ to 20 $\mathrm{m}$ to match the flood hazard data, using nearest-neighbour resampling and assignment to 20 $\mathrm{m}$ cells based on a proportional cell method, following Lloyd et al. (2017). WorldPop has been extensively validated and compared to other datasets (Reed et al., 2018; Leyk et al., 2019; Tuholske et al., 2021), including for flood exposure applications (Mazzoleni et al., 2020; Smith et al., 2019). Smith et al. (2019) found that WorldPop produces larger exposure estimates in comparison to the High Resolution Settlement Layer (HRSL) (Tiecke et al., 2017), likely due to a combination of coarser resolution and assignment of populations to buildings. Recently, Tuholske et al. (2021) identified the importance of conducting a sensitivity assessment of gridded population products to capture the inherent uncertainties in the use of gridded population estimates. However, HRSL, high-resolution population density maps (HRPDMs) (Bonafilia, 2019), and WorldPop are likely to give different estimates in our case, not least due to the different dates of the datasets before and after Hurricane Maria, where approximately 8 % (230 000) of the population is estimated to have emigrated following the event (Audi et al., 2018). Total population estimates for the main island using HRPDM and HRSL are 4.87 million and 3.66 million, which is considerably higher than the UN-adjusted WorldPop estimate of 2.70 million, resulting in higher population exposure values. The future population was not considered due to a lack of available high-resolution datasets ($<$ 100 $\mathrm{m}$ grid size) estimating changes in the future population. For consistency, population exposure exceedance was calculated for each event using the same method as the hurricane rainfall as 1 divided by the annual exceedance probability (Emanuel and Jagger, 2010; Feldmann et al., 2019; Vosper et al., 2020).

2.4 Model validation

To determine the skill of our flood hazard estimation, we assessed the model performance using high-water-mark (HWM) data collected by USGS following Hurricane Maria (available at https://stn.wim.usgs.gov/FEV/{#}MariaSeptember2017, last access: 24 October 2023). For more information about the suitability assessment of the HWM data for validation, see Sect. S1 and Table S2 in the Supplement. See Fig. S1 in the Supplement for the HWM locations used in this study. Ideally it would be better to validate the event set with a lower-magnitude flood considering the focus of this work is primarily on low-magnitude, high-frequency events. However, there is no known validation data for small hurricane-rainfall-driven flood events in Puerto Rico. As a result, Hurricane Maria was chosen as the event to validate against, despite its high magnitude.

Firstly, to produce flood hazard estimates of Hurricane Maria for validating the model and event set, we ran the hydrodynamic model using two observational rainfall products (IMERG and NCEP Stage IV) that provide space- and time-varying estimates of Hurricane Maria rainfall through the flood inundation model. We use an identical hydrodynamic model set-up to the event set, only changing the input rainfall data. IMERG (NASA Global Precipitation Measurement Mission, 2023) was run at $\sim$ 10 $\mathrm{km}$ spatial resolution and in 30 $\min$ intervals, whilst NCEP Stage IV (NCEP/EMC 4KM Gridded Data (GRIB) Stage IV Data, 2023; Du, 2011) was run at $\sim$ 4 $\mathrm{km}$ spatial resolution, with an hourly temporal resolution. NCEP Stage IV was used instead of the higher-resolution Multi-Radar Multi-Sensor (MRMS) rainfall product as the landfall year of Hurricane Maria (2017) falls outside of the MRMS archive period (2020–onwards) (University of Oklahoma Cooperative Institute for Mesoscale Meteorological Studies et al., 2023).

We compare the flood hazard produced using IMERG and NCEP Stage IV to understand the uncertainty in flood hazard estimates using the different observation inputs.

IMERG has been widely compared to gauge-based rainfall data over many locations globally, demonstrating good performance in the estimation of total rainfall (Freitas et al., 2020; Pradhan et al., 2022), as well as good representation of temporal (Yu et al., 2021) and spatial event structures (Omranian et al., 2018; Rios Gaona et al., 2018; Pradhan et al., 2022). For example, Rios Gaona et al. (2017) show IMERG has a low relative bias over the Netherlands ($-$1.51 %), and Tan et al. (2017) report a correlation coefficient of 0.78 against radar and gauge-based observations in the US. IMERG has also been shown to perform well at capturing rainfall from tropical cyclones (Rios Gaona et al., 2018; Yu et al., 2021). For example, Omranian et al. (2018) found IMERG correctly predicted 62 % of rainfall from Hurricane Harvey. Nonetheless, some studies have identified a tendency for IMERG data to underpredict rainfall intensity during extreme rainfall events (Freitas et al., 2020; Mazza and Chen, 2023; Tian et al., 2018; Yu et al., 2021). For example, Yu et al. (2021) found that extreme precipitation rates from IMERG were 7.53 % lower than those of gauge data for Typhoon Lekima in 2019.

NCEP Stage IV is a ground-based gauge and radar observation product that is often used in multi-product comparison studies as the baseline observed dataset (Nelson et al., 2016). These studies have demonstrated that NCEP Stage IV produces good representation of overall rainfall rates across the US (Nelson et al., 2016; Prat and Nelson, 2015), as well as the spatial and temporal structures of rainfall (Habib et al., 2009), including for tropical cyclones (Gao et al., 2020; Villarini et al., 2011). Prat and Nelson (2015) compare annual rain rates for the conterminous US using NCEP Stage IV against gauge data, finding a correlation coefficient of 0.93 ($R^{\mathrm{2}}$). Gao et al. (2020) show that NCEP Stage IV only overestimated rainfall from Hurricane Harvey by 2 %. However, underestimation of extreme rainfall has been shown in some studies to be due to an increase in the number of missed events as the rain rate increases (Habib et al., 2009; Prat and Nelson, 2015). For example, Prat and Nelson (2015) report that NCEP Stage IV has a tendency to underestimate rainfall in comparison to surface observations across the conterminous US ($-$14 % to $+$1 %, depending on the location). This is likely a product of the inherent limitations of radar-based precipitation products (see Nelson et al., 2016).

The model used to produce the synthetic hurricane rainfall event set utilized in this study has previously been compared to NCEP Stage IV data over Puerto Rico, showing very good agreement (Feldmann et al., 2019). This demonstrates the suitability of the use of NCEP Stage IV as an observation dataset for comparison against in this study. Omranian et al. (2018) showed IMERG was able to represent 62 % of rainfall from Hurricane Harvey in comparison to NCEP Stage IV, thus suggesting that IMERG is also likely capable of adequately representing extreme rainfall associated with Hurricane Maria. However, the performance of IMERG and NCEP Stage IV data can be dependent on the number of gauge-based observations available (Tang et al., 2018; Tian et al., 2018). A total of 14 out of 24 USGS gauges were damaged during Hurricane Maria in Puerto Rico (Bessette-Kirton et al., 2020). As a result, this is a key limitation of using observed data products to estimate tropical cyclone rainfall that should be considered when drawing conclusions about the accuracy of flood hazard associated with these rainfall products.

Next, we compared the performance of the event set with the HWM data and the estimates from the observed rainfall products to sense check the model. Hurricane Maria-like events were identified across all model scenarios, first by the maximum total rainfall and then by the spatial characteristics of the hurricane track. Maximum total rainfall is defined as the highest total rainfall accumulation at a point on the island. This metric was used as opposed to mean total rainfall, as studies that have investigated Hurricane Maria rainfall describe the maximum total rainfall as the most significant anomaly in the historical record associated with the event (Ramos-Scharrón and Arima, 2019; Keellings and Hernández Ayala, 2019; Pokhrel et al., 2021). Maximum total rainfall is also the metric used to estimate the return period of Hurricane Maria rainfall: at least a 1-in-115-year rainfall event (Keellings and Hernández Ayala, 2019). Studies use different metrics to derive maximum total rainfall, including interpolation of rain gauge data and observation products such as NCEP Stage IV. This means that the maximum total rainfall for Hurricane Maria varies between studies, ranging between 733–1029 $\mathrm{mm}$ (Pasch et al., 2018; Keellings and Hernández Ayala, 2019; Ramos-Scharrón and Arima, 2019; Pokhrel et al., 2021). There is a limited number of events in our event set with a $>$ 100-year return period magnitude maximum total rainfall (mean: 3.46 samples per climate model scenario) due to the comparatively short simulated time record of our event set (range: 332–427 years). However, Puerto Rico experiences on average one hurricane each year and has a mean annual rainfall of over 4000 $\mathrm{mm}$ in some locations (Hernández Ayala and Matyas, 2016). There are therefore many events in the event set with a total mean rainfall (total accumulated rainfall averaged across the island) in the range of Hurricane Maria (range: 375–380 $\mathrm{mm}$; Pokhrel et al., 2021; Keellings and Hernández Ayala, 2019; Ramos-Scharrón and Arima, 2019). However, these events have widely varying spatial characteristics and associated flood hazard and are therefore not all Maria-like. Thus, it is also important to consider the spatial characteristics of the hurricane rainfall events so that events with similar rainfall and spatial characteristics to Hurricane Maria can be identified. Similarity to Hurricane Maria based on track location was assessed based on four criteria: (i) direct landfall on the main island; (ii) south-western trajectory; (iii) whether it makes landfall on the eastern portion of the main island; and (iv) a similar track trajectory across the island, whereby the event track and Hurricane Maria track intersect at at least one point on the island.

Figure 3

Map showing the differences between flood hazard estimates of Hurricane Maria produced using IMERG and NCEP Stage IV precipitation data from the island to local scale.

[Figure omitted. See PDF]

Figure 4

Graph showing the relationship between the root mean square error (RMSE) and maximum total rainfall for all simulated events under all climate scenarios (4909 events in total). All simulated events are indicated in blue. Events identified with Hurricane Maria maximum rainfall totals and spatial characteristics (20 events) are shown in red. The range of reported Hurricane Maria rainfall is indicated with the red band. The orange square is the NCEP Stage IV model estimate, and the brown triangle is the IMERG model estimate.

[Figure omitted. See PDF]

3 Results

3.1 Hurricane Maria model validation

Figure 3 shows the flood hazard estimates produced by simulating the IMERG and NCEP Stage IV rainfall products spatiotemporally through the flood inundation model from the island to local scale. The RMSE between the modelled flood hazard and the HWM is 1.18 $\mathrm{m}$ for IMERG and 1.22 $\mathrm{m}$ for NCEP Stage IV (see Fig. 4). This is comparable to post-event HWM validation done in other locations (Wing et al., 2021) (see Sect. 4.1 for a discussion of this). There is a significant difference in the flood extents produced using IMERG and NCEP Stage IV, with larger areas flooded using NCEP Stage IV than IMERG. This highlights the uncertainty in so-called “observed” flooding from Hurricane Maria.

In the event set, when the spatial characteristics of the hurricane rainfall events are considered in addition to the maximum total rainfall, events we select as Hurricane Maria-like events have some of the lowest RMSEs between the observed and modelled water surface elevations (range: 1.13–1.33 $\mathrm{m}$), as demonstrated in Fig. 4. The track locations of these events are shown in Fig. S2. The relationship between maximum total rainfall and RMSE for all events is expected, whereby as the intensity of the event increases, the sensitivity to the flood depths decreases as the floodplain fills and thus becomes less responsive to additional increases in rainfall (Wing et al., 2021). However, there are events in the event set with both much higher and much lower rainfalls than Hurricane Maria that have RMSEs that are both similar to and very different from those of the Maria-like events. This demonstrates the importance of the spatial characteristics of the events beyond just the rainfall.

When comparing the flood estimates using IMERG and NCEP Stage IV against the high-water-mark data, the Maria-like events have similar RMSE scores (Fig. 4). However, both observational rainfall products have different maximum total rainfalls than those found in the literature. In particular, the IMERG maximum total rainfall is considerably lower. This is likely because satellite products such as IMERG often underestimate orographic rainfall such as that exhibited over Puerto Rico (Dinku et al., 2008).

3.2 Design return period flood hazard maps

The probability of inundation was calculated for each pixel in the model domain, calculating how many times each pixel would be inundated above a 10 $\mathrm{cm}$ depth in each climate model temperature scenario. The return period of inundation in each pixel was then determined by calculating how many times we expect a pixel to flood based on the number of years of data simulated (range: 332–427 years, depending on the climate model). Using this, we derived a set of return period flood hazard maps, which provide a spatially explicit representation of a given return period flood event under present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ warming. This supersedes any currently available hurricane-rainfall-driven flood risk information in Puerto Rico, both under current and future climate change. This approach also moves beyond the traditional uplift approach often used in flood risk assessment under climate change, as it provides spatially explicit flood hazard information for a given return period at the island scale and at high resolution.

Figure 5 highlights the scale and detail of flood hazard information using this approach, from the island scale (Fig. 5a) to the local scale (Fig. 5c). For example, Fig. 5c shows flooding at the street level in Levittown, Toa Baja – a town significantly impacted by flooding from Hurricane Maria in 2017 (National Weather Service, 2017).

Based on this example for a 20-year return period flood hazard event using the ECHAM6-3-LR climate model, several schools and hospitals would likely be impacted under present-day and 1.5 ${}^{\circ }\mathrm{C}$ climate change. The estimated flooded area of the 20-year return period flood increases under 1.5 ${}^{\circ }\mathrm{C}$ climate change in comparison to present-day climate change (2006–2015) (Fig. 5c), meaning areas currently not at risk are affected at 1.5 ${}^{\circ }\mathrm{C}$ climate change. Changes at 2 ${}^{\circ }\mathrm{C}$ are similar to 1.5 ${}^{\circ }\mathrm{C}$ but are not shown in Fig. 5 for presentation purposes.

Figure 5

Map showing the 20-year return period flood based on the probability of inundation under present-day and 1.5 ${}^{\circ }\mathrm{C}$ climate change for the ECHAM6-3-LR climate model. (a) Flooding at the island scale. (b) Flooding in the Toa Baja and Cataño districts. (c) Flooding in Levittown, Toa Baja. For presentation purposes, only inundation probabilities under present-day and 1.5 ${}^{\circ }\mathrm{C}$ climate change are shown here.

[Figure omitted. See PDF]

Flooding in the northwest of the island shown in Fig. 5a (latitude–longitude location: 18.3, $-$67.0 to 18.4, $-$66.5) is a feature of the topography and model structure, not data error. This area is dominated by karst hydrology (Hughes and Schulz, 2020). Therefore, these areas of pooled water would likely not feature if karst processes were explicitly represented in the model set-up. The inclusion of karst processes was beyond the scope of this study, and as this area is sparsely populated it is unlikely to impact the presented estimates of population exposure.

3.3

Characterizing changes in population exposure under present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change

This research estimates changes in population exposure to hurricane-rainfall-driven flooding for the island of Puerto Rico under present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change. The climate change scenarios are analysed for each individual climate model, as opposed to the aggregate results, as there are important differences between models that are obscured when using the mean. This is a way of investigating uncertainty explicitly, by understanding the differences between models. Studies such as Daron et al. (2021) have highlighted the importance of assessing individual model performance when climate models give a wide range of projections.

Figure 6

Graph showing population exposure exceedance for present-day, 1.5 ${}^{\circ }\mathrm{C}$, and 2 ${}^{\circ }\mathrm{C}$ climate change, as well as the number of samples in each climate model at a given return period (dotted line). Population exposure above the 30-year return period is faded to represent the uncertainty associated with the limited number of samples at these return periods.

[Figure omitted. See PDF]

Figure 6 shows the return period of a given exceedance of population exposure from hurricane-rainfall-driven flooding in Puerto Rico under present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change. Return periods of population exposure exceedance above the 30-year return period are not considered and are thus faded in Fig. 6. The number of samples for each climate model scenario above the 30-year return period is too small (mean: 12.7 samples) to determine accurate estimates of population exposure above the 30-year return period (see Fig. 6). Thus, changes in population exposure above the 30-year return period in this event set are subject to significant uncertainty resulting from limited samples at these event magnitudes and are therefore not considered further in this analysis. A much longer event set would be required to simulate robust changes in population exposure at higher-magnitude return periods.

Figure 7

Plot showing the percentage of the population exposed to flooding under present-day, 1.5 ${}^{\circ }\mathrm{C}$, and 2 ${}^{\circ }\mathrm{C}$ climate change and the difference between the three scenarios for each HAPPI climate model. The green dot represents present-day population exposure (as a percentage of the total population), with the orange and purple dots representing the population exposure (%) at 1.5 and 2 ${}^{\circ }\mathrm{C}$. The difference between the population exposure between the different scenarios is represented by the line between the dots.

[Figure omitted. See PDF]

Three of the four climate models show agreement in the direction of change between present and future climate change, with increases in population exposure associated with a given return period at 1.5 and 2 ${}^{\circ }\mathrm{C}$ compared to the present day. However, one climate model (CanAM4) shows the opposite trend above the 10-year return period (see Fig. 7). One key reason for this is likely to be the differences in the resolution of the underlying global climate model (GCM) data: the CanAM4 GCM has a coarser resolution (2.81${}^{\circ }$ $\times$ 2.81${}^{\circ }$) than the next most coarse GCM, ECHAM6-3-LR (1.88${}^{\circ }$ $\times$ 1.88${}^{\circ }$) (see Table 1). As a result, the underlying variables driving extreme hurricane rainfall are less likely to be well represented in CanAM4 compared to the other three climate models. It is well understood that higher-resolution GCMs are better able to simulate the underlying conditions important for the development of extreme rainfall and tropical cyclones (Knutson et al., 2020).

Present-day population exposure to flooding from hurricane rainfall in Puerto Rico is approximately 2 %–5 % at the 2-year return period, rising to 8 %–10 % at the 5-year, 9 %–12 % at the 10-year, and 11 %–14 % at the 20-year return periods (see Fig. 6). These are the first published estimates of present-day population exposure to flooding in Puerto Rico. It is difficult to corroborate the population exposure estimates with those for previous events in Puerto Rico due to a lack of data; however these estimates are plausible in light of the universal islandwide flash flood warning given to Puerto Rico during Hurricane Maria (Pasch et al., 2018).

As shown in Fig. 7, the estimated number of people exposed to flooding from hurricane rainfall on average every 2 years would increase by the largest percentage across the different return periods (20 %–140 % at 1.5 ${}^{\circ }\mathrm{C}$; $-$3 %–85 % at 2 ${}^{\circ }\mathrm{C}$). The lower bound here represents the results from the CanAM4 model, which has the lowest GCM resolution (see Table 1). The reason for the widest range at the 2-year return period could be because of the different bed elevation sizes at the historical 2-year return period for each climate model. For a return period population exposure of 5 years as shown in Fig. 8, the percentage increase in population exposure at 1.5 and 2 ${}^{\circ }\mathrm{C}$ ranges from 2 %–15 % and 1 %–20 %, respectively. This is a considerably lower range than the population exposure exceedance at the 2-year return period but also shows more agreement between the climate models.

Figure 8

Bar graph showing the number of people exposed to flooding under present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change for the 5-year population exposure exceedance for each HAPPI climate model.

[Figure omitted. See PDF]

As shown in Fig. 7, there is a notable difference in population exposure exceedance between the present day and 1.5 ${}^{\circ }\mathrm{C}$ in three of the four climate models but a less clear difference between 1.5 and 2 ${}^{\circ }\mathrm{C}$. In two of the four climate models (CAM5-1-2-025degree and ECHAM6-3-LR), the percentage of the population exposed at a given return period is higher at 1.5 ${}^{\circ }\mathrm{C}$ compared to 2 ${}^{\circ }\mathrm{C}$ and in one climate model (NorESM1-HAPPI) higher at 2 ${}^{\circ }\mathrm{C}$ compared to 1.5 ${}^{\circ }\mathrm{C}$. In the CanAM4 climate model, depending on the return period, the percentage of population exposure varies between the three climate scenarios, and no consistent pattern is shown between the three across different return periods.

Figure 8 demonstrates that the range in the absolute population exposure numbers estimated for a given return period between the four climate models is the same as or greater than the percentage uplift in population exposure associated with 1.5 and 2 ${}^{\circ }\mathrm{C}$, highlighting the range of possible absolute population exposure estimates. For the 5-year return period, present-day absolute population exposure ranges from 217 000 (ECHAM6-3-LR) to 264 000 (CAM5-1-2-025degree). This is a 21 % difference, whereas the highest population exposure increase is 22 % between the present day and 2 ${}^{\circ }\mathrm{C}$ for the NorESM1-HAPPI climate model. This underlines the difficulty in estimating current population exposure to flooding. This is not only the case in data-sparse areas such as Puerto Rico, but also in data-rich areas such as the conterminous US (Bates et al., 2021). However, the direction of change between present-day and future climate change (1.5 and 2 ${}^{\circ }\mathrm{C}$) is robust across three of the four climate models, meaning the signal in population exposure to flooding is observable when comparing present-day and future climate change, despite the uncertainty in absolute terms.

4 Discussion

Our estimates of flood hazard and population exposure driven by hurricane rainfall under current and future climate change supersede previous efforts to estimate hurricane-rainfall-driven flood risk in Puerto Rico. Previous estimates rely on local-scale FEMA fluvial assessments or global large-scale assessments that most often neglect small islands through the choice of scale, although the FEMA models will likely be more accurate locally where they exist, depending on the local river channel and flood defence information that was available to the model developers. This study is one of the first known published studies which propagates spatially and temporally explicit hurricane rainfall data to the impact modelling of flood hazard and population exposure estimates and is the first to do so for a small island. Utilizing hydrodynamic flood models to understand changes in flooding under climate change is a critical gap in the literature, despite the widespread use of hydrodynamic models to assess current flood risk. The latest IPCC AR6 Working Group I report has demonstrated that changes in rainfall were still the dominant method used to assess changes in pluvial flooding under climate change (Seneviratne et al., 2021). However, here we find that the changes in population exposure between present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change do not correspond linearly with changes in hurricane rainfall using the HAPPI climate models (Vosper et al., 2020) analysed here, and therefore this rainfall proxy method may not be appropriate when investigating changes in flooding from hurricane rainfall.

4.1 Validating an event-based model

We present the first estimates of rainfall-driven flooding from Hurricane Maria using IMERG and NCEP Stage IV precipitation data. Comparison against HWM data from Hurricane Maria showed that the RMSE of these estimates was reasonable given the typical uncertainties in data of this type (IMERG: 1.18 $\mathrm{m}$, NCEP Stage IV: 1.22 $\mathrm{m}$). There is uncertainty associated with the HWM vertical datum transformation using VDatum ($+$-0.92 $\mathrm{m}$), which is likely to impact the RMSE. However, these RMSEs have a similar magnitude to studies conducted in data-rich regions with similar-quality HWM data, such as the conterminous US ($\sim$ 1 $\mathrm{m}$) (Wing et al., 2021). This shows that the model is capable of realistically simulating flood depths, demonstrating the suitability of the model for estimating flood hazard under current and future climate change. Inevitably, this finding should be considered alongside the inherent limitations when comparing flood estimates to high-water-mark data. For example, RMSEs in this study are higher than in studies such as Neal et al. (2009) (RMSE: 0.28 $\mathrm{m}$). Yet, the HWM data in this study are arguably lower-quality data due to the catastrophic nature of the hurricane, which limited accessibility for post-event assessment due to widescale infrastructure failure (Main et al., 2021). HWMs in this study are concentrated in populated areas and were probably constrained to where it was safe to travel immediately post-event. The performance of the model is likely biased towards these coastal, more populated areas. However, this is also where a considerable portion of the risk is on the island, as this is where the majority of the population resides.

Moreover, there are limitations of the observation precipitation datasets used, which propagate into the flood estimates. Many studies have compared the performance of NCEP Stage IV and IMERG rainfall data (Li et al., 2022; Mazza and Chen, 2023; Omranian et al., 2018; Villarini et al., 2011). Tropical cyclone precipitation in the conterminous US between 2002–2019 was much higher in NCEP Stage IV than in satellite products such as IMERG (Mazza and Chen, 2023). Other studies support this conclusion and find that the explanation for this difference is more likely an underestimation of other products and not an overestimation bias in NCEP Stage IV itself (Villarini et al., 2011). For example, IMERG is likely to underestimate orographic rainfall, which could explain why the flood extent using IMERG is lower than that using NCEP Stage IV (see Fig. 3). This provides an incentive for the event set approach outlined in this study, as it allows for a consideration of a wider range of plausible events to get a greater understanding of uncertainty than just the observations.

Based on the events selected as Hurricane Maria-like, highlighted in Fig. 4, we find that our event set contains events like Hurricane Maria and that these events have amongst the lowest RMSEs in comparison to observed HWMs from Hurricane Maria (range: 1.13–1.33 $\mathrm{m}$). It was expected that given the extreme magnitude of Hurricane Maria ($\sim$ 115-year return period hurricane rainfall event; Keellings and Hernández Ayala, 2019), there would be a limited number of events in our event set with this magnitude due to the comparatively short simulated time record of our event set (range: 332–427 years per climate model scenario). In our event set across all climate model scenarios, we find 20 events that we classify as Hurricane Maria-like based on maximum total rainfall and spatial characteristics. This finding firstly reinforces just how extreme Hurricane Maria was, following the devastating impact on both the population and the infrastructure (Audi et al., 2018; Michaud and Kates, 2017; Main et al., 2021), as well as the literature examining the event in the context of the historical record (Keellings and Hernández Ayala, 2019; Ramos-Scharrón and Arima, 2019). This also indicates that the model has the capacity to replicate events such as Hurricane Maria when both maximum total rainfall and spatial characteristics are considered. Two key conclusions can be taken from this. Firstly, this highlights the importance of variables other than rainfall when estimating rainfall-driven flooding, such as the spatial characteristics of the hurricane, including landfall location and trajectory. Just considering the rainfall was not sufficient to identify Maria-like events. As a result, simulating the spatial and temporal distributions of the rainfall in an event set is a crucial step needed to accurately represent the relationship between hurricane rainfall and flood hazard in Puerto Rico. This finding reinforces previous research, which has identified the importance of hurricane landfall and spatial location for the generation of floods in Puerto Rico (Hernández Ayala et al., 2017; Hernández Ayala and Matyas, 2016; Smith et al., 2005). Secondly, considering that there is uncertainty in so-called observed flooding from Hurricane Maria (see Fig. 3), the event set provides the opportunity to assess many more realizations of events with similar characteristics to Hurricane Maria than those available just using observations. This may allow for a better understanding of uncertainty in rainfall-driven flooding for a given event and thus a greater understanding of risk. Future research investigating changes in flooding from hurricane rainfall should thus take an event-based approach as outlined in this study.

4.2 Current population exposure to flooding from hurricane rainfall

Our results highlight the first published estimates of population exposure to flooding in Puerto Rico under the present-day climate, with approximately 8 %–10 % of the population currently exposed to flooding from hurricane rainfall at the 5-year recurrence interval. This level of population exposure has important implications for resilience to floods. It also underlines the exposure to hydrometeorological hazards already experienced in SIDS, which is a key reason for their high risk to climate change and disasters (Thomas et al., 2020). It is also worth noting that these population exposure estimates are for the present-day (2005–2016) climate at around 0.9 ${}^{\circ }\mathrm{C}$ of global mean warming and therefore do not represent a pre-industrial climate. This means population exposure estimates for the present day identified in this study are likely to be already influenced by climate change, given the significant impact of climate change on recent hurricane rainfall events, such as Hurricane Maria, in Puerto Rico (Keellings and Hernández Ayala, 2019; Patricola and Wehner, 2018).

4.3

Population exposure to flooding from hurricane rainfall under 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change

The results presented in this research estimate that population exposure to flooding from hurricane rainfall will be amplified under 1.5 and 2 ${}^{\circ }\mathrm{C}$ in all but one of the four HAPPI climate models analysed. The Paris Agreement includes the 1.5 ${}^{\circ }\mathrm{C}$ target as the higher ambition goal and is often touted as our best chance to limit the impacts of climate change to within a “safe limit”. However, our analysis contributes to the discourse SIDS have been highlighting for some time now, which is that even a 1.5 ${}^{\circ }\mathrm{C}$ temperature rise above pre-industrial levels leads to a serious threat to adaptive capacity (Ourbak and Magnan, 2018; Mycoo, 2018; Hoegh-Guldberg et al., 2018; Mycoo et al., 2022). Here, we find that even at 1.5 ${}^{\circ }\mathrm{C}$, the increase in population exposure associated with hurricane-rainfall-driven flooding in Puerto Rico is enhanced for events with a return period below 30 years. This may have wide-reaching implications for the resilience of Puerto Rico's population. Moreover, although the 1.5 ${}^{\circ }\mathrm{C}$ goal is technically feasible (IPCC, 2018, 2021), it is not currently the most likely temperature rise based on existing policy pledges. At the time of writing, global temperature increase had already reached $\sim$ 1.1 ${}^{\circ }\mathrm{C}$ above pre-industrial levels (World Meteorological Organization, 2021). Based on our analysis, it is likely that flood hazard and population exposure would increase further under higher-warming scenarios. These changes are likely to vary between GCMs.

Due to the range in both absolute population numbers and the relative changes in population exposure between present-day as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ across the four climate models in this event set, there is uncertainty in both how many people might be exposed to a particular flood event and how much this may change in the future. Moreover, the range in present-day absolute population numbers is often larger than the climate signal, which underlines the difficulty in understanding current population exposure (Bates et al., 2021). This demonstrates the importance of assessing a range of different climate model projections to understand the range in uncertainties, which taking an event set approach enables because it allows for many more realizations of a given event magnitude than is likely to have occurred in the historical record to be considered. Overall, three of the four climate models utilized in this study show that there is a difference in the percentage of the population exposed at a given return period under 1.5 or 2 ${}^{\circ }\mathrm{C}$ climate change in comparison to the present day. It is likely that the difference between 1.5 and 2 ${}^{\circ }\mathrm{C}$ is too small to determine a robust directional change above variability, particularly as only 4 of the $>$ 50 HAPPI ensemble members are utilized in this analysis. Other studies have also shown a spread around the median in precipitation, flood hazard, and population exposure estimates under future scenarios (Bates et al., 2021; Swain et al., 2020; Lopez-Cantu et al., 2020), as well as uncertain differences between 1.5 and 2 ${}^{\circ }\mathrm{C}$ given the influence of underlying uncertainty in the GCM and precipitation data (Uhe et al., 2019).

Other reasons for uncertainty in absolute population exposures likely stems from the choice of population data and the corresponding methodology used to assign population to pixels, as well as the underlying population data used to inform the population totals. This is evidenced by the differences in total population among WorldPop, HRSL, and HRPDM, as discussed in Sect. 2.3. Moreover, flood defences are not included in the model due to a lack of available data, meaning the absolute population exposure numbers – particularly for the lower return periods where flood defences are most likely to provide protection – will probably be an overestimate in some locations. If flood defence information was available, the standards of protection could be applied to the exposure estimates provided in this dataset to estimate population exposure when flood defences are included. On the other hand, as this study does not include estimates of coastal flooding, the population exposure estimates may also be an underestimate. This means that it is important to consider that the exposure estimates outlined in this study are for inland rainfall-driven flooding only.

4.4 Limitations of event set size

Population exposure estimates above the 30-year return period are subject to significant uncertainty due to the limited number of samples (mean of $<$ 12.7 samples across the four climate models) available in the event set with these return periods. As a result, the changes in population exposure between current as well as 1.5 and 2 ${}^{\circ }\mathrm{C}$ above the 30-year return period were not considered in this study. This was an acceptable trade-off based on the current work, as this study was mostly focused on understanding changes in lower-magnitude, higher-frequency events. Flood events with a $>$ 30-year return period are often valley-filling, and therefore the impact of such events is already likely to be very significant for the population, as demonstrated during Hurricane Maria (Pasch et al., 2018). Larger events also often lead to a greater domestic and international response. However, smaller, more frequent events lead to the erosion of resilience in communities over time and do not receive the same level of relief or response (Hamdan, 2015; Bull-Kamanga et al., 2003; Allen et al., 2017; United Nations Office for Disaster Risk Reduction, 2019). Research to date has also mostly focused on changes in the 100-year return period event (Arnell and Gosling, 2016; Lehner et al., 2006; Hirabayashi et al., 2013). Therefore, assessing changes in lower-magnitude, higher-frequency events was a key aim of this study. To detect changes in the 100-year return period population exposure, a much longer event set would be required to detect a significant change between 1.5 and 2 ${}^{\circ }\mathrm{C}$. Although we have shown that 20 events like Hurricane Maria do occur in the event set overall, preferably there would be at least 30–50 events to have confidence in relative changes, as is shown in Fig. 6. This would require at least 1000 years of synthetic data per climate model as a minimum. This should be considered in the future when producing event sets derived from GCMs with the intention to utilize these in flood impact modelling. Inevitably, running a much larger ensemble comes at the expense of computational cost; therefore a trade-off, particularly with inundation model resolution, is likely to be necessary.

5 Conclusions

We present the most detailed estimates of present-day and future (1.5 and 2 ${}^{\circ }\mathrm{C}$) hurricane-rainfall-driven flood hazard and population exposure estimates in Puerto Rico to date. This analysis quantifies present-day population exposure to flooding in Puerto Rico for small- to medium-sized events ($<$ 30-year return period). Population exposure to flooding is likely to increase under both 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change. Estimates here suggest that for the present-day, 8 %–10 % of the total population of Puerto Rico would be exposed to flooding (defined as residing at a location with inundation depth $>$ 10 $\mathrm{cm}$) from hurricane rainfall every 5 years, increasing by 2 %–15 % and 1 %–20 % at 1.5 and 2 ${}^{\circ }\mathrm{C}$, respectively. Increases in the number of people exposed to small- to medium-sized flood events ($<$ 30-year return period) could have a cumulative negative impact on the long-term resilience of the Puerto Rican population without appropriate adaptation. Uncertainty in absolute population exposure estimates, as well as the range in estimated percentage increases in flooding under 1.5 and 2 ${}^{\circ }\mathrm{C}$, should be considered when using these estimates to inform appropriate adaptation.

Through validation of our model in comparison with observed high-water-mark data for Hurricane Maria ($\sim$ 115-year return period rainfall event), we find that our model is able to replicate similar levels of flooding to that which occurred and that there are events like Hurricane Maria in the event set when events with both similar maximum total rainfall and similar spatial track characteristics are considered. This has important implications for future research, as an event-based approach allows for the assessment of many more plausible scenarios than is available in the observed historical record.

Puerto Rico is predicted to experience increased population exposure to flooding associated with hurricane rainfall in the future under 1.5 and 2 ${}^{\circ }\mathrm{C}$ climate change. These findings add to the growing body of research that highlights the critical and disproportionate risk climate change poses to Small Island Developing States, amidst the uncomfortable irony that they are among the lowest contributors to greenhouse gas emissions responsible for anthropogenic climate change (Hoegh-Guldberg et al., 2018; Thomas et al., 2020). This simultaneously highlights the impact of every increment of global temperature increase on Small Island Developing States and thus the importance of high-ambition mitigation efforts, as well as the urgent need for increased climate change adaptation and disaster risk reduction in the region.

Data availability

The HAPPI climate model data described in Mitchell et al. (2017), 10.5194/gmd-10-571-2017, can be found and downloaded under an Attribution-NonCommercial-ShareAlike 2.0 generic license at https://www.happimip.org/happi_data/ (HAPPI, 2024). The lidar data can be found on the USGS Data Access Viewer: https://coast.noaa.gov/dataviewer/#/lidar/search/where:ID=8630 (USGS, 2024). The LISFLOOD-FP 8.0 hydrodynamic model (version 8.0) (software) is available at 10.5281/zenodo.4073011 (LISFLOOD-FP Developers, 2020). The WorldPop population data can be found in Bondarenko et al. (2020), 10.5258/SOTON/WP00684, under a Creative Commons Attribution 4.0 international license. The high-water-mark data can be found on the USGS Flood Event Viewer: https://stn.wim.usgs.gov/FEV/#MariaSeptember2017 (USGS, 2024). IMERG data can be downloaded from the Global Precipitation Measurement database at https://gpm.nasa.gov/data/imerg (NASA, 2024). NCEP Stage IV data can be downloaded from Du (2011). The probability-of-inundation (event set) flood hazard maps from Archer et al. (2023) are available via data.bris at 10.5523/bris.2qtinf5lw52u52snyl5ruwekef under a Creative Commons CC BY-NC 4.0 license.

The supplement related to this article is available online at: https://doi.org/10.5194/nhess-24-375-2024-supplement.

Author contributions

LA conceptualized the idea, conducted the analysis, developed the methodology, validated the results, and wrote the paper; JN and PDB contributed to conceptualization, analysis, methodology, and validation and supervised the research; EV, DC, and JS were involved in the data curation and methodology; DM contributed to conceptualization. All authors were involved in reviewing and editing the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

Paul Bates is supported by the Royal Society Wolfson Research Merit Award. Emily Vosper is supported by the UKRI ERSPC Centre for Doctoral Training.

Financial support

This research has been supported by the UK Research and Innovation GW4$+$ doctoral training partnership (grant nos. NE/S007504/1, NE/S003061/1, and NE/S006079/1).

Review statement

This paper was edited by Yves Bühler and reviewed by two anonymous referees.