## Historical Observations

TCs that made landfall in the United States from 1949–2018 were used for historical trend analysis. The information on historical TCs was obtained from the International Best Track Archive for Climate Stewardship (IBTrACS38), which provides six-hourly TC locations and intensities. To quantify and examine the change of impacts from historical landfalling TCs, observations of hourly water levels, daily rainfall, and daily maximum wind at 9 locations across the U.S. coastal areas (location shown in Fig S1) were obtained from the Center for Operational Oceanographic Products and Services (CO-OPS) and National Centers for Environmental Information (NCEI). If the tidal gauge sites had no information for rainfall and wind, we found the closest weather stations to the tidal gauge location and used the rainfall and wind observations in these stations to represent the wind and rainfall hazards at the tidal gauge station.

We found that there was more missing data of wind observation from 1949–1979 than in later decades while the missing data for surge and rainfall observations had no significant decadal variations. To eliminate the possibility that the trends were caused by the more frequent missing wind observations in the earlier decades, similar analysis was conducted for only surge and rainfall hazards, and the same conclusion was reached. The fact that using only surge and rainfall hazards can obtain a similar conclusion as using the triple hazards does not imply that the wind hazard is not important; rather, it implies that the three hazards are physically correlated.

**Synthetic TCs**

Synthetic TCs generated with a statistical-deterministic TC model15,25 were used to study the change of sequential hazard-producing TCs in the United States. The synthetic storms used in this study are the same as in ref6. The storms were generated under the environment of six CMIP6 climate models (CanESM5, CNRM-CM6-1, GFDL-CM-4, EC-Earth3, IP-SL-CM6A-LR, and MIROC6) for the control (1984-2005) and Shared Socioeconomic Pathway 5 8.5 scenarios (SSP 5 8.5; 2070-2100); 4400 and 6200 U.S. landfalling storms were generated for each model under each scenario, respectively.

**Wind Modeling**

The complete wind profile model that merges the inner core and outer radii TC wind profiles (C1519) was used to simulate the wind hazard associated with the synthetic TCs. To obtain the surface wind affected by environment winds, a correction was added to the simulated wind profile following ref39. The C15 model was also used to prepare the required wind input for the TC rainfall simulation.

**Hydrodynamical Modeling**

The ADCIRC model was used to simulate the storm tides produced by the synthetic TCs. The unstructured computational mesh developed by ref40, which has a spatial coverage of the entire North Atlantic basin, is used in this study. Eight tidal constituents41 are applied as boundary conditions at the ocean boundary of the mesh. The wind and pressure fields associated with synthetic TCs are needed for the storm tide simulations and were obtained using physics-based parametric models42. This wind model was used in ref6 to drive surge simulation as it is simpler and has similar performance to the C15 wind model. Further details regarding ADCIRC model and simulation setups can be found in ref40 .

**Rainfall Modeling**

The physics-based TC rainfall model (TCR12) was used to simulate rainfall associated with the synthetic TCs. Detailed formulation on the TCR can be found in ref (18). The simulation setup of the model followed ref43, and the environmental field needed for the TCR simulation was obtained following ref12. The C15 model was used to drive the TCR rainfall simulation following ref6 .

**Sea Level Rise (SLR) Projection**

Localized probabilistic SLR projection20 under the RCP 8.5 emission scenario was incorporated in this analysis. The projection of SLR20 was developed for tidal gauge locations; for each point along the coastline, we selected the nearest tide gauge to represent the projected change of SLR in the future. We followed ref20 to generate multiple SLR values from 2070 to 2100 and randomly drew TC events in each year to combine them with the random SLR value in the corresponding year.

**Joint Hazard Analysis**

Statistical analysis on the triple of the modeled maximum storm maximum water level (L), maximum daily rainfall (R), and maximum wind speed (W) across each coastal segment of Texas, Louisiana, Mississippi-Alabama, West Florida, East Florida, Georgia, South Carolina, and North Carolina11 were performed. The marginal distributions of rainfall, surge, and wind were fitted by generalized Pareto distributions (GPDs) to characterize the long tails that corresponding to the extreme events6,29. The GPDs were fitted at each coastal segment, and the threshold was set by minimizing the mean squared error between empirical quantiles and the theoretical quantiles5. There was no parametric probabilistic distribution that described the trivariate GPD distributions, so nested Gumbel Copulas44 were used to represent the dependent structure of the three individual hazards. Gumbel Copulas are used here because previous studies6,7,45 show that each pair of hazards (surge-rainfall, surge-wind, and wind-rainfall) are correlated especially at tails, and Gumbel Copulas are often used to quantify the tail-dependent structure46.

After fitting the marginal and joint distributions of the triple of TC hazard, we calculated the probability of three hazards jointly exceeding respective thresholds (joint exceedance probability, or JEP) and the probability of at least one of the three hazards exceeding the respective threshold (at least one exceedance probability, OEP). The thresholds \({L}_{T}\), \({R}_{T}\), and \({W}_{T}\) are the marginal return levels of surge, rainfall, and wind hazard with a return period of T years in the control simulation. Thus, JEP and OEP are mathematically defined as:

$$JEP\left(T\right)=\mathbb{P}\left(L>{L}_{T}\cap R>{R}_{T}\cap W>{W}_{T}\right)$$

1

and

$$OEP\left(T\right)=1-\mathbb{P}(L\le {L}_{T}\cap R\le {R}_{T}\cap W\le {W}_{T})$$

2

To explicitly discuss the effect of SLR, the JEP and OEP under the SSP 5 8.5 scenario were calculated with and without the consideration of SLR. To explore the causation of the change of JEP and OEP, the changes of the exceedance probability of the single hazards (the marginal cumulative density function) were also investigated. The return period (\({T}_{JE}\)) of the three hazards jointly exceeding respective return level thresholds (under return level \(T\)) can be calculated as

$${T}_{JE}\left(T\right)=\frac{1}{\lambda \bullet JEP\left(T\right)}$$

3

where \(\lambda\) is the arrival rate of the storms.

Similarly, the return period \({(T}_{OE})\)of at least one of the three hazards exceeding the respective threshold can be calculated as

$${T}_{OE}\left(T\right)=\frac{1}{\lambda \bullet OEP\left(T\right)}$$

4

**Probabilistic Model for Sequential Hazard-Producing TCs**

The Poisson-Gaussian model for sequential landfalling TCs developed in ref11 was extended to capture the storms’ hazard-producing ability and hazard durations. The arrival of hazard-producing storms was modeled as a nonstationary Poisson process as

$$\nu \left(t,s\right)={\lambda }_{hazard}\left(t\right){S}_{hazard}\left(s\right)$$

5

$${S}_{hazard}\left(s\right)= \frac{1}{{\sigma }\sqrt{2{\pi }}}{\text{e}}^{-\frac{1}{2}{\left(\frac{\text{s}-{\mu }}{{\sigma }}\right)}^{2}}$$

6

where \({\lambda }_{hazard}\left(t\right)\) is the annual frequency of hazard-producing TCs in year t. \({S}_{hazard}\left(s\right)\), representing the seasonal variation, is the likelihood of a hazard starting to occur on days relative to the likelihood of occurrence during the season. The difference between the Poisson-Gaussian model in this study and ref11 is that this study uses the annual frequency and seasonality of hazard-producing storms instead of landfalling storms, and the ratio between the annual frequency of hazard-producing storms and landfalling storms can be viewed as a metric of storm hazard-producing ability. The Poisson-Gaussian model was shown to be capable of capturing the relationship between TC landfall climatology and the minimal landfall intervals11, and larger frequency (\({\lambda }_{hazard}\)) and smaller seasonal span (\({\sigma }\)) favors sequential hazard-producing events.

We define the MII mathematically as:

$$\text{M}\text{I}\text{I}\left(\text{t}\right)=\underset{i>j}{\text{min}}(\text{m}\text{a}\text{x}({B}_{i}\left(\text{t}\right)-{E}_{j}\left(\text{t}\right),0\left)\right)$$

7

where \({B}_{i}\)(t) is the hazard beginning time of the i-th TC in year t, and \({E}_{j}\left(\text{t}\right)\) is the hazard end time of the j-th TC in year t; if the impact of two TCs overlap, we set the MII equal to 0. The TCs in a single year were ordered by the starting time of hazards associated with the TCs. For each individual TC, the beginning times of hazards were randomly drawn from the probabilistic model of Eq. (5) and Eq. (6), and the beginning and end times of hazards were connected as:

$${E}_{i}\left(\text{t}\right)= {B}_{i}\left(\text{t}\right)+{D}_{i}\left(\text{t}\right)$$

8

where \({D}_{i}\left(\text{t}\right)\) is the total hazard duration of the i-th TC in year t. As the impact duration of TC-related hazards may change in the future, a probability distribution that describes the duration of TC hazard impact at each selected location is needed. The probability distributions of duration in different coastal locations under different climate scenarios do not share the same parametric probability distribution. Thus, in our analysis and simulation, we applied kernel density estimation to fit the non-parametric probability distribution to the duration. For climate simulation analysis, the explicit simulation of each hazard component was performed before the probabilistic model was fitted with the simulated synthetic events.

**Definition of Hazard-Producing TCs and Hazard Duration**

There is no universal definition for “hazard-producing” as different infrastructure systems and communities may respond differently to sequential TC events11. A statistically reasonable, physically meaningful, and engineeringly applicable definition of “hazard-producing” should be able to both categorize the landfalling storms into “hazard-producing” and “non-hazard-producing” categories and separate the days of impact from a single “hazard-producing” TC into “hazard days” and “non-hazard days” (such as the days when TCs are too weak or too far away from the point of interest to produce hazards).

To do so, specific percentiles of daily maximum water level, total rainfall, and maximum wind speed were used as the thresholds to define “hazard-producing” TCs. In the historical analysis, for each individual point of interest, the daily maximum water level, total rainfall, and maximum wind from every TC that ever approached 250 km from the location were collected. and the percentiles of each hazard component were calculated based on the observation of all TCs impacting this location. For the climate simulation using CMIP6 models and the synthetic storm model, the same method was applied to calculate percentiles for the control simulation (historical period) of each climate model (so the percentiles differ for the climate models).

The specific percentile chosen as the threshold should be both high enough to eliminate nuisance TC-events and low enough to include some “non-extreme” events that are still hazardous. In this main text of this study, the 95th percentile of each hazard was chosen as the threshold, and the days when at least one hazard component exceeds the threshold were defined as “hazard days.” The threshold was spatially varied by using the 95th percentiles specific to each coastal location to account for spatial variation of the “preparedness” or “awareness” of hazards. The storms that cause at least one “hazard day” for a point of interest were defined as “hazard-producing” storms for this location. The selection of the 95th percentile of threshold fulfills the requirements mentioned in the beginning of this paragraph. For example, the mean of the 95th percentile of 6 climate model simulations of daily accumulated TC rainfall (maximum tide level) in the control climate is 98.8 mm (1.4 m) near New Orleans, which is approximately ¼ (1/6) of the total rainfall (maximum tide level) that Hurricane Katrina (2005) produced in this location. Admittedly, the selection of the 95th percentile was ad hoc, so the results of climate simulations of sequential TC-related hazards under other thresholds (90th and 99th percentiles) are shown in the Supplementary Material as a sensitivity test of the results of this study.