a, TG locations; colour represents the number of overlapping years between storm surge and precipitation time series, numbers in brackets denote the number of precipitation stations in a radius of 25 km around the TGs. b,c, τ for Cases I (b) and II (c); markers denote the selected copula type (square is Frank; triangles pointing right, left, up and down are Gumbel, Clayton, Galambos and Hüsler-Reiss, respectively). Insets show linear trends for the storm surge (red; unit is cm yr−1) and precipitation (blue; unit is mm yr−1) time series for Cases I (b) and II (c) for all sites used in the non-stationary analysis (see text). Open circles denote insignificant correlations or trends (90% confidence).
Data selection and processing.
Hourly sea-level records14 from 30 tide gauges (TGs) with at least 30 years of data between 1900 and 2012 are used for the study. If two or more TGs are located within 55 km, only the longer record is used, to avoid the same precipitation stations being allocated to more than one TG. The data at Mayport started in 1928 and ended in 2000, and was merged with Fernandina Beach (located only 31 km north of Mayport) TG data for the periods from 1897 to 1928 and 2001 to 2012. For the overlapping period (1985 to 2000) the correlation between the two raw data sets is 0.99. After this correction all records end in December 2012. The MATLAB T_Tide package31 is used for a year-by-year harmonic tidal analysis (with 67 constituents), which also effectively removes the MSL influence (including the seasonal MSL cycle). Years with less than 75% of data were omitted from the analysis (this affects less than 3% of the overall data). A visual inspection of the storm surge time series (tidal signal removed) revealed obvious outliers at several TGs. For Gulf and east coast stations these could be related to strong hurricanes or extratropical storms; on the west coast tsunamis were responsible for the highest water levels at several sites. Such events are driven by different mechanisms and often triggered thousands of kilometres away. Therefore we should not assume any relationship with rainfall events. We removed the following days from the records of all TGs located on the west coast: 23 May 1960 (Great Chilean Earthquake), 28 March 1964 (Great Alaskan Earthquake), 11/12 March 2011 (2011 Tōhoku earthquake in Japan). We note that measured water levels were possibly affected by tsunamis at other times throughout the past century, but did not always result in large sea-level anomalies at the sites of interest, and can therefore be neglected for the purpose of the present study.
Cumulative daily precipitation records15 from stations in a radius of 25 km around the TG sites are used. There are two exceptions to this rule: for Grand Isle we select stations in a radius of 50 km, as the TG is located offshore and there is only one precipitation station within 25 km, which covers a period much shorter than the TG record; for Neah Bay we use a radius of 35 km, as there are only few precipitation stations within 25 km of the TG and the time period covered by them is shorter than the TG record. All precipitation records allocated to a TG are averaged into a single time series and compared to the storm surge records. Generally, the precipitation time series are longer and more complete than the TG records; the entire precipitation time series are used to identify the overlapping periods which are then considered for the dependence study. A sensitivity test using a 50 km radius instead of 25 km reveals that the τ values derived for the 30 sites are not considerably affected by this definition (Supplementary Fig. 2). For Australia, significant correlation between precipitation and storm surge was found for much larger spatial distances (up to 1,000 km), but decreases further away from the coast12.
The significance of the linear trends of the raw data sets—that is, the annual maximum and coincident event time series (insets in Fig. 1b, c)—is tested with a modified Mann–Kendall Test for autocorrelated data32.
Dependence measure and modelling.
Dependence between storm surge and precipitation is measured through Kendall’s rank correlation coefficient τ (ref. 16) (as opposed to Pearson’s correlation coefficient, which captures only linear relationships). Storm surges sometimes occurred without precipitation, and this leads to ties (several zero values) in the Case I data sets. These ties are broken as suggested in ref. 23, by assigning ranks randomly. This may have a small effect on the rank correlation when the number of ties is large; therefore we repeat the procedure 100 times and calculate the average rank correlation. Copula parameters are derived with the maximum pseudo-likelihood estimator23. Copula models capable of simulating the prevailing dependence structures are selected by visually inspecting the observed and simulated (1,000 data pairs) pairs of ranks (or pseudo-observation clouds; see Supplementary Fig. 1 for examples) and using tail dependence analysis by comparing non-parametric TDCs (ref. 20) above a threshold of 0.6, derived from observed and simulated rank pairs.
Temporal changes in the overall rank correlation and tail dependency.
The rank correlation is calculated for 50-year moving windows (shifted by one year each time step). In each window we required at least 30 data pairs to assure robust estimates of τ. From this it follows, to be consistent, that at the beginning (end) of the time series, τ is calculated from the first (last) 30 data pairs if all values are available (that is, we shift the 50-year windows beyond the edges of the available time series and calculate τ values as long as 30 data pairs are available).
The significance of the temporal changes in τ is assessed in two ways. First, for six selected sites and cases (among Cases I and II) the range of the natural variability (10% and 90% levels) of the rank correlation is derived by resampling τ (from 50 data pairs) 10,000 times (Fig. 2a–f). Second, we calculate linear trends of the running τ time series for the common period since 1948 and test the significance with a similar resampling approach: the ranks of the raw data sets are randomly rearranged (1,000 times) and the running τ time series and linear trends are recalculated; the 30% and 70% quantiles are used to identify significant trends. Here we chose these relatively low levels because the number of degrees of freedom is small as a result of the running window approach, the linear model might not be ideal if τ increased rapidly recently after having been relatively stable before (for example, Fig. 2b), and most importantly the results from the multivariate statistical analysis reveal that the observed changes in τ—even if the trends are not significant at high confidence levels (for example, 90%)—greatly affect joint return periods; the estimated slopes at sites with significant trends range from ~0.002 to ~0.008, which equals an increase in τ of 0.1 to 0.4 over a 50-year period (or 0.2 to 0.8 over 100 years).
TDCs for the first and second halves of the raw data sets and above a threshold of 0.6 are compared to investigate whether similar changes as found in the overall rank correlation are also evident in the upper tails (Supplementary Fig. 3).
Tropical/extratropical and compound/non-compound events in NYC.
We use the HURDAT database24 to distinguish between tropical and extratropical storm surge events that were recorded at the Battery TG in NYC. When the time stamp of the day of a particular event occurs in the HURDAT database and the respective storm entered the region east of 70° W longitude and north of 30° N latitude at some point in its life we assume that it was a tropical event, otherwise it is flagged as an extratropical event.
We identify compound and non-compound events based on the empirical non-exceedance probabilities of storm surge (s) and precipitation (p). The events with grey squares in Figs. 3a and 3b fulfil the criteria s > 0.6, p < 0.4 and p > s + 0.4 (that is, s is high, p is small, and the events lie far above the diagonal), events with grey circles satisfy s > 0.6, p > 0.6, p < s + 0.15 and p > s − 0.15 (that is, s and p are high and events lie close to the diagonal).
Weather patterns associated with compound/non-compound events in NYC.
SLP, PWC and wind fields from the 20CR (ref. 26) and the days when the selected compound/non-compound events (see previous paragraph) occurred are averaged into composites, representing reference synoptic-scale weather situations favouring the occurrence of storm surges with much (Fig. 3c, d) or little (Fig. 3e, f) precipitation; the PWC composites are shown in Supplementary Fig. 5. We apply the concept of centred pattern correlation (CPC; ref. 27) as a widely used pattern similarity statistic in climate studies to correlate, for each spatial point, SLP and wind fields from all ‘event days’ in NYC with the composites used as reference situations (shown in Fig. 3c–f). We focus on the area between 70° to 40° W and 20° to 60° N, as the reference composites reveal that the prevailing SLP and wind situation in this area determines whether a storm surge is accompanied by large or small amounts of precipitation. The summed CPC coefficients (SLP + winds) with the low-precipitation reference situation (CPC2; Fig. 3e, f) are subtracted from the summed CPC coefficients with the high-precipitation reference situation (CPC1; Fig. 3c, d); if the difference is >0.5 the respective event is flagged as a compound event, if it is <−0.5 the event is flagged as a non-compound event, and if it is in between it is not allocated to either of the two. We then calculate ratios between compound and non-compound events from the Case I and II data sets for the entire records and for the first and last 30-year periods separately.
Multivariate statistical analysis for NYC.
We perform a full multivariate statistical analysis for the data set from NYC to demonstrate the effect of a changing rank correlation (that is, changes in the joint distribution) on multivariate return periods of variables relevant to flood risk analysis and infrastructure design/adaptation. The marginal distributions for Cases I and II were selected from a pool of five distributions commonly used in hydrology and (coastal) engineering, namely the LogNormal distribution, normal distribution, exponential distribution, Weibull distribution, and the generalized extreme value (GEV) distribution. A goodness-of-fit (GoF) test comparing the empirical and theoretical non-exceedance probabilities using the root mean squared error was applied to identify the distributions fitting best to the underlying data sets; distribution parameters were estimated with the maximum likelihood method. We use two Gumbel copulas—identified with the GoF test(s) employed here—to analyse the Case I and Case II data sets separately. To quantify the effect of changes in the dependency, the analysis is performed with the highest and lowest values of the running τ time series in Fig. 2e. The relevant ‘AND’ joint return periods29 are calculated for a selected compound event.