Recent increases in tropical cyclone precipitation extremes over the US east coast

All seven of the adjusted latewood (Materials and Methods) chronologies were significantly positively correlated with regional instrumental ΣTCP but, importantly, did not respond strongly to regional non-ΣTCP (SI Appendix, Table S1). These results support the hypothesis that late-season wood production (i.e., latewood width) responds primarily to the large, short-duration inputs of precipitation produced by TCs (27, 28). Due to different record lengths of the chronologies, the overall reconstruction model is composed of a series of nested principal-component reconstruction models (Materials and Methods) that explain between 34 and 41% of the variance in regional instrumental ΣTCP (SI Appendix, Fig. S1 and Table S2). After combining the nested models, our reconstruction spans 1700 through 2015 CE (SI Appendix, Table S2). The reconstruction overestimated low (1- to 40-mm) ΣTCP values and underestimated high (>100-mm) values, and we used a quantile-mapping bias-correction procedure (29) (Materials and Methods) to reduce these systematic differences (Fig. 2).

Using quantile regression to estimate linear trends for the entire bias-corrected reconstruction (Materials and Methods), the upper tail (≥0.95 quantile) of ΣTCP shows a significant increase over time, whereas other quantiles exhibit either negative or nonsignificant trends (Fig. 3). This increasing trend starts at the 0.90 quantile (∼1 mm/decade) with the 0.95 and greater quantiles (years with ΣTCP of at least 174 mm) showing stronger trends (∼2 to 4mm/decade) (Fig. 3B). A total 6 of the 10 yrs with ΣTCP in the 0.95 quantile happened after the year 1900. Further, 5 of the 7 yrs of highest ΣTCP occurred in the late 20th/early 21st centuries (1955, 1996, 1999, 2016, and 2018). For years of extreme (≥0.95 quantile) ΣTCP in the 21st century, the estimated linear trend of 2 to 4 mm/decade translates into an increase of 64 to 128 mm compared to years of extreme ΣTCP in the early 1700s.

The trend in the reconstruction toward increasingly more extreme ΣTCP is supported by the instrumental record, which also shows its largest ΣTCP extremes (>300 mm) in more recent years (Fig. 2A). The two most recent years of extreme ΣTCP include 2016, with Tropical Storm Julia and Hurricane Matthew producing a combined 334 mm of TCP, and 2018, when Hurricane Florence produced 292 mm of TCP, causing $10.3 and $24.2 billion in damages, respectively (30, 31). We also examined whether the seasonal timing of TC occurrence differentially impacted latewood growth, potentially influencing the reliability of the reconstruction but found no significant differences in adjusted latewood growth between different months of TC occurrence (SI Appendix, Fig. S2).

Only in 2 yrs in the early to mid-18th century, 1703 and 1741, did extreme ΣTCP rival these recent extremes (Fig. 3A). These earlier extreme ΣTCP are unrelated to increased variance in the nested tree-ring models (SI Appendix, Fig. S3). Further, historical documents substantiate that a strong TC affected a large area of the East Coast from Virginia to Pennsylvania in October 1703 (32). Our results suggest that TCP from this storm also affected the North Carolina and South Carolina coasts. Our reconstruction further suggests an undocumented storm (or multiple storms) in 1741 that affected the North Carolina and South Carolina coasts, a year currently without historical documentation of heavy rainfall or flooding from TCs, to the best of our knowledge.

ΣTCP in a given year can be influenced by the annual number of TCs that make landfall, the average duration of these TCs, and their translation speed. Therefore, an increase in ΣTCP extremes can indicate an increase in one or more of these factors. For our region, we found that instrumental (1948 to 2018) seasonal ΣTCP was significantly positively correlated with average (r = 0.807, P < 0.01) and summed (r = 0.893, P < 0.01) seasonal TC duration but negatively correlated with seasonal average TC translation speed (ρ = −0.602; P < 0.01; SI Appendix, Table S3). Correlations with these TC metrics were of the same sign and significant, albeit weaker, for our reconstructed ΣTCP over the instrumental period (SI Appendix, Table S3) and with the longer Hurricane Database (HURDAT2) dataset time period (1851 to 2018) (SI Appendix, Table S3; correlation with TC translation speed is not significant).

The relationship between ΣTCP with TC duration and TC translation speed, with a particular emphasis on extremes in both metrics, is further demonstrated by the results of our SEA analyses (Fig. 4 D and E). Years with extreme (≥0.95 quantile) average TC durations or TC translation speeds (HURDAT2; 1851 to 2018) are associated with significantly larger ΣTCP (Fig. 4D). Similarly, years with extreme (≥0.95 quantile) reconstructed ΣTCP are associated with significantly longer average TC durations and slower TC translation speeds (Fig. 4E). Our findings thus suggest that the TC slowdown in recent decades, expressed by lower TC translation speeds and longer durations, is an important control on the trend of increasing extreme ΣTCP.

While our findings support a relationship between ΣTCP and annually averaged TC duration and speed, we also find a relationship between ΣTCP and the number of TCs (SI Appendix, Table S3), making it difficult to determine whether the increases in ΣTCP extremes are due to increases in the number of TCs or slower TCs or both. However, we find that the number of TCs making landfall in our study region has not changed using Mann–Kendall trend test (−0.002 TCs per year; τ = −0.067; P = 0.8821) since the start of the HURDAT2 (33) dataset in 1851, whereas their translation speed has decreased (−0.026 kn per year; τ = −0.137; P < 0.05) and their average duration has increased significantly (0.01 h per year; τ = 0.093; P < 0.05) (Fig. 4 A–C). To determine the role that the number of storms in a given year could have on our results, we examined storm-level statistics (Materials and Methods). We found a significant relationship (ρ = −0.480; P < 0.01; 1948 through 2018) between the translation speed per TC and the TCP per TC (Materials and Methods). We examined trends in TCP from TCs that have similar translation speeds (binned in groups of 5 kn) to further demonstrate the role of TC speeds on the increasing trend of ΣTCP extremes. In the three groups (5 to 9 kn, 10 to 14 kn, and 15 to 19 kn) that had a sufficiently large sample size (≥∼20 TCs), TCP did not have a significant trend through time (SI Appendix, Fig. S4). However, if we examine TCP per TC from every storm from 1948 to 2018, there is a weak increasing trend (τ = 0.111; P = 0.08; SI Appendix, Fig. S4). These results support that the TCP of a given storm is related to the speed of that TC and those speeds have been slowing through time (Fig. 4 A and B), resulting in larger extremes in ΣTCP.

While instrumental ΣTCP are weakly related to climatic drivers along the US Atlantic coast, particularly the Bermuda High and El Niño Southern Oscillation (ENSO) (10, 34), reconstructed ΣTCP have minimal to no connection to these climate drivers (SI Appendix, Table S3). Sea surface temperatures (SSTs) in the tropical and North Atlantic have increased in recent decades and have been linked to increased TC activity (35). However, ΣTCP were not significantly correlated with either gridded SSTs for any portion of the North Atlantic (1850 through 2016; Hadley Center) or with reconstructed western Atlantic tropical SSTs (36) (SI Appendix, Table S4), suggesting that it is unlikely that climatic/oceanic interactions are driving the increase in ΣTCP extremes.

To generate seasonal ΣTCP, we used daily instrumental gridded (0.25° × 0.25°) TCP data from TCPDat (10) (1948 through 2018), which is a publicly available product (https://github.com/jbregy/TCPDat). TCPDat combines gridded precipitation data from the Climate Prediction Center (48) with the TC best track data from HURDAT2 (33, 49), available from the World Meteorological Organization’s International Best Track Archive for Climate Stewardship initiative version 03r09 (49). TCPDat provides robust spatial coverage of TCP allowing the observed record to not be dominated by a single weather station. While the hurricane season lasts officially from June 1 to November 30, most longleaf pine trees stop growing by October 15 (50). We therefore defined seasonal ΣTCP for each year as the sum of daily TCP for the period June 1 to October 15, the period best captured by latewood ring width (27). We created a regional (spatial) average of ΣTCP based on 43 contiguous grid points that received at least an average of 40 mm/year and ensured that every grid point that contained a tree-ring site was included (Fig. 1). We also calculated seasonal non-ΣTCP by subtracting the ΣTCP from seasonal total (June 1 to October 15) precipitation for each of the 43 grid points and then spatially averaging for each year.

We used five previously published latewood-width chronologies (27, 28) along with two newly sampled longleaf pine latewood width chronologies for a total of seven sites across North Carolina and South Carolina (Fig. 1). We further extended two of the previously published chronologies with remnant wood cross-section data (SI Appendix, Table S5). All seven sites are located on Carolina bay sand rims and support longleaf/wiregrass woodland communities. Seasonal ΣTCP in these sand rims, which rise 0.5 to 1.5 m above the surrounding area and drain water efficiently, account for approximately half of the variance in latewood growth from resident longleaf pine (27, 28). At each site, we targeted every canopy-dominant tree with old-growth morphology (13 to 36 trees per site). For each tree, we extracted two samples from opposite sides at ∼1.3 m above the ground. To extend the chronology, we also extracted cross-sections from stumps and remnant wood using a chainsaw.

We measured the earlywood, latewood, and total ring width using WinDendro (version 2017a) for all core samples and statistically cross-dated the samples using the program COFECHA (51). Similarly, we measured two transects for each cross-section. We detrended individual tree-ring time series to remove the biological growth trend and nonclimatic influences of growth (e.g., disturbance) by using a two-thirds spline (52) with the adjustment for endpoints (53) in the program ARSTAN (54). Fire is relatively infrequent at these sites compared to other pine savannas, but it is still part of the ecosystem. Thus, we used the two-thirds spline as the detrending method to maximize the retained climate signal and remove nonclimatic noise. Lastly, we created chronologies for each site using a biweight average of the individual series (i.e., core and cross-section measurements; SI Appendix, Table S5).

Climate and other environmental conditions in the early growing season influence latewood width in other Pinus species via an inertia effect (55, 56). When removing the influence of climate in the early growing season in the southwestern United States, the resulting adjusted latewood width has a stronger relationship to monsoon precipitation compared to the unadjusted latewood (55, 56). So far, only unadjusted latewood widths have been used to reconstruct ΣTCP (27, 28). To remove the influence of early growing season climate on latewood width, we regressed earlywood width onto latewood width and retained the residuals for each individual measurement series (55, 56). By adding one to the zero-mean residuals, we produced an adjusted latewood index (LWa) that, like standard ring-width measurements, has a mean of one over the entirety of the time series (55). The sensitivity of P. palustris latewood growing on Carolina bay sand rims is related to microtopographic differences (57). Therefore, the sensitivity of individual trees to ΣTCP varies. By calculating LWa at the series level, we can potentially remove some of the influence of landscape position on latewood and better retain the ΣTCP signal (SI Appendix, Fig. S5).

We used a nested principal component (PC) regression (58) to reconstruct regional ΣTCP from our seven LWa chronologies. We first conducted a PC analysis of the seven LWa site chronologies. We retained only the first PC as a predictor, as it was the only PC with an eigenvalue greater than one and consistently explained over ∼80% of the variance in the LWa chronologies (SI Appendix, Table S2). To account for decreasing sample depth back in time, we used a nested reconstruction approach (59). The first reconstruction nest covers the common period of all chronologies (1880 through 2015). We then dropped the youngest chronology as determined by the expressed population signal (60) (SI Appendix, Fig. S3) and reran the PC regression with the remaining chronologies as predictors and instrumental ΣTCP as the predictand. We repeated this process until there was only one chronology remaining. The final reconstruction is the combination of the nested models, with each nested model being calibrated from 1948 to 2015.

To evaluate the reconstruction skill of each nested reconstruction model, we used a standard split-sample calibration and validation approach. We divided the instrumental period into two periods of approximately equal length (1948 through 1978 and 1979 through 2015) to cross-validate the ΣTCP reconstruction model. Model validation was performed over the later period (1979 through 2015), and verification statistics were calculated for the early period. The verification was based on the coefficient of determination (R²), the Reduction of Error (RE) (61), and the Coefficient of Efficiency (CE) (62). We then reversed the calibration period to the latter half and validated on the early half using the same metrics. Positive RE and CE values indicate adequate skill of the model at predicting ΣTCP through time. We did this validation procedure for each reconstruction nest. The earliest nest (based only on the Lewis Ocean Bay chronology) had low explained variance, and we therefore excluded that nest from the reconstruction. To assess model uncertainty, we show the inner 95% (for each year) of 999 reconstruction replicates derived from a maximum entropy bootstrapping (MEBoot) (63). MEBoot preserves the dependence structure (i.e., the autocorrelation and partial autocorrelation functions) of the time series and does not assume stationarity of the time series.

To examine how the timing of the TCs may have influenced LWa, we examined the month of occurrence for each TC that generated TCP from 1948 to 2015. To better attribute ΣTCP to a specific month, we removed years that did not have TCP and removed years that had TCs from multiple months contributing to the ΣTCP. Using this subset of data, we used an ANOVA with a Tukey honestly significant difference post hoc test to compare the LWa grouped by the month of TC occurrence and found that there were no significant differences between groups (SI Appendix, Fig. S2). We also performed an ANOVA on TCP for each month of occurrence and found no significant differences (SI Appendix, Fig. S2). It is important to note that the sample size of these groups was relatively small (n ranged from 4 to 13); however, August and September account for over 80% of TC activity in the north Atlantic Basin, which is the peak season for latewood growth. Further, the timing of TC occurrence did not exhibit a trend through time, and every TC that produced >100 mm occurred in August and September, indicating that the timing of TCs had a minimal impact on the interpretation of the reconstruction.

Because the probability distributions of tree-ring reconstructions often differ from those of instrumental data, especially in the tails, we bias-corrected (29) our ΣTCP reconstruction. We bias-corrected each nested model separately, then combined the bias-corrected values for each nested model to get the final bias-corrected reconstruction of ΣTCP. Bias-correction is used here as a postprocessing procedure, similar to its use with climate model outputs (64). We used quantile mapping with the “RQUANT” method in the qmap package (65) in R (66). RQUANT approximates the quantile–quantile relationship by using local linear regressions on regularly spaced quantiles. The bias of a given quantile in the reconstruction is estimated by calculating the difference between the local regression and the quantile derived from the observations. These different bias estimates for each quantile can then be added to the respective quantiles from the reconstruction. Another benefit of bias-correcting the reconstruction is that we can directly compare the reconstructed ΣTCP values with very recent observations, such as those during 2016, 2017, and 2018, when we do not have well-replicated tree-ring data. Once systematic error within the reconstruction has been reduced, a direct comparison to observed values is more appropriate. Therefore, we include the three most recent years from the instrumental data in the bias-corrected reconstruction of ΣTCP. To determine how the postprocessing bias-correction procedure influenced the reconstruction, we calculated the Nash–Sutcliffe CE to assess overall model error (1948 through 2015) before and after bias-correction. We also examined the ratio of the total error that was systematic or related to bias by calculating the systematic mean square error and dividing it by the total mean square error (29, 67). Kernel density estimates of the probability density functions demonstrate that the bias-corrected reconstruction better matches the empirical distribution of the observed data compared to the uncorrected reconstruction (Fig. 2D). The bias-corrected reconstruction also addresses the problem of the reconstruction model producing negative ΣTCP estimates. To assess the bias-corrected reconstruction model uncertainty, we applied the same bias-correction procedure to each of the 999 members of the ensemble from the uncorrected reconstruction derived using MEBoot. We then estimated the percentiles for the 95% CIs from the bias-corrected ensemble. While the fit of each nested model is not included in the calculation of the CIs, the RMSEs for each nest only differ by ∼2 mm. Further, for the majority of the reconstruction (from 1740 onward), the RMSE values for the various nested reconstructions only differ by ∼1 mm. Overall model evaluation statistics are slightly improved by the bias-correction procedure, with calibration CE = 0.255 before versus 0.279 after bias-correction for the most replicated nest. The proportion of total error that was systematic error slightly improved after bias-correction, with 20% of the error being systematic in the original reconstruction and 18% following bias-correction.

We used the HURDAT2 data set (1851 through 2018) to determine the number, translation speed, and duration of landfalling TCs in our study area. We tabulated any TC that tracked within 223 km—the average size of the moisture plume from a landfalling TC and the diameter that was used in TCPDat (10)—of the center of our study region to generate a time series of the number of landfalling TCs. We calculated translation speed (14) but only for the portions of the TC tracks that fell within 223 km of our study region. We used the distances between the 6-h locations along each TC track using a great circle arc to calculate translation speed. We then averaged the translation speed for all segments of all TC tracks within 223 km of our study region for each year to generate a time series of seasonal average translation speed of landfalling TCs. To calculate seasonal average TC duration, we summed the number of 6-h periods from HURDAT2 during which a TC was located within 223 km of our study region and then divided it by the number of TCs that were in our study region for a given year. Because translation speed as well as the number of TCs influence ΣTCP, we also summed the number of 6-h periods from HURDAT2 that a TC was located within 223 km of our study region for each year to create a time series of seasonal total TC duration that combines the number of TC and the speed of those TCs. We then use seasonal average translation speed, seasonal average TC duration, and total TC duration to examine the relationship of TC speed and duration with instrumental and reconstructed ΣTCP using Spearman’s correlation.

To determine potential relationships between extremes in seasonal average TC durations and ΣTCP, we conducted a series of superposed epoch analyses (SEA). In a first SEA, we used reconstructed extreme (≥95th quantile; n = 12) ΣTCP years as event years and seasonal average TC duration as a chronology for the current year and over the previous year and following 5-y lags. We then conducted a second SEA in which we used extremes (≥95th quantile; n = 12) in seasonal average TC duration as event years and reconstructed ΣTCP as the chronology. We then repeated the same procedure (creating two more SEAs) but replaced seasonal average TC duration with the seasonal average TC translation speed. For the SEA using transition speed as the event series, we had to use the 90th quantile to reach a number of events high enough for the analysis (n = 11). All SEAs were conducted for the common period (1851 through 2015). We used 1,000 Monte Carlo simulations to develop bootstrapped CIs to determine statistical significance.

To examine the relationship between TCP and TC speed for individual storms, we computed the TCP per TC and the translation speed per TC. We calculated the TCP for each storm by first summing the daily values of TCP for those days that a storm impacted our study region (from TCPDat) for each of the 43 contiguous grid points used in the ΣTCP reconstruction (Fig. 1). We then averaged the TCP data across the 43 grid points to obtain a regional TCP estimate for a given storm. We calculated translation speed per TC for the portions of the TC tracks that fell within 223 km of our study region, but we averaged the speeds for each storm instead of averaging for the entire season as was done for the seasonal average TC translation speed. We then used Spearman’s correlation to calculate the correlation between TCP and translation speed per TC. Lastly, we binned individual TCs in groups that had similar speeds (5-kn groups). For each group that had a sample size of ∼20 or more TCs (5 to 9 kn, 10 to 14 kn, and 15 to 19 kn), we tested for trends in TCP per TC with a Mann–Kendall trend test with Sen’s slope using the “MannKendall” function in the R package “Kendall” (68). We then tested for trends in TCP per TC for all individual TCs that influenced our study region to determine whether trends in TCP were related to the translation speeds of TCs.

To quantify trends in extremes of reconstructed ΣTCP over time, we use quantile regression (69), using R package “quantreg” (70), to estimate the slope of a regression using time as the predictor and ΣTCP as the response. Whereas ordinary least squares (OLS) regression provides an estimate of the conditional mean (i.e., an estimate of the mean ΣTCP for a given year), quantile regression provides trend estimates for ΣTCP across the full probability distribution. Quantile regression is formulated as an optimization problem that estimates the conditional quantiles by minimizing the sum of absolute errors that are weighted according to their quantile position (i.e., it does not partition the data into quantiles and then estimate parameters). Here, we estimate the linear relationship between quantiles of ΣTCP and time for quantiles ranging from 0.01 to 0.99. Given that we are particularly interested in changes in extreme ΣTCP, we focus on trends in the upper tail of the distribution, such as the 0.95 quantile and higher, and compare them to trends in the center of the distribution (Fig. 3). Like OLS, quantile regression provides estimates of the uncertainty of different quantiles. In our application (using a zero-limited variable), the higher quantiles clearly have higher uncertainties compared to the lower quantiles (Fig. 3B). We used the “rank” method to determine the 95% CI around each quantile regression slope. While our bias-correction procedure did not inflate the variance during the instrumental period (variance of 5,162.3 in the original reconstruction compared to 5,271.7 in the bias-corrected reconstruction), bias-correction increased the variance during the preinstrumental period, with the variance being 43% higher over the full reconstruction in the bias-corrected time series. The estimates of uncertainty that are derived from the quantile regression trends (Fig. 3) became larger after bias-correction. As a result, the changes in uncertainty that come from the distribution-altering properties of bias-correction are, to a large extent, captured in the quantile regression-based trend estimates. To determine whether the annual number, the seasonal average duration, and the seasonal average translation speed of TCs (1851 through 2018) showed a significant trend, we used a Mann–Kendall trend test with Sen’s slope using the “MannKendall” function in the R package “Kendall” (68).

To determine the role of large-scale atmospheric-oceanic modes of variability on ΣTCP variability, we compiled climate indices with a well-known influence on TC genesis, development, and tracks. Specifically, we used the bivariate ENSO index (BEST) (71) time series (June through November; 1950 through 2016) from the National Oceanic and Atmospheric Administration Physical Science Laboratory. We also used the summer North Atlantic Oscillation index (NAO; 1948 through 2018) (72) and June through November averages of the Atlantic Multidecadal Oscillation (AMO; 1856 through 2018) (73) and the Bermuda High Index (1948 through 2018) (74, 75). Lastly, we gathered gridded SSTs for the North Atlantic from the Hadley Center (76) for the period 1850 through 2016. We then calculated Pearson correlation coefficients between all climate drivers and the ΣTCP reconstruction over their period of overlap. To compare over the full reconstruction period, we also compared our ΣTCP reconstruction with AMO (May through September) (77) (annual) (78), NAO (winter) (79, 80), NAO (summer) (81), ENSO (annual) (82) (winter) (83), and tropical North Atlantic SSTs (April through March tropical year) (36) reconstructions.