Tropical cyclone climatology change greatly exacerbates US joint rainfall-surge hazard

11 Tropical cyclones (TCs) are among the largest drivers of extreme rainfall and surge, but 12 current and future TC joint hazard has not been well quantified. We utilize a physics-based 13 approach to simulate TC rainfall and storm tides and quantify their joint hazard under 14 historical conditions and a future (SSP5 8.5) climate projection. We find drastic increases in 15 the frequency of exceeding joint historical 100-yr hazard levels by 2100, with a 10-36 fold 16 increase along the southern US coast and 30-195 fold increase in the northeast. The joint 17 hazard increase is induced by sea-level rise and TC climatology change; the relative 18 contribution of TC climatology change is higher than that of sea-level rise for 96% of the 19 coast due to large increases in rainfall. Increasing storm intensity and decreasing 20 translation speed are the main TC change factors that cause higher rainfall and storm tides 21 and up to 25% increase in their dependence.

6 changes drive large increases in joint hazard across all locations. The TCC impact on JRP is 157 larger than the SLR impact for 96% of locations along the coastline. 158 The change in the dependence between hazards also causes a small to moderate 159 decrease in JRP for most locations in Figure 3, indicating that the extremes of the two 160 hazards are projected to become more dependent in the future climate. To further examine 161 the change in hazard dependence, Figure 4a shows the conditional probability of 24-hour 162 rainfall exceeding the 90 th percentile given a storm tide that exceeds the 90 th percentile, 163 calculated for the historical period. The conditional probability is a representation of the 164 tail dependence between the hazards, as higher conditional probability corresponds to 165 higher tail dependence. The eastern Gulf of Mexico and Chesapeake Bay exhibit the 166 strongest dependence between hazards, the western Gulf of Mexico and Southeast Atlantic 167 have moderate hazard dependence, and the Mid-Atlantic and New England have relatively 168 low dependence. Figure 4b shows the composite change in the conditional probability from 169 the historical to future climate, with areas of red (blue) indicating increases (decreases) in 170 dependence. With the exception of the eastern Gulf of Mexico and Chesapeake Bay, most 171 regions are projected to have higher dependence between extreme rainfall and storm tides 172 in the future. Specifically, the lower Texas, Georgia, North Carolina, and New Jersey 173 coastlines are projected to experience the largest strengthening of hazard dependence in 174 the future, resulting in up to 0.2 increase in the conditional probability (Fig. 4b). Along the 175 eastern Gulf of Mexico there is almost no change in the dependence strength, which is 176 because the two hazards are already highly correlated in the historical climate (Fig 4a) and 177 will remain similarly correlated in the future climate. The Chesapeake Bay stands as an 178 outlier, and it is the only location where the dependence strength between hazards 179 decreases in the future climate (discussed below). 180 181

Changes in dominant TC storm characteristics 182
Since TC climatology change is the dominant contributor to JRP change, we investigate how 183 projected changes in TC storm characteristics drive changes in rainfall accumulations, peak 184 storm surges, and their dependence at the coast. After investigating correlations between 185 each hazard and storm intensity, approach angle, translation speed, and landfall location 7 intensity and translation speed are both projected to change significantly in the future ( Fig  188   5a and 5b, respectively) and are significantly correlated with rainfall and/or storm tide (Fig  189   5c-f). For the vast majority of the coastline, both the peak storm tide and 24-hour rainfall 190 are significantly correlated with TC intensity, although the strength of correlation is higher 191 for rainfall (Fig 5c-d). The 24-hour rainfall is also strongly negatively correlated with storm 192 translation speed (Fig 5f), as slower moving storms will drop more rainfall in a given 193 coastal location than faster moving storms. The peak storm tide is not strongly correlated 194 with translation speed (Fig 5e), since both slow and fast moving storms can generate high 195 surges, and the additional wind speed contribution from fast moving storms is generally 196 small compared to the cyclonic wind speed. Under future storm climatology, the 90 th 197 percentile of TC intensity is projected to increase by 15-30% along the Gulf of Mexico and 198 Southeast Atlantic, 30-50% along the Mid-Atlantic, and 20-30% along the New England 199 coastline (Fig 5a). The vast majority of previous studies also project an increase in North 200 Atlantic TC intensity, and many predict an increase in the frequency of high intensity 201 (category 3-5) TCs 26 . We also find a large future reduction in the translation speed of 202 storms that exceed the 90 th percentile intensity (Fig 5b). For all regions except New 203 England, storms that exceed 90 th percentile intensity are likely to move 20-30% slower in 204 the future than in the historical period. The decrease in translation speed found here is 205 consistent with previous work examining changes in translation speed in the historical 206 record 27 and projections of TC translation speed under future climate conditions 28-30 . The 207 increase in storm intensity coupled with the decrease in translation speed drives an 208 increased likelihood to observe both extreme rainfall and extreme storm tide in the future 209 and increases the upper tail dependence between the hazards. By comparing Figure 4b  210 with Figures 5a-b it is clear that most regions experiencing a significant increase in the 211 hazard dependence also experience significant increases in storm intensity and decreases 212 in translation speed. The Chesapeake Bay is a notable exception, since the hazards are 213 projected to become less dependent in the future even though there is an increase in TC 214 intensity and decrease in translation speed. In the future a larger number of intense storms 215 but they still induce extreme rainfall. Thus, the increase in the number of these types of 218 storms causes a decrease in the hazard correlation at this location in the future climate. 219 220

Discussion 221
The results presented here provide new insight into how the spatial pattern of TC hazards 222 and their co-occurrence may evolve in the future under a combination of SLR and changing 223 storm climatology. In particular, we project large increases in joint hazard across the US 224 East and Gulf coasts, with the most extreme increases (up to 195-fold decrease in JRP) 225 occurring along the Mid-Atlantic and New England coastlines. We also find that the impact 226 of storm climatology change on joint hazard is larger than the SLR impact for 96% of the 227 coastline. This is an important finding, since previous work that examined the influence of 228 storm climatology change on extreme sea level change found large impacts at low latitudes 229 but small impacts at high latitudes 31 . Here, we find that storm climatology change is a 230 dominant contributor to future joint surge-rainfall hazard at northern coastal locations, 231 mostly due to projected increases in rainfall hazard. Lastly, we find that the statistical 232 dependence between extreme rainfall and storm tide increases in the future for large 233 swaths of the coastline, resulting in a higher probability to observe multi-hazard extremes 234 during future storm events. This finding is significant since many previous studies of future 235 compound flooding have neglected potential increases in hazard dependence 9,11,12,32 , which 236 could underestimate compound flood risk. 237 The findings presented here are associated with inevitable uncertainties. We utilize 238 a single TC model to downscale all GCMs and reanalysis data, and the model predicts a 239 significant increase in future TC frequency for five of the eight GCMs. Although a few other 240 studies 33,34 have also predicted increases in TC frequency, the majority of studies predict a 241 decrease or no change in North Atlantic storm frequency 26 . However, the main findings of 242 our study are unchanged even if we assume no change in future TC frequency. The future 243 JRP change calculated by holding TC frequency constant at the historical level is slightly 244 lower at each coastal location (up to 149-fold decrease in JRP; see comparison in Fig S6), 245 but the spatial trends (i.e. higher JRP change in the north compared to the south) are 246 unchanged. The relative importance of TC climatology change compared to SLR also 9 a larger JRP change than SLR for 84% of the coastline. The reason our results are relatively 249 unchanged if we neglect frequency change projections is because the increase in TC 250 hazards and their joint occurrence is largely driven by projected increases in TC intensity 251 and decreases in translation speed. 252 Although the findings from this work cannot directly predict future compound flood 253 hazard, which must be quantified using high-resolution coastal inundation models, we 254 provide evidence that joint rainfall-surge extreme events could become an increasing 255 threat to coastal communities in the future. Our modeling results and characterization of 256 joint rainfall-surge probability distributions can be used to develop flood mapping 257 scenarios 35 for regional 9,36  which is an axisymmetric vortex model coupled to a 1D ocean model 45 . For each TC the 281 outer radius at which the cyclonic wind speed goes to zero (henceforth outer radius) is 282 randomly drawn from an empirical lognormal distribution 46 . We neglect the variation in 283 outer radius size over the TC lifetime 47 since previous work has shown the outer radius 284 variation to be relatively small 48 . Using the CHIPS-estimated intensity and outer radius, we 285 estimate the radius to maximum winds based on a theoretical wind model that links the 286 outer descending region of the TC with the inner ascending region 48 . We assume no change 287 in the distribution of TC outer size for the future climate since historical trend analysis for 288 the North Atlantic basin found no statistically significant changes in TC size over time 49 . 289 Moreover, an analysis of dynamically-downscaled TCs based on RCP 4.5 end of century 290 forcing found nearly constant outer radius compared to the historical period 50 . 291 292

Bias Correction 293
The downscaled TCs from each GCM may be biased compared to the NCEP-downscaled 294 TCs, and biases within the TC characteristics can propagate to become biases in the hazard 295 estimation. TC intensity and annual frequency are both important drivers of coastal flood 296 risk, and both variables are likely to be biased due to the GCM projections and the method 297 of downscaling. Therefore, we perform bias correction at the storm level based on the 298 difference between the NCEP TC intensity distribution and the GCM-predicted intensity 299 distribution for the historical period. We then bias correct the annual TC frequency 300 (independently from the intensity bias correction) at each location based on the NCEP-301 downscaled historical frequency and the GCM-downscaled historical frequency. Using our 302 method of bias correction, we avoid multivariate bias correction on the modeled storm 303 tides and rainfall, which often fails to preserve the entire dependence structure between 304 hazards 51 . Additionally, bias correction at the storm level is computationally efficient, while 305 bias correction at the hazard level requires performing intensive hydrodynamic 306 simulations for thousands of historical period GCM TCs. 307 2100) and historical  downscaled Vmax quantiles is added to the NCEP-311 downscaled historical quantiles to create a corrected future Vmax distribution for each 312 GCM model at each location. Then by the principle of importance sampling 53 the GCM-313 projected storms are weighted and re-sampled with weights corresponding to the ratio of 314 the corrected Vmax probability density to the GCM-projected Vmax probability density. By 315 doing weighted re-sampling of the storms at each location we are able to match the 316 corrected future Vmax distribution and consequently generate a storm set at each location 317 that is unbiased with respect to the intensity distribution. Figure S7 shows the bias 318 correction procedure applied at a sample location for a sample GCM, demonstrating that 319 after weighting/re-sampling the target Vmax distribution is matched accurately. After 320 correcting the Vmax distribution, we bias correct TC frequency by adding the GCM-321 predicted frequency change to the NCEP-derived frequency at each location. Nevertheless, a recent study modeled compound flooding using TCR-predicted rainfall 353 fields for several historical events and found that TCR rainfall produced similar flood 354 depth/extent compared to using radar rainfall forcing 22 . In our study, we utilize area-355 averaged TCR rainfall over each coastal catchment delineated according to USGS hydrologic 356 units (HUs) 58 , and each coastline point is paired with its upstream coastal catchment. We 357 utilize the maximum 24-hour rainfall accumulation (over the catchment) from each storm 358 event as our rainfall metric because the 24-hour storm duration is frequently used for 359 rainfall risk assessment studies 59 . 360 We treat TC climatology change and SLR as independent, although they may be 372 significantly correlated. Ref 60 found a significant correlation between SLR and changes in 373 power dissipation index (an integrated measure of TC intensity, frequency, and duration) 374 for the North Atlantic, suggesting that large increases in mean sea level are more likely to 375 co-occur with larger increases in TC hazard. By neglecting correlations between SLR and 376 climatology changes our results may underestimate the composite (weighted-average) 377 change in climatology and SLR, and consequently represent a conservative estimate of joint 378 hazard change. 379 380

Statistical Analysis of Joint Hazard 381
We conduct statistical analysis on the pairs of modeled storm tides (or storm tides plus 382 SLR) and 24-hr rainfall at each location along the coastline to quantify their marginal and 383 joint hazard. 384 The marginal distributions of both rainfall and storm tides are often characterized 385 by a long tail representing the rare but extreme events 41,42 . The heavy tail can be modeled 386 with a Peaks-Over-Threshold approach, where the probability above a high threshold is 387 estimated by a Generalized Pareto (GP) distribution 61 . We fit marginal GP distributions 388 using the maximum likelihood method 61 for the rainfall and storm tides at each location, 389 and the threshold is set by numerically minimizing the root mean square error between the 390 empirical quantiles and the theoretical quantiles. According to bivariate extreme value 391 theory, a logistic model can be used to estimate the joint distribution of two GP variables 392 such that their bivariate CDF takes the form 61,62 : 393 (1) 394 Where and are the Fréchet-transformed versions of the variables and , and a is a 395 parameter that quantifies the strength of the dependence between the variables (a à 0 396 signifies complete dependence and a=1 complete independence). At each location we fit 397 the bivariate distribution of extreme storm tide and rainfall, based on their respective GP 398 thresholds, using a censored maximum likelihood approach 62 (within the "evd" R-399 package 63 ). The bivariate logistic model employed here has previously been utilized to 400 model dependence between rainfall and storm surges 2,62,64,65 . 401 the historical period (Fig. S8). To additionally account for SLR impacts, we randomly draw 425 from the SLR probability distribution specified for each coastline location and add it to the 426 modeled peak storm tide for each event. We apply a standard bootstrapping approach with 427 500 iterations at each location and for each GCM storm set to quantify the sampling 428 uncertainty from the synthetic storm set and the SLR distribution. 429

Attribution of Changes in Joint Hazard 431
To quantify the isolated impact of various climate factors on changes in joint rainfall-surge 432 hazard, we adjust a single factor at a time and then re-calculate JRP. To quantify the 433 isolated impact of SLR on changes in JRP, we randomly draw SLR values from location-434 specific probability distributions 13 and add them to the historical rainfall-storm tide pairs. 435 The impact of changes in future storm frequency is quantified by simply changing the value The data generated from this study are deposited to the NSF DesignSafe-CI and can be 456 accessed online (https://www.designsafe-ci.org/}; link to the dataset will be provided 457 upon publication).