Main

Coastlines across the globe are vulnerable to the joint occurrence of high sea levels and intense rainfall1,2,3, which can increase flooding beyond the level predicted by considering either hazard alone and result in compound floods4,5. Coastal compound floods are most often triggered by cyclonic storm events, either tropical cyclones (TCs) or extra-tropical cyclones (ETCs)3, which are both low-pressure systems that can generate hazardous storm surges and rainfall5. The future incidence of coastal rainfall and storm tides may be affected by the combination of sea-level rise (SLR) and changes in storm climatology. Recent projections of storm climatology change suggest an increase in the probability of joint rainfall–surge events along much of the global coastline, mostly driven by an increase in rainfall hazard6,7. Previous studies of US compound flood potential have considered changes in the joint hazard resulting from changes in a subset of climate-induced variables, such as SLR8 and changes in river flow9 or rainfall10.

Along the US Atlantic and Gulf coasts, TCs are one of the largest drivers of coastal flood losses11,12. Although less frequent than ETCs at mid–high latitudes, TCs typically dominate the upper tail distribution (>50 year return period) of both storm surges13,14 and rainfall-induced flooding15,16, and TCs have been responsible for many extreme compound floods1,17. However, few regional studies of compound flood hazards have explicitly accounted for TC events10 due to their sparse occurrence in the historical record and challenges in representing TCs within reanalysis datasets and typical global circulation models (GCMs)7. It is unclear how future changes in TC climatology and SLR will alter the severity and spatial variation of extreme rainfall–surge hazard across the US Atlantic and Gulf coasts, what will be the relative contribution of storm climatology change and SLR to changes in the joint hazard and how changes in TC characteristics are related to changes in rainfall hazard, storm surge hazard and their dependence.

To address these questions, we apply a full probabilistic joint hazard analysis framework to investigate the current and future joint rainfall–surge hazard from TC events impacting the US Atlantic and Gulf coasts under the combined influence of an end-of-twenty-first century high-emissions scenario SLR (representative concentration pathway (RCP) 8.5)18 and storm climatology change (shared socio-economic pathway (SSP) 5–8.5)19. We generate synthetic TCs from a statistical-deterministic TC model20 forced with reanalysis or GCM output. We downscale 5,018 synthetic TCs consistent with the historical (1980–2005) climate (equivalent to 1,500 simulation years) from National Centers for Environmental Prediction (NCEP) reanalysis data to represent the historical storm climatology. To represent the future storm climatology, 6,200 projected future (2070–2100) TCs are downscaled from each of eight Coupled Model Intercomparison Project phase 6 (CMIP6)19 GCMs, bias corrected and combined into a single weighted average composite projection (for 800 simulation years; Methods). We simulate storm tides (storm surge plus astronomical tide) for each event with the advanced circulation (ADCIRC) hydrodynamic model21,22 using a high-resolution mesh that spans the entire North Atlantic basin and has been previously validated23 (Methods). We estimate rainfall fields using the physics-based Tropical Cyclone Rainfall (TCR) model, which has previously been used to assess historical rainfall climatology24,25, project changes in rainfall hazard26 and to simulate flood impacts27,28 (Methods). To evaluate the impact of SLR, we incorporate spatially varying, probabilistic SLR projections for 2100 from Kopp et al.18, which are based on projections from a suite of CMIP5 GCMs (Methods).

To focus on a particular metric to measure the joint hazard, we define a joint extreme event as one that exceeds both the historical 100-year storm tide (relative to the historical sea level) and the historical 100-year 24-hour rainfall at a given coastal location. On the basis of the TC simulations and bivariate extreme value analysis, we quantify the joint extreme event return period (JRP) in the historical and future climates (Methods) and show that SLR and TC climatology change cause drastic increases in the frequency of joint extreme events. We quantify the relative importance of the change of different climatological variables (that is, sea level, storm frequency, rainfall, storm tides and hazard dependence) in driving the changes in JRP (Methods) and find that TC climatology changes drive larger increases in the joint hazard compared with SLR. We further investigate the effect of TC characteristic changes and find that increases in intensity and decreases in translation speed cause increases in rainfall and surge hazards and their dependence. Our findings motivate explicit consideration of TC climatology changes in compound flood hazard analysis.

Spatial pattern of current and future joint hazard

For each location along the coastline, we calculate the peak storm tide and maximum 24-hour rainfall accumulation occurring anywhere in the upstream catchment for each storm event. On the basis of the NCEP simulations, we quantify the univariate 100-year storm tide (that is, the storm tide level that has a 1% annual probability to be exceeded) and univariate 100-year 24-hour rainfall for the historical period (Supplementary Fig. 1). Using the thresholds of historical 100-year storm tide and rainfall, we quantify the probability of joint extreme event occurrence through JRP in the historical climate (Fig. 1a) and in the future climate (Fig. 1b). We also show the most dominant driver of the JRP change in Fig. 1c. There are large variations in JRP across the US coastline under historical conditions (Fig. 1a). The coastlines of the Gulf of Mexico and southeast Atlantic (up to the Chesapeake Bay) have lower JRPs, typically ranging from 200 to 500 years, signifying a higher probability of joint extreme occurrence compared with other regions. JRP increases along the northern Mid-Atlantic (up to Connecticut) due to a decrease in the statistical dependence between storm tide and rainfall. Along the New England coastline, JRP is much larger than other regions (>1,000 years) because in this region, the two hazards occur almost independently. The low correlation between rainfall and storm tides in New England is due to the large tidal constituents that dominate total extreme sea levels compared with TC-induced storm surges23.

Fig. 1: Joint rainfall–surge hazard in the current and future period and largest driver of joint hazard change.
figure 1

a,b, JRP of 100-year rainfall and 100-year sea level for the NCEP historical period (a) and GCM composite future projection (2070–2100) with 2100 SLR (b). Grey dots in a show representative locations that are analysed further in Figs. 2 and 3. Red tick marks in a show boundaries of the Gulf of Mexico, southeast Atlantic, Mid-Atlantic and New England regions. c, Largest single factor contributing to the overall increase in joint hazard or NA if no single hazard is larger than others. US state outlines were obtained from the US Census Bureau47.

Source data

Due to the combination of future storm climatology change and SLR, future JRPs may decrease to 3–30 years with higher JRP values along the Gulf of Mexico and southeast Atlantic (10–30 years) and lower JRP values along the Mid-Atlantic and New England regions (3–10 years; Fig. 1b). The reason for higher future JRP values along the southern coastline is because these regions are already prone to extreme rainfall and surges in the historical climate (Supplementary Fig. 1) and the percent increase in the hazard there is smaller than the percent increase for northern regions. Thus, across the entire coastline, JRP decreases drastically compared with its historical values. Also, the change in JRP generally increases moving from south (7-36-fold change) to north (30-195-fold change), with the largest decreases in JRP occurring in northern locations. However, even the locations with the smallest JRP changes still correspond to a 7-fold increase in the frequency of joint events. The southeast Florida coast (the Miami region) is an exception to the spatial trend of future JRP. There the historical JRP is 600 years and the future JRP is three years, resulting in a JRP change that is much greater than the JRP change for the rest of the southeast Atlantic. The reason for the large change in JRP in the Miami region is because modelled extreme storm tides and TC rainfall are not highly correlated in the historical period, but large increases in rainfall hazard and SLR in the future cause the joint extreme sea level and rainfall thresholds to be exceeded frequently.

The projection of JRP is associated with statistical and physical modelling uncertainties; Fig. 2 depicts the median JRP estimate (as discussed above) and 95% bootstrapped sampling uncertainty bounds under historical (grey) and composite future (blue) conditions and the JRP estimates from individual GCMs for representative coastal locations. The sampling uncertainty ranges of the composite future JRP (blue boxes) are much smaller than the historical uncertainties because joint exceedances are more frequent in the future period, and consequently, JRP can be estimated with less sampling uncertainty. The variations in JRP estimates among different models are primarily due to differences in the projected future TC frequency and intensity. The Max Planck Institute Earth System Model (MPI), Meteorological Research Institute Earth System Model (MRI) and Geophysical Fluid Dynamics Laboratory Climate Model (GFDL) consistently predict smaller decreases in JRP since these GCMs project low/no increase in storm frequency (Supplementary Fig. 2) and low–moderate increases in storm intensity (Supplementary Fig. 3). Conversely, the EC-Earth Consortium Model (EC-Earth) and Institute Pierre Simon Laplace Climate Model (IPSL) consistently predict large decreases in JRP because both models project the highest increases in storm frequency and intensity. The variations among the GCMs are consistent for the entire coastline (Supplementary Fig. 4). Although there is a relatively large intermodel range of future JRP estimates, especially for locations in the Gulf of Mexico, even the most conservative GCM (that is, MPI) projects large increases in future joint hazard.

Fig. 2: JRP sampling uncertainty and model ranges for representative coastal locations.
figure 2

Median JRP estimate and 95% bootstrapped uncertainty bounds under the NCEP historical (grey) and GCM future composite (blue) forcing. The GCM model ensemble spread at each location for the future period (2070–2100) is shown as coloured dots. CANESM refers to the Canadian Earth System Model, CNRM refers to the Centre National de Recherches Météorologiques and MIROC refers to the Model for Interdisciplinary Research on Climate.

Source data

Drivers of joint hazard change

The change in JRP can be driven by three mechanisms: (1) changes in storm frequency, (2) marginal changes in rainfall totals and/or extreme sea level driven by TC climatology changes and SLR and (3) changes in the statistical dependence between extreme rainfall and storm surges. To understand the relative contribution to changes in JRP from each mechanism, we calculate the isolated impact of changes in storm frequency, rainfall hazard, storm tide hazard, hazard dependence and SLR (Methods). In Fig. 1c we plot the single variable that causes the largest decrease in JRP at each coastal location. Across the Gulf of Mexico and Florida coastlines, the increase in rainfall is the largest driver of changes in JRP, while the increase in storm frequency has the largest impact on JRP changes for parts of the southeast and Mid-Atlantic. Along the upper Mid-Atlantic and New England coastlines, SLR causes the largest decrease in future JRP. For the select locations, we show the relative impact on JRP change of each variable and the combined impact of all storm climatology variables (Fig. 3). Across all locations in Fig. 3 the change in marginal rainfall distribution is among the largest contributor to the change in JRP because all GCMs project large increases in rainfall totals (Supplementary Fig. 5) due to both the increased saturation specific humidity of the warmed environment and the projected increase in TC intensity. In contrast to the large rainfall impact, the change in marginal storm tide distribution has a small impact on the change in JRP for northern locations and a small to moderate impact on JRP change for locations along the Gulf of Mexico. The relative impact of SLR on JRP change generally increases moving south to north, with the largest impact at Portland, ME. Importantly, the storm climatology changes drive large increases in joint hazard across all locations. The combined impact of storm climatology changes on JRP is larger than the SLR impact for 96% of locations along the coastline.

Fig. 3: Relative impact of each climate factor on JRP change and impact of total changes in TC climatology or SLR.
figure 3

Zero indicates no change in JRP compared with the NCEP historical JRP, and one indicates that the factor causes the entire change between historical and future JRP. Negative impact values indicate that the factor increases the JRP compared with the historical best estimate (vertical black lines in Fig. 2). Note that the combined impact of all climate factors on JRP is highly nonlinear and thus the sum of the relative impact of each single factor does not equal one.

Source data

The change in the dependence between hazards also causes a small to moderate decrease in JRP for most locations in Fig. 3, indicating that the extremes of the two hazards are projected to become more dependent in the future climate. To further examine the change in hazard dependence, Fig. 4a shows the conditional probability of 24-hour rainfall exceeding the 90th percentile given a storm tide that exceeds the 90th percentile, calculated for the historical period. The conditional probability is a representation of the tail dependence between the hazards, as higher conditional probability corresponds to higher tail dependence. The eastern Gulf of Mexico and Chesapeake Bay exhibit the strongest dependence between hazards, the western Gulf of Mexico and southeast Atlantic have moderate hazard dependence, and the Mid-Atlantic and New England have relatively low dependence. Figure 4b shows the change in the conditional probability from the historical to future climate, with areas of red (blue) indicating statistically significant increases (decreases) in dependence. With the exception of the eastern Gulf of Mexico, Chesapeake Bay and Maine coastlines, most regions are projected to have higher dependence between extreme rainfall and storm tides in the future. Specifically, the lower Texas, Georgia, North Carolina and New Jersey coastlines are projected to experience the largest strengthening of hazard dependence in the future, resulting in an increase of up to 0.2 in the conditional probability (Fig. 4b). Along the eastern Gulf of Mexico, there is almost no change in the dependence strength because the two hazards are already highly correlated in the historical climate (Fig. 4a) and will remain similarly correlated in the future climate. Along the coast of Maine, there is a small projected increase in hazard dependence, although this increase is not statistically significant. The Chesapeake Bay stands as an outlier, and it is the only location where the dependence strength between hazards decreases in the future climate (discussed below).

Fig. 4: Historical and future change in tail dependence between 24-hour rainfall and peak storm tide.
figure 4

a, Conditional probability (P) of extreme rainfall (exceeding 90th percentile; R) given extreme storm tide (exceeding 90th percentile; S) in the historical period. b, Change in conditional probability (ΔP) of extreme rainfall given extreme storm tide due to future storm climatology. Positive (negative) values indicate increases (decreases) in conditional probability. Areas of grey indicate that the projected change in conditional probability is not significant compared with the range of natural variability in the historical period (set as the 10th–90th percentiles of the tail dependence estimated through bootstrapping). US state outlines were obtained from the US Census Bureau47.

Source data

Changes in dominant TC storm characteristics

Since TC climatology change is the dominant contributor to JRP change, we investigate how projected changes in TC storm characteristics drive changes in rainfall accumulations, peak storm surges and their dependence at the coast. After investigating correlations between each hazard and storm intensity, approach angle, translation speed and landfall location and quantifying projected changes in each storm characteristic, we find that storm intensity and translation speed are both projected to change greatly in the future (Fig. 5a,b, respectively) and are significantly correlated with rainfall and/or storm tide (Fig. 5c–f). For the vast majority of the coastline, both the peak storm tide and 24-hour rainfall are significantly correlated with TC intensity, although the strength of correlation is higher for rainfall (Fig. 5c,d). The 24-hour rainfall is also strongly negatively correlated with storm translation speed (Fig. 5f) as slower-moving storms will drop more rainfall in a given coastal location than faster-moving storms. The peak storm tide is not strongly correlated with translation speed (Fig. 5e) because both slow- and fast-moving storms can generate high surges, and the additional background wind contribution is generally small, even for fast-moving storms, compared with the cyclonic wind speed. Under future storm climatology, the 90th percentile of TC intensity is projected to increase by 15–30% along the Gulf of Mexico and southeast Atlantic, 30–50% along the Mid-Atlantic and 20–30% along the New England coastlines (Fig. 5a). The vast majority of previous studies also project an increase in North Atlantic TC intensity, and many predict an increase in the frequency of high-intensity (category 3–5) TCs29. We also find a large future reduction in the translation speed of storms that exceed the 90th percentile intensity (Fig. 5b). For all regions except New England, storms that exceed 90th percentile intensity are projected to move 20–30% slower in the future than in the historical period. The decrease in translation speed found here is consistent with previous work examining changes in translation speed in the historical record30 and projections of TC translation speed under future climate conditions31,32,33. The increase in storm intensity coupled with the decrease in translation speed drives an increased likelihood to observe both extreme rainfall and extreme storm tide in the future and increases the upper tail dependence between the hazards. By comparing Fig. 4b with Fig. 5a,b, it is clear that most regions experiencing a significant increase in the hazard dependence also experience large increases in storm intensity and decreases in translation speed. The Chesapeake Bay is a notable exception because the hazards are projected to become less dependent in the future even though there is an increase in TC intensity and decrease in translation speed. In the future, a larger number of intense storms are projected to approach the coast north of the Bay opening. These storms do not induce high storm surges within the Bay since the cyclonic winds are pointed away from the coast but they still induce extreme rainfall. Thus, the increase in the number of these types of storm causes a decrease in the hazard correlation at this location in the future climate.

Fig. 5: Change between future composite TC characteristics and historical characteristics and correlation between rainfall/storm tide and TC characteristics.
figure 5

a,b, Change in 90th percentile TC intensity (Vmax) (a) and median translation speed (Vt) (b) of storms that exceed 90th percentile intensity. cf, Kendall correlation (Ck) between Vmax and storm tide (c) or rainfall (d) and between Vt and storm tide (e) or rainfall (f). US state outlines were obtained from the US Census Bureau47.

Source data

Discussion

The results presented here demonstrate that TC climatology change and SLR may cause large increases in joint rainfall–surge hazard across the US East and Gulf coasts. The projected increase in extreme rainfall hazard (considering the maximum 24-hour rainfall accumulation over the catchment in the above analysis) is often the largest driver of the increase in the extreme joint hazard. Our projections of extreme rainfall are consistent with the work of Emanuel26, who found a 100–120% increase in the 100-year storm total rainfall at a single point location in Houston, TX, while we project a 123% increase (Supplementary Table 2). Our projections are also consistent with previous studies focusing on mean rainfall changes. Using the RCP 4.5 scenario and a suite of CMIP5 models, most previous studies found a 10–22% increase in mean end-of-century TC rain rates within 100 km of the storm centre34,35,36. A recent study using a high-resolution GCM projected a larger increase of 29%37. Here we project a slightly higher 32% increase in inner core mean TC rain rate (Supplementary Table 1), which is reasonable given our use of the SSP5–8.5 high-emission scenario.

We also find that the overall impact of storm climatology change on the change in the extreme joint hazard is larger than the SLR impact for 96% of the coastline. The contribution of TC climatology change is also dominant for lower joint TC hazard thresholds, such as 25-year or 50-year levels (Supplementary Table 3 and Supplementary Fig. 6). Although we find that TC climatology change is more dominant than SLR in driving changes in TC joint hazard, SLR also impacts other types of compound flooding arising from, for example, ETCs or two unrelated meteorological events, especially for return periods shorter than 50 years13,14,15,16. Moreover, 2017 work by Kopp et al.38 that incorporated a physical model for ice sheet hydro-fracturing, a mechanism that is deeply uncertain, found significantly higher SLR by 2100 than Kopp et al. in 201418 (which we use here). Therefore, the overall role of SLR on total compound flood hazard may still be dominant compared with TC climatology change.

The findings presented here are associated with inevitable uncertainties. We use a single TC model to downscale all GCMs and reanalysis data, and the model predicts an increase in future TC frequency for five of the eight GCMs. Although a few other studies39,40 have also predicted increases in TC frequency, the majority of studies predict a decrease or no change in North Atlantic storm frequency29. However, the main findings of our study are unchanged even if we assume no change in future TC frequency. The future JRP change calculated by holding TC frequency constant at the historical level is only slightly lower at each coastal location (up to 149-fold decrease in JRP; Supplementary Fig. 7), and the spatial trends (that is, higher JRP change in the north compared with the south) are unchanged. The relative importance of TC climatology change compared with SLR also remains similar when assuming constant frequency, and TC climatology change still causes a larger JRP change than SLR for 84% of the coastline. The reason our results are relatively unchanged if we neglect the projected frequency change is because the increase in TC hazards and their joint occurrence is largely driven by projected increases in TC intensity and decreases in translation speed.

This study cannot directly predict the overall compound flood hazard, which is driven by a combination of ETC events (especially at lower return periods) and TCs. Moreover, compound flood depths must be quantified using high-resolution inundation models. Nevertheless, we provide evidence that joint rainfall–surge extreme events could become an increasing threat to coastal communities in the future. We also find that the statistical dependence between extreme rainfall and storm tide increases in the future for portions of the coastline, resulting in a higher probability of multi-hazard extremes during future storm events. This finding is important because many previous studies of future compound flooding have neglected potential increases in hazard dependence8,9,10,41, which could underestimate compound flood risk. Our projections of joint TC rainfall–surge hazard can be combined with ETC hazard distributions42 to develop overall flood-mapping scenarios43 for regional-10,44 or local-scale17,45,46 flood models to assess the impact of joint rainfall–surge occurrence on coastal flooding in a changing climate.

Methods

To characterize the present and future joint rainfall–surge hazard, we implement a physics-based modelling framework that is driven by the large-scale atmospheric and ocean climatology of reanalysis (historical period) or GCM (future period) data. First, we construct monthly climatologies of relevant environmental variables48 based on the reanalysis/GCM data. Next, we generate thousands of synthetic TCs that are consistent with the large-scale environment using a statistical-deterministic TC model. These synthetic TCs represent around 1,000 simulation years for each climate condition. For each TC we model the coastal storm tides using a high-resolution hydrodynamic model, and we model the rainfall fields using a computationally efficient physics-based rainfall model. On the basis of the modelled storm tides and rainfall accumulations for the thousands of synthetic TCs, we conduct bivariate statistical analysis to quantify the probability of joint extreme events.

Data

We generated 5,018 synthetic TC tracks for the historical time period (between 1980 and 2005) based on the NCEP reanalysis49. We then generated 4,400 synthetic TCs for the historical period (1984–2005) and 6,200 TCs for the future period (2070 to 2100) under the SSP5–8.5 emissions scenario19 based on each of eight CMIP619 climate models: Canadian Earth System Model (CANESM), Centre National de Recherches Météorologiques (CNRM), EC-Earth Consortium Model (EC-Earth), Geophysical Fluid Dynamics Laboratory Climate Model (GFDL), The Institute Pierre Simon Laplace Climate Model (IPSL), Model for Interdisciplinary Research on Climate (MIROC), Max Planck Institute Earth System Model (MPI) and Meteorological Research Institute Earth System Model (MRI).

Synthetic TC model

The statistical-deterministic TC model20, which has been widely applied for TC hazard assessment27,50,51,52,53, generates synthetic events based on data about the large-scale environment and can be forced with either reanalysis data or projections from GCMs. Vortices are randomly seeded in space and time and are moved according to the large-scale environmental winds plus a beta-drift correction54. TC intensity is estimated at each time step based on the Coupled Hurricane Intensity Prediction System (CHIPS), which is an axisymmetric vortex model coupled to a 1D ocean model55. Storms are retained only if their intensity exceeds 21 m s−1 (40 knots). Thus, only seed vortices that encounter favourable large-scale environmental conditions will strengthen into TCs, and the timing of TC development is consistent with the environmental climatology. For each TC, the outer radius at which the cyclonic wind speed goes to zero (henceforth outer radius) is randomly drawn from an empirical lognormal distribution56. We neglect the variation in outer radius size over the TC lifetime57 because previous work has shown the outer radius variation to be relatively small58. We also assume no change in the distribution of TC outer size for the future climate because historical trend analysis for the North Atlantic basin found no statistically significant changes in TC size over time59. Moreover, an analysis of dynamically downscaled TCs based on RCP 4.5 end-of-century forcing found nearly constant outer radii compared with the historical period36. Using the CHIPS-estimated intensity and outer radius, we estimate the radius to maximum winds based on a theoretical wind model that links the outer descending region of the TC with the inner ascending region58. Each simulated storm is characterized by a time series of storm parameters (time, centre position, maximum wind speed, pressure deficit and radius to maximum wind) for every two hours.

Bias correction and model combination

The downscaled TCs from each GCM may be biased compared with the NCEP-downscaled TCs, and biases within the TC characteristics can propagate to become biases in the hazard estimation. TC intensity and annual frequency are both important drivers of coastal flood risk, and both variables may be biased due to biases in GCM projections. Therefore, we perform bias correction at the storm level based on the differences between the NCEP TC frequency and intensity distribution and the GCM-predicted frequency and intensity distribution for the historical period. Using our method of bias correction, we avoid multivariate bias correction on the modelled storm tides and rainfall, which often fail to preserve the entire dependence structure between hazards60. Additionally, bias correction at the storm level is computationally efficient, while bias correction at the hazard level requires performing intensive hydrodynamic simulations for additional thousands of GCM TCs for the historical period.

Specifically, at each location we bias correct the TC frequency by multiplying the GCM-predicted future frequency by the ratio of the NCEP-derived historical frequency and GCM-predicted historical frequency. To correct the GCM-projected TC intensity (Vmax) of each storm set, we first use the quantile delta mapping approach described by Cannon, Sobie and Murdock61 applied to each location along the coast. Essentially, the change between the GCM-projected future (2070–2100) and historical (1984–2005) downscaled Vmax quantiles is added to the NCEP-downscaled historical quantiles to create a corrected future Vmax distribution for each GCM model at each location. Then, by the principle of importance sampling62, the GCM-projected storms are weighted and re-sampled with weights corresponding to the ratio of the corrected Vmax probability density to the GCM-projected Vmax probability density. By doing weighted re-sampling of the storms at each location, we are able to match the corrected future Vmax distribution and consequently generate a storm set at each location that is unbiased with respect to the intensity distribution. Supplementary Fig. 8 shows the bias correction procedure applied at a sample location for a sample GCM, demonstrating that after weighting/re-sampling, the target Vmax distribution is matched accurately. We also create a composite projection for the future climate using a weighted average across all GCM storm sets, where the weights of each GCM are based on their Willmott skill63 in matching the NCEP TC intensity return level curve in the historical period (Supplementary Fig. 9).

Hydrodynamic modelling

We simulate TC storm tides using the 2D depth-integrated version of the ADCIRC model21,22. We use an unstructured computational mesh developed by Marsooli and Lin23 that spans the entire North Atlantic basin and has resolution varying from >50 km in the deep ocean to ~1 km near the coastline. Eight tidal constituents are incorporated as periodic boundary conditions at the ocean boundaries of the mesh, and tidal data are obtained from the global model of ocean tides, TPXO8-ATLAS64. The timing of the tide is matched to the timing of the synthetic storm (simulated according to the climatology). Wind and pressure fields are developed based on the Vmax and radius to maximum wind (Rmax) of each synthetic TC and physics-based parametric models65,66. Further details regarding the mesh formulation, tidal forcing and wind/pressure models are documented by Marsooli and Lin23. Simulated storm tides from the model configuration used in this study were compared against observed water levels for 191 historical TCs impacting the US East and Gulf coasts, and the model was found to satisfactorily reproduce peak storm tides (with an average root mean square error of 0.31 m and Willmott skill of 0.90)23. In this study we do not account for wave setup because the computational expense of coupled wave-surge model would prevent a large-scale Monte Carlo risk assessment. For each TC we extract peak storm tides at nodes along the coastline that are spaced roughly 25 km apart.

Rainfall modelling

We estimate rainfall fields from each synthetic TC using the TCR model described by Zhu, Quiring and Emanuel24. TCR is a physics-based model that simulates convective TC rainfall by relating the precipitation rate to the total upward velocity within the TC vortex. Vertical velocity is estimated by taking into account frictional convergence, topographic forcing, vortex stretching, baroclinic effects and radiative cooling. TCR has been previously used in risk assessment studies26,27 and was recently compared against observed TC rainfall across the United States25. Xi, Lin and Smith found25 that TCR simulates the rainfall climatology of coastal regions with relatively good accuracy, although it underperforms in inland and mountainous regions. The performance of the model for inland regions has been addressed and improved in subsequent work53. TCR does not simulate outer TC rain bands, which are 3D in nature and cannot be directly simulated with an axisymmetric model. Nevertheless, a recent study modelled compound flooding using TCR-predicted rainfall fields for several historical events and found that TCR rainfall produced similar flood depth/extent compared with using radar rainfall forcing27. In our study, we use TCR rainfall over each coastal catchment delineated according to US Geologic Survey hydrologic units67. We pair each coastline point with its upstream coastal catchment, and for the coastal point, we use the maximum 24-hour rainfall accumulation occurring anywhere in the upstream catchment as our rainfall metric for each storm event. The 24-hour storm duration is frequently used for rainfall risk assessment studies68, and rainfall occurring anywhere in the immediate upstream catchment will drain to the same coastal point and can impact compound hazard.

Validation of integrated modelling of TC surge–rainfall hazard

Previous studies have independently evaluated the TC model20,48, rainfall model2553 and storm tide model23 by comparing against historical observations. Here we additionally evaluate the ability of our models to reproduce observed dependence between TC rainfall and storm tides. We compare the Kendall rank correlation69 computed from modelled rainfall and storm tides (derived from reanalysis data) against the Kendall correlation computed from observed storm tides and observed daily rainfall at 31 gauge locations across the coastline (Supplementary Fig. 10). The Kendall correlation coefficient can capture nonlinear dependence between two variables by using the relative ranks of each observation rather than the magnitude, and Kendall correlation has been used extensively as a metric to assess dependence between rainfall and storm tides1,70,71. If the modelled rainfall and storm tides from the NCEP synthetic TCs produce a correlation coefficient similar to the observations, this suggests that the models produce joint high (and joint low) events with similar likelihood as the real observed TCs and thus increase our confidence in the use of our models to project current and future joint hazard.

On the basis of Supplementary Fig. 10, the model-based correlations match well with the observed correlations with an overall root mean square error of 0.09 and bias of 0.02 (indicating slight overestimation of rainfall–surge dependence). For the majority of locations, the difference between modelled and observed correlations is within ± 0.1. The model overestimates the correlation for the region between Mississippi and the Florida panhandle. The discrepancy between modelled and observed correlation in this region is probably due to the occurrence of other observed rainfall mechanisms such as extra-tropical transition or interaction with fronts that are not simulated by the TC model and cause lower correlation between observed rainfall and storm tides.

SLR projections

We incorporate probabilistic, localized SLR projections from Kopp et al.18 for 2100 considering the RCP 8.5 emissions scenario. In that study18, SLR probability distributions are developed for tide gauge locations across the globe by taking into account ice sheet components (Greenland, West Antarctic and East Antarctic), glacier and ice cap surface mass balance, thermal expansion and oceanographic processes, water storage on land and other non-climatic factors. Sea-level changes due to thermal expansion and oceanographic processes are based on ensemble mean projections from a suite of CMIP5 GCMs. For each point along the coastline, we select the nearest tide gauge and adopt the probability distribution specified by Kopp et al.18.

We calculate total sea level for each TC by randomly drawing from the SLR distributions and superimposing on the modelled storm tides for computational efficiency. The assumption of linearity between SLR and storm tides is a reasonable approximation of extreme sea levels, but nonlinear interactions between SLR and storm tides can be significant in complex local areas, particularly small bays and estuaries72,73. We also treat TC climatology change and SLR as independent, although they may be significantly correlated. Little et al.74 found a significant correlation between SLR and changes in the power-dissipation index (an integrated measure of TC intensity, frequency and duration) for the North Atlantic, suggesting that large increases in mean sea level are more likely to co-occur with large increases in TC hazard. By neglecting correlation between SLR and climatology changes, our results may underestimate the composite (weighted average) change in climatology and SLR and consequently represent a conservative estimate of joint hazard change.

Statistical analysis of joint hazard

We conduct statistical analysis on the pairs of maximum modelled storm tides (or storm tides plus SLR) and maximum 24-hour rainfall accumulation at each location along the coastline to quantify their marginal and joint hazard.

The marginal distributions of both rainfall and storm tides are often characterized by a long tail representing the rare but extreme events50,51. The heavy tail can be modelled with a peaks-over-threshold approach, where the probability above a high threshold is estimated by a Generalized Pareto (GP) distribution75. We fit marginal GP distributions using the maximum likelihood method75 for the rainfall and storm tides at each location, and the threshold is set by numerically minimizing the root mean square error between the empirical quantiles and the theoretical quantiles. According to bivariate extreme value theory, a logistic model can be used to estimate the joint distribution of two GP variables such that their bivariate cumulative distribution function (CDF) represented as G(x,y) takes the form75,76:

$$G(x,y) = \exp \{ - (\tilde x^{ - 1/\alpha} + \tilde y^{ - 1/\alpha })^\alpha \}$$
(1)

where \(\tilde x\) and \(\tilde y\) are the Fréchet-transformed versions of the variables x and y, and α is a parameter that quantifies the strength of the dependence between the variables (α → 0 signifies complete dependence and α = 1 signifies complete independence). At each location we transform the rainfall and storm tide pairs based on their respective marginal distributions and GP thresholds to obtain Fréchet versions of the variables. Then we fit the bivariate distribution using a censored maximum likelihood approach76 that considers pairs that jointly exceed their GP thresholds (within the ‘evd’ R-package77). We additionally ensure that there are at least 20 pairs of joint exceedances to fit the bivariate model. The bivariate logistic model employed here has previously been used to model dependence between rainfall and storm surges2,76,78,79.

After characterizing the marginal and joint distributions of rainfall and storm tides at each coastal location, we quantify the return period (inverse of the annual exceedance probability) of joint extreme events. For each location, we model TC occurrence as a Poisson process with arrival rate λ per year. The basin arrival rate is a parameter of the TC model20 and is calibrated to match the observed occurrence rate in the North Atlantic basin for the historical period. The location-specific arrival rate (λ) is an adjustment of the basin arrival rate according to the proportion of storms passing within 200 km of each location. We define xT, yT as the marginal 100-year storm tide and 100-year rainfall defined in the historical period. Then the return period of an event that jointly exceeds xT and yT (JRP) is calculated as follows:

$$\mathrm{JRP} = \frac{1}{{1 - e^{ - \lambda P}}}$$
(2)

where P is the joint exceedance probability:

$$P = 1 - \Pr (X \le x_\mathrm{T}) - \Pr (Y \le y_\mathrm{T}) + G(x_\mathrm{T},y_\mathrm{T})$$
(3)

Pr(.) is the probability operator and G is defined in equation (1).

We quantify JRP under the current and future storm climates by fitting marginal and joint distributions to storm tide and rainfall pairs from NCEP or each GCM-derived storm dataset. We estimate the sampling uncertainty bounds of the JRP estimates by implementing a bootstrapping approach with 500 iterations for each location and each GCM. For each iteration, we re-sample (with replacement) pairs of modelled storm tides and rainfall, fit the univariate and joint distributions and re-calculate JRP.

Attribution of changes in joint hazard

To quantify the isolated impact of various climate factors on changes in joint rainfall–surge hazard, we adjust a single factor at a time and then re-calculate JRP. To quantify the isolated impact of SLR on changes in JRP, we randomly draw SLR values from location-specific probability distributions18 and add them to the historical rainfall–storm tide pairs. The impact of changes in future storm frequency is quantified by simply changing the value of λ in equation (2) to reflect the future period frequency. Because storm tide and rainfall are dependent, we quantify the impact of changes in (1) marginal rainfall distribution, (2) marginal storm tide distribution and (3) dependence between hazards through quantile matching. Specifically, we calculate Fr,h and Fs,h, which are the historical rainfall (rh) and storm tide (sh) CDFs, and Fr,f and Fs,f, which are the future CDFs. Given historical pairs of rainfall and storm tide (rh, sh), we can evaluate the impact of changes in rainfall hazard by changing rh values to \(r_h\!^\ast = F_\mathrm{r,f}\!^{ - 1}(F_\mathrm{r,h}(r_\mathrm{h}))\) so that the magnitude of rainfall is increased according to the future period rainfall distribution but the sh values and dependence between hazards are unchanged. We similarly calculate the storm tide values (\(s_\mathrm{h}\!^ \ast\)) while keeping the rainfall values (rh) constant to evaluate the impact of increases in storm tide on the JRP change. The methodology above guarantees the rank correlation between TC rainfall and surge is unchanged. To measure the impact of changes in hazard dependence (α in equation (1)), we adjust the future rainfall and storm tide pairs (rf, sf) as follows: \(r_\mathrm{f}\!^ \ast = F_\mathrm{r,h}\!^{ - 1}(F_\mathrm{r,f}(r_\mathrm{f}))\), \(s_\mathrm{f}\!^ \ast = F_\mathrm{s,h}\!^{ - 1}(F_\mathrm{s,f}(s_\mathrm{f}))\). The adjusted values of rainfall and storm tide are reduced according to their historical distributions, but the dependence between hazards is based on the future period climatology.