A general pattern of trade-offs between ecosystem resistance and resilience to tropical cyclones.

Tropical cyclones drive coastal ecosystem dynamics, and their frequency, intensity, and spatial distribution are predicted to shift with climate change. Patterns of resistance and resilience were synthesized for 4138 ecosystem time series from n = 26 storms occurring between 1985 and 2018 in the Northern Hemisphere to predict how coastal ecosystems will respond to future disturbance regimes. Data were grouped by ecosystems (fresh water, salt water, terrestrial, and wetland) and response categories (biogeochemistry, hydrography, mobile biota, sedentary fauna, and vascular plants). We observed a repeated pattern of trade-offs between resistance and resilience across analyses. These patterns are likely the outcomes of evolutionary adaptation, they conform to disturbance theories, and they indicate that consistent rules may govern ecosystem susceptibility to tropical cyclones.

INTRODUCTION

Tropical cyclones, including hurricanes and typhoons, are among the most powerful natural phenomena on Earth. Even weak storms can bring devastating rainfall and storm surges that cause catastrophic loss of life, damage property, and disrupt ecosystem services (). Predicting the socioecological effects of tropical cyclones is critically important as human coastal populations rise (, ), the spatial distribution of storms extends into higher latitudes (), and the frequency and intensity of storms increase (–). However, variation among storms, the diversity of affected ecosystems and responses, and the necessarily opportunistic nature of most hurricane research (, ) have generated a plethora of seemingly unique case studies. Synthesis is needed to reveal common predictor and response variables that are translatable across storms, systems, and processes to foster prediction of ecosystem susceptibility to future storms ().

Several disturbance frameworks describe how system state (e.g., abiotic conditions and biota life history) and disturbance mechanism (e.g., identity and intensity) drive the ecosystem responses (, ). These frameworks predict that naturally dynamic and frequently disturbed ecosystems exhibit higher intrinsic resistance and/or resilience, whereas variation in intrinsic resistance and resilience among individual species is expected to be a function of traits including generation time and mobility (). Intrinsic resistance and resilience are dictated by different processes and need not be positively correlated with one another (). For example, short generation time in a species is hypothesized to enhance intrinsic resilience through rapid population growth, but this trait does not necessarily enhance resistance (). Concordantly, traits that confer resistance need not also enhance resilience. Although we may intuitively expect intrinsic resistance and resilience to be positively correlated with one another, empirical evidence from multiple case studies suggests that the opposite may be true (–). This would imply that managers seeking to simultaneously enhance both resistance and resilience in coastal ecosystems face an impossible task (, ), but there is a need to test these ideas with a large-scale data synthesis. While new data collection through coordinated global research networks is the ideal path forward (), geospatial analyses and existing data allow us to immediately address this need ().

Here, we provide a comprehensive synthesis of coastal ecosystem susceptibility to tropical cyclones rooted in disturbance theory. We describe coastal ecosystem susceptibility to tropical cyclone disturbance as a combination of intrinsic resistance, the degree to which an ecosystem can remain unchanged despite the presence of disturbance, and intrinsic resilience, the ability of an ecosystem to return to the reference state after a temporary disturbance (). While these intrinsic ecosystem properties are unmeasurable latent variables (), we can measure their indicators: observed resistance, the unit increase in effect size per unit increase in the stressor (e.g., maximum wind speed or total rainfall), and observed resilience, the rate of change per unit time for ecosystem return to baseline conditions after disturbance (see Materials and Methods for formulas). Scaling our observed measures of resistance by stressor intensity and resilience by effect size differentiates them from the classic engineering definitions of resistance and resilience () and allows us to compare them across ecosystems and storms.

We assembled a dataset of 4138 ecosystem time series from 26 different storms primarily in the Atlantic basin affecting 118 study sites (Fig. 1) between 1985 and 2018. Using a disaggregation approach, we took each tropical cyclone event and divided it into multiple, translatable components including storm stressor attributes such as wind speed and total rainfall, ecosystem type such as coastal wetlands and open water, response variable categories including biogeochemical and organismal responses, and aspects of the response such as observed resistance and resilience (). We report repeated patterns in ecosystem responses to hurricanes, most notably a pattern of negative covariation, or trade-offs, between resistance and resilience among ecosystems and response variables. These consistent patterns across ecosystems and measured variables from such a large dataset can be treated as confirmed generalizations that support an empirically validated general theory of tropical cyclone disturbance in coastal ecosystems.

Fig. 1.

Map of sites with storm tracks and ecosystem types studied.

Middle panel is all site (dots) and storm track (lines) except data from Taiwan. Bordering center panels are cumulative rainfall values (bottom, right) and max wind speed (left, top) by longitude and latitude.

RESULTS

Data distribution

We assembled monitoring data from before and after 26 tropical cyclones across 118 different study sites from 1985 to 2018. Using time series (n = 4138) classified into five variable categories—biogeochemical, hydrographic, mobile biota, sedentary fauna, and vascular plants—we quantified the resistance and resilience to maximum estimated wind speed, ranging from 58 to 253 km/hour. and maximum estimated rainfall, ranging from near zero to 28 cm/day. (Fig. 1). Time series measurements included 151 unique response variables that were grouped into five variable categories: biogeochemical, hydrographic, mobile biota, sedentary fauna, and vascular plants. The greatest proportion of time series data fell into the biogeochemical (62%, n = 2561 series) and hydrographic (29%, n = 1217 series) categories, with the remaining data distributed among the mobile biota, sedentary fauna, and vascular flora. There were eight different ecosystems that were grouped into four ecosystem categories: fresh water (lotic and lentic; 6%, n = 233 series), salt water (estuaries and offshore; 88%, n = 3643 series), wetland (mangroves, salt marsh, and freshwater wetlands; 5%, n = 190 series), and terrestrial (forests and dunes; 2%, n = 72 series). The wide range of represented ecosystems, response variables, and stressor intensity metrics allows us to characterize variation in observed resistance and resilience among ecosystems and response variables and quantify cross-ecosystem patterns in relationships among these variables.

Ecosystems, organisms, and processes all exhibit a resistance-resilience trade-off

We observed a negative relationship between observed resistance and resilience associated to both wind speed and rainfall among time series across the entire dataset (Fig. 2 and tables S1 and S2), suggesting that coastal ecosystem functions, flora, and fauna all exhibit trade-offs in response to tropical cyclones. Simulation models demonstrated that the observed relationships are unlikely to have occurred by chance (see the Supplementary Materials). Observed resilience was more tightly correlated with observed resistance to wind speed than with observed resistance to rainfall. There was wide variation in the observed resistance to wind speed, observed resistance to rainfall, and observed resilience, and these responses differed among ecosystems and variable categories (Fig. 3 and table S3). However, we again observed patterns of negative covariation between resistance and resilience. Terrestrial ecosystems had higher biogeochemical resistance than all other ecosystems to both wind and rain, but the lowest resilience (Fig. 3 and table S4). In contrast, freshwater ecosystems had the highest biogeochemical and hydrographic resilience and lowest biogeochemical resistance to wind and rain among systems. Mobile biota in freshwater systems also had the highest resistance to rainfall among systems, followed by wetland and then saline ecosystems. Among variable categories within ecosystems, mobile biota and biogeochemical responses tended to have the lower resistance to wind and rain but the higher mean resilience. In contrast, hydrographic variables tended to have higher resistance and lower resilience. Other variables showed more system-specific patterns. Vascular plants had low resilience and high resistance to wind within wetlands, but not in saline systems. Sedentary fauna also showed no difference from other variables in saline systems, likely due to small sample size and high variation within that category.

Fig. 2.

Relationships between resistance to wind or rain and resilience among ecosystem and variable types.

Ecosystem types, ecosystem processes, and flora and fauna all demonstrate negative covariation between resilience and (A) resistance to wind speed and (B) resistance to rainfall. Metrics are on a natural log scale and should be considered unitless for simple interpretation and viewed in relative terms from low to high. Each graph displays the stability summary metrics calculated from each time series labeled by ecosystem (color) and variable category (shape). Best fit lines shown include the relationships for each group-variable category combination in the varying slope and intercept model with line color indicating ecosystem type.

Fig. 3.

Comparison of resistance and resilience among variable categories within ecosystem types.

Columns correspond to ecosystem type, rows correspond to response metric, and colors denote variable category. Violin plots show the distribution of data within each group, while inset box plots show the median, quartiles, and outliers (points). Bars over each plot indicate which variables are not significantly different from one another within ecosystems. Alphanumeric labels under each violin plot indicate which variable categories were significantly different among ecosystems (columns). Letter denotes variable category (e.g., b = biogeochemistry, h = hydrography, m = mobile fauna, s = sedentary fauna, v = vascular plants). Number denotes rank, with different numbers across plots (e.g., b1 versus b2) indicating a significant difference. Red text indicates that for the response variable (row) and the indicated variable category (letter), we observed significant differences among ecosystem types (table S4).

We suspect that variation in the evolutionary history of species’ adaptation to disturbance exhibits scale-up behavior and regulates the susceptibility of populations, communities, and the ecosystem processes they drive. For individual species, there is selection to excel at either resistance or resilience, but rarely both simultaneously (, ). This has been observed in individual case studies, but not as a widespread pattern among field observations spanning storms, systems, and organism types. For example, tree species in Jamaican forests following Hurricane Gilbert showed strong among-species covariation between damage exhibited (low resistance) and fast regrowth rates (high resilience) (). Dominant trees tended more toward resistance, and those that exploit light gaps (disturbance specialists) tended more toward resilience (). Certain sets of life history characteristics such as long life, delayed sexual maturity, or low reproduction, and organism traits such as low adult mobility dictate a need for greater resistance to the disturbance events that are likely to be experienced within the organism’s lifetime, as we observed in sedentary and vascular plants.

In contrast, the relatively low resistance but high resilience of mobile biota in our synthesis supports the idea that there is less advantage for resistance to cyclone effects in highly mobile species that have generation times substantially shorter than the local return interval of tropical cyclones. These species can rapidly recolonize and proliferate in the wake of disturbances that are infrequent relative to their life span, and strong selection would require substantial lasting changes to the gene pool as a consequence of the event (). Mobile biota can be locally extirpated or emigrate in advance of the storm and then rapidly return or recolonize the affected habitats (). A similar argument can be made for single-celled organisms with short generations such as phytoplankton and microbes (), although they may be more likely to undergo rapid evolution in response to the event and/or via general mechanisms of dispersal (). Areas of coastline or particular ecosystem types such as streams experiencing more frequent disturbances may be dominated by species that have evolved to be successful within that regime (). Intermediate cases such as fast-growing weedy plants (low adult mobility, but high fecundity and short life span) or elasmobranchs (high adult mobility with low fecundity and delayed reproduction) may be phylogenetically constrained to evolve to maximize one strategy or the other ().

The concept of evolutionary trade-offs may also explain negative covariation between observed resistance and resilience in the physical and chemical variables that are controlled by biota. For example, ecosystem efflux of carbon dioxide is driven by heterotrophic microbial activity (), erosion rates are constrained by the stability of the plant roots stabilizing the soil (), nutrient dynamics are influenced by microbial uptake rates (), and watershed nutrient export is constrained by the uptake rates of vegetation on the landscape (). Given the biological control over these variables, the observed resistance of a physical or biogeochemical process should be at least partially explained by the observed resistance of the dominant biota responsible for the process.

Furthermore, the observed resilience of a process dependent on biota is reliant on both the degree of functional redundancy among species supporting the process () and the variation in intrinsic resilience to disturbance among those species (, ). Compensatory dynamics could result in high observed resilience of a process if secondary taxa with high observed resistance and/or resilience are present and able to fill the missing role in the absence of the dominant taxa. However, in cases where the entire guild of taxa responsible for a process has low intrinsic resilience, the dependent process will exhibit low observed resilience. For example, following Hurricane Hugo (1989) in Puerto Rico, stream nitrate concentrations remained elevated for years, likely due to the naturally long duration required for downed vegetation to completely break down and for relatively slow-growing tree communities to recover (). In principle, the tighter the link between the measured ecosystem function and the biota, and the greater the similarity of intrinsic resilience among taxa driving the function, the more closely the function should adhere to the pattern of a resistance-resilience trade-off. Together, the average observed resistance and resilience of entire ecosystems should reflect the evolutionary forces on the resident biota, the dependence of requisite processes on the biota, and the degree of redundancy in the biological communities ().

Responses that are completely or largely independent of control by biota should tend toward lower resistance, because the lack of biological regulation allows an efficient energy transfer between the force applied and the physical response. For example, a drop in estuarine water conductivity would be simply a function of the amount of storm-related freshwater inflow into the basin relative to the volume of the receiving estuary (). The resilience of these variables must then be positively related to the potential energy of the new state. For example, estuaries vary in their water residence time as a function of mean daily watershed inflow, tidal range, and basin morphology (). Estuaries with higher residence time will thus take longer to return to baseline water conductance after massive storm-related freshwater inflow. This pattern was observed in estuaries affected by Hurricane Harvey, in which those with high residence times returned to baseline within 7 months, while those estuaries with low residence time returned within a few weeks ().

Storm characteristics affect ecosystem resistance and resilience

We observed overall positive relationships between measures of resistance and the total rainfall, maximum rainfall, and maximum wind speed, the slope of which varied among response categories within ecosystems (Fig. 4 and tables S1 and S2) with clear splits between aquatic (saline and fresh water) and nonaquatic (wetland and terrestrial) habitats. In the aquatic systems, biogeochemical responses were enhanced with wind speed and reduced with rainfall. A similar pattern occurred with hydrographic variables, although freshwater responses were reduced to both wind and rain (Fig. 4 and table S2). Terrestrial systems showed an opposing pattern, with biogeochemical responses enhanced with rainfall and reduced with wind speed (Fig. 4 and table S2). The amplified sensitivity of biogeochemical and hydrographic variables in aquatic systems to rainfall likely reflects the dependence of these responses, largely measured in the water column and inclusive of microbial and phytoplankton responses, to freshwater inflow and loads of particulate and dissolved materials carried into the ecosystem (). In contrast, rainfall is less likely to have large effects on biogeochemical responses, primarily litter and leaf chemistry in this dataset, in terrestrial systems. Conversely, direct effects of wind on wave energy or water column stability are less likely to create lasting changes in aquatic systems, such as upwelling-induced phytoplankton blooms, to the water column constituents or microbial and algal communities beyond a few weeks (), whereas wind-induced damage to terrestrial vegetation can cause lasting impacts on uptake and biogeochemical cycles on the landscape ().

Fig. 4.

Relationship between storm metrics (total rainfall, max rainfall, and max wind speed) and resistance among ecosystems within each variable category.

Points display the conditional random effect (x axis) of each storm metric (shape) on each variable category (row) in each ecosystem (color) from the respective mixed effects models for each storm metric. Error bars are the conditional standard deviations of the random effect.

Vascular plant response was enhanced in all cases in aquatic ecosystems and reduced in all cases of terrestrial and wetland ecosystems (Table S2). This divergence is most likely explained by the buffering effects of water to rainfall and high winds on submersed aquatic vegetation. While high winds can be locally damaging, severe impacts on these organisms is more strongly linked to sedimentation () and prolonged reduced water clarity following a storm (, ). These effects are more strongly dependent on factors like basin morphology and residence time and overall storm impacts as the watershed scale than the specific characteristics of the storm event directly over the submersed plants ().

In the saline ecosystems, mobile biota, primarily fish, showed an increase in resistance as rainfall increased and a decrease in resistance with wind speed (Fig. 4 and table S2). The positive relationship between observed resistance of mobile biota and rainfall intensity could be a function of drops in barometric pressure that may serve as a warning for mobile biota to prepare for a storm by emigrating from the region or seeking shelter (). For example, Strickland et al. () observed that juvenile bull sharks (Carcharhinus leucas), driven by barometric pressure cues, evacuated a South Florida estuary 24 hours before Hurricane Irma. Common snook (Centropomus undecimalis) in Florida rivers also exhibited similar behavior, moving downriver into deeper waters as Hurricane Irma approached (). Freshwater systems did not show this same pattern, possibly because there is no similar downstream deep-water refuge habitat for obligate freshwater fish in the shallow coastal rivers that were included in the dataset. This explanation may also apply to mobile biota in wetlands, primarily fiddler crabs and snails, and sedentary fauna in saline systems, primarily infaunal bivalves and invertebrates, which all have limited ability to seek refuge against freshwater inundation and salinity changes because of high rainfall. The different spatial patterns of wind and rain within a storm offer one possible explanation for why increasing wind velocity does not show a similar effect on mobile biota in saline systems (, ). Specifically, the highest winds are relatively localized, and thus, there is high spatial variability in wind speed across the entire affected region. Topographic variation can create additional variation in local wind speed. While rainfall also has some spatial variability, watersheds remove some of that variability by integrating rainfall across the entire basin before the storm flows affect downstream ecosystems.

In contrast to observed resistance, observed resilience demonstrated no relationship with stressor (i.e., wind speed and rainfall) intensity (table S1). The lack of a relationship suggests that resilience, or recovery rate, within the range and types of effects included in these analyses is independent of stressor magnitude. This result indicates that while a larger perturbation to an ecosystem may require a longer recovery period than a smaller perturbation, the actual rate of recovery would be similar in both cases. We interpret this to mean that the limits on recovery rates are primarily dictated by intrinsic resilience and ecosystem conditions that are independent of perturbation magnitude such as the available food or nutrient resources during the recovery. However, stressor-invariant resilience may not always be the case. For example, a very severe disturbance that completely extirpates a taxon from a region would result in a slower rate of recovery (e.g., low resilience) if there are time lags until the species recolonizes from another source population, if there is displacement by a new opportunistic species, or if the habitat has changed ().

DISCUSSION

Gaps and directions in global hurricane research and coastal ecosystem management

Our data assembly efforts were largely focused on North America. A broader geographic scope, including both the Southern Hemisphere and the Pacific basin, would enhance future efforts. We also observed biases in the types of data gathered in the synthesis. Estuarine and offshore ecosystems were exceptionally well represented, and coastal terrestrial ecosystems had the poorest representation. Although our dataset was not exhaustive, the relative abundances of data types indicate that a general research bias toward coastal marine ecosystems may be pervasive and, thus, limit the scope of our interpretations. There were also biases in the types of response variables in our dataset, with poorest representation in animals, plants, and ecosystem processes such as decomposition rates, plant growth, or nutrient cycling. Because of data limitations, we were also unable to differentiate microbial populations (including phytoplankton) as distinct from the biogeochemical processes they control. Our study highlights the need for more studies of the effects of tropical storms on other parts of the coastal plain (e.g., streams, forests, and freshwater marshes) and more research on diverse organisms spanning the domains of life and ecosystem processes.

Beyond gathering additional data, our synthesis raises several important research questions that are critical for understanding the factors that drive ecosystem susceptibility to tropical cyclone disturbances.

What sets of species traits best explain patterns of resistance and resilience at the population-level?

For instance, generation time, reproductive potential, time until maturity, dispersal mode, dispersal distance, and physiological tolerance limits at different life stages could be quantified and placed into an analysis. A trait-based approach to responses of biota could yield advances in our mechanistic understanding and our ability to predict population-level effects of a stormier future.

How do past disturbance regimes of a region influence the response of ecosystems to subsequent disturbances?

In other words, are coastal ecosystems that have a long history of frequent hurricane effects more resistant to future storms, or does frequent disturbance generate ecosystems with low resistance but high resilience to disturbance? Are these relationships linear, or are there tipping points? Disturbance operates in both a metacommunity context, filtering out species lacking the requisite traits for long-term survival, and in an evolutionary context, driving populations in high disturbance regions to exhibit population-level differences from congeners in less frequently disturbed regions. Our synthesis data could be further classified in terms of disturbance regime (frequency and magnitude of prior events) to explore patterns. Field experiments using genetic tools and common garden experiments could further test these ideas in a more mechanistic fashion (). These questions are particularly important for forecasting the effects of the increased frequency of hurricanes in higher latitudes predicted with climate change.

How do antecedent environmental conditions immediately before a tropical cyclone influence its impact?

Ecosystems are inherently dynamic, changing throughout the year and between years in accordance with seasons, demographic cycles, weather conditions, and climatic oscillations (e.g., El Niño Southern Oscillation, North Pacific Oscillation). Do these factors interact with storm characteristics to dictate the impact of a tropical cyclone? For example, is a storm more damaging in a drought year or an exceptionally wet year? Is the impact of early season storms comparable to late season storms? To illustrate, Hurricane Agnes was an exceptional rainfall event that wiped out the majority of submersed aquatic vegetation in the Chesapeake Bay in 1972 but is thought to have been particularly disruptive because it arrived in June before seed set had occurred for the majority of the aquatic plant species (). Information about antecedent conditions and timing could be incorporated into future analyses to quantitatively evaluate the influence of recent weather and timing on ecosystem sensitivity.

How applicable are these concepts to understanding the stability of coupled human-natural coastal ecosystems in the face of tropical cyclones?

The resistance and resilience of natural ecosystems must depend, in part, on the local-scale effects of humans on those ecosystems. Similarly, the sensitivity of the natural ecosystems may influence the sensitivity of both human infrastructures by moderating effects such as flooding and storm surges and human economic activity dependent on ecotourism and natural resource harvesting. Incorporating metrics of anthropogenic system sensitivity to tropical storms and quantifying the effect of anthropogenic activities on the resistance and resilience of coupled natural ecosystems should enhance our ability to understand these complex system dynamics. Furthermore, these investigations will enhance our ability to understand the relative balance of local-scale anthropogenic stressors and global stressors (e.g., climate change) on resistance and resilience among ecosystems, geographic localities, and response variables and help target effects to enhance coastal resilience or resistance ().

Although much remains to be learned, our synthesized data show consistent patterns of trade-offs in ecosystem sensitivity to tropical cyclones among coastal ecosystems, which may extend more broadly to applied ecology such as management responses to wildfires and ecological disturbance theory. Our findings emphasize that managing for increased strength in one axis (e.g., resistance or resilience) may result in decreases in the other, an importance consideration for coastal management. For example, in the Chesapeake Bay, seagrass meadows dominated by fast-growing widgeon grass (Ruppia maritima) have higher resilience but lower resistance to water quality declines than meadows dominated by eelgrass (Zostera marina) (). This trade-off is an important consideration that must be evaluated when making coastal management decisions such as choosing which species of seagrass to restore. In areas with infrequent disturbance, resistance may be better to manage for, whereas resilience may be preferable in areas experiencing frequent perturbations.

These results enhance our cross-ecosystem understanding, highlight new questions, and suggest that an empirically validated general theory of controls on ecosystem susceptibility to tropical storms is achievable. Our synthesis only represents a fraction of published and unpublished data on ecosystem responses to tropical cyclones, and there remains much to learn from historical data. Continuing to build on this approach by learning from the past will enhance our ability to predict the effects of tropical cyclones on coastal ecosystems in the future.

MATERIALS AND METHODS

Meta-analysis

A multidisciplinary group of participants in the Ecosystem Responses to Hurricanes Synthesis Meeting held in Corpus Christi, Texas, in April of 2019 assembled the dataset. The meeting included 52 participants representing a diverse range of expertise, geographic regions of the United States and Taiwan, demographics, and career stages, who all had time series data on ecosystem responses to hurricanes to contribute to the synthesis. Participants submitted their own data, identified and contacted collaborators and colleagues with additional data, and helped identify publicly available data sources that could be included in the dataset.

Each submitted datum was required to come from a time series of observations of a single response variable (e.g., dissolved oxygen, species abundance, and salinity) taken repeatedly in a fixed location before and after an individual storm event. The frequency of observation (e.g., 15-min, daily, and monthly intervals) had to be high enough to reasonably capture the temporal pattern, as determined by best professional judgment of the investigator who collected and submitted the data. A single, decades-long time series might be broken into multiple data points, each observation representing a different storm event that occurred during the overall series. Likewise, a given storm event typically had multiple data points representing different spatial locations where data were collected and different response variables were measured at each location. For each datum, we recorded the name of the response variable, definition of the response variable, units of the response variable, latitude, longitude, storm name, date of the event, mean prestorm value, value of the greatest change caused by the storm, and time in days until observations returned to typical values for that site at that time of year. In addition to covariate and site information, we also created meta-data to describe the source of the data and the funders to aid future researchers accessing this dataset. We converted measurements of the same response variable to a common set of units among entries.

To understand how storm characteristics affected different coastal ecosystems, processes, and their components, we grouped all data by ecosystem type and response variable type. Ecosystem types included wetlands, fresh water, salt water, and terrestrial. We defined wetlands as systems dominated by hydrophilic plants living in hydric soils and included plants influenced by fresh water and salt water such as grass-dominated marshes and mangroves. We defined freshwater ecosystems as coastal rivers and streams. We defined saltwater ecosystems as estuarine and marine ecosystems. We defined terrestrial ecosystems to include coastal forests as well as several minor habitats such as coastal dune vegetation. We classified response variables into five categories: biogeochemical, hydrographic, vascular plants, sedentary fauna, and mobile biota. We defined biogeochemical variables as stocks and fluxes of dissolved constituents such as nutrients and trace elements, as well as photosynthetic prokaryotic and eukaryotic microbes that influence the cycling of those materials. The decision to combine geochemical stocks and fluxes (i.e., nutrients) with biotic drivers of those processes (i.e., phytoplankton biomass and microbial populations) was, in part, driven by an overrepresentation of bulk phytoplankton biomass (as chlorophyll-a) in the assembled dataset and limited data on abundances of distinct species or groups of microbes. We recognize that this approach limits our ability to resolve the complex interrelationships between these factors and believe these limitations represent important research gaps for future microbially focused efforts. We defined hydrographic variables as measurable properties of waters such as salinity and turbidity that are influenced by both physical and biogeochemical drivers. We divided measures of the abundance, biomass, and diversity among the mobile biota, animals able to actively disperse during their adult life stage, sedentary fauna, animals with limited or no mobility during adult life stages, and vascular plants.

Storm characteristics

To describe how tropical cyclone characteristics affect coastal resilience and resistance, we chose two key aspects of storm disturbance known to affect population, community, and ecosystem dynamics: wind and rainfall.

Maximum wind speed

As most weather stations are damaged before hurricane force winds are registered, we relied on the relationship between wind speed and distance from the storm center to estimate wind speeds at all locations. We built this relationship from surface wind field data from the U.S. National Hurricane Center (www.nhc.noaa.gov/gis), which provide the maximum extent of sustained wind speeds at thresholds of 35, 50, and 64 knots. We used all available wind fields to model the relationship between wind speed and distance from the storm’s center, with different models constructed for all categories of the Saffir-Simpson cyclone wind scale. Then, using the IBTrACS dataset of all tropical cyclone tracks (), we used the distance of each storm to a study site location to estimate the maximum sustained wind speeds at each of the roughly 1-hour interval locations. Also, from the IBTrACS dataset, we constructed separate models relating a storm’s sustained winds to maximum wind gusts, and this relationship was used to estimate maximum wind gusts at each location. Analyses were performed in the statistical program R () using sf() library ().

Maximum and total rainfall

To ensure a uniform coverage, precision, and accuracy of rainfall data for each storm and location, we estimated rainfall metrics from multiple remote sensing–derived datasets. Datasets included the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) dataset and Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks-Climate Data Record (PERSIANN-CDR). From each rainfall dataset and unique location/storm combination, we extracted rainfall data for a 4-day period temporally centered on the time when the storm’s center track came closest to the location. At each location, we created a 1000-m radial buffer (2000-m diameter) around the center point, and grid cell values for average precipitation, maximum precipitation, and temperature from the various datasets (e.g., PERSIANN and CHIRPS) for the identified dates were extracted and averaged using the location buffer. We defined the hurricane event as 2 days before to 2 days after the storm “event” date. We then calculated rainfall within the 24-hour period of the storm passing, the highest daily estimate of rainfall, and total rainfall, the sum of cumulative rainfall over the time period. Estimates of each storm metric from each data source were then averaged for each site/storm combination. Analyses were performed using Google Earth Engine ().

Analyses

For each time series of ecosystem responses to tropical cyclones, we calculated metrics indicative of observed resistance and resilience that were modified from Pimm () to normalize for the variation in effect size and stressor magnitude as well as to create intuitive axes (). We calculated observed resistance as the ability for an ecosystem response variable to withstand change, as the change in effect size per unit increases the stressor (e.g., wind speed). Effect size is normalized maximum amplitude of the response, calculated as the natural log response ratio (LRR) of the absolute observed change relative to the baseline valuewhere is the average prestorm value of the response variable, and r are poststorm observations taken up until the return to baseline. We calculated Observed Resistance by taking the natural log of the ratio of LRR to the measured value of the stressor (i.e., wind speed and total rainfall) and multiplying by negative one to flip the axis for ease of interpretabilitywhere str is the stressor. We performed this calculation for each measurement of wind speed (maximum gust) and rainfall (total rainfall, rainfall within 24-hour period) and averaged resistance values across measurements within stressor classes.

To estimate Observed Resilience, the rate at which ecosystem response variables returned to baseline (i.e., rate of return), we calculated the natural log of the ratio of percent absolute change of each ecosystem response variable divided by the time to return to baselinewhere d is the number of days until return to baseline.

Given that the formulas for observed resistance and resilience have a common numerator and one is multiplied by negative one, these metrics are biased to be negatively related to one another. However, before statistical quantification of the relationship between observed resistance and resilience, we evaluated this bias using simulated data and determined that relationships in our data are unlikely to have been produced by chance and can be quantified and interpreted using standard statistical approaches (Supplementary Materials).

To quantify the relationship between resistance and resilience, for each resistance metric (one per storm metric), we fit an exponential decay linear mixed effects model using maximum likelihood with the natural log of resistance as our response, resilience as a fixed effect, and variable category within ecosystem type as a random slope and intercept effect to account for differences between variable types and ecosystems.

To understand variation in metrics of ecosystem susceptibility among systems and variable response categories, we ran a one-way analysis of variance (ANOVA) with type III error distributions to account for unequal sample sizes followed by a Tukey–post hoc test for each combination of stability metrics (n = 3) among variable categories for each system type (n = 4) and among systems for each variable type present in multiple systems (n = 4). ANOVA models were evaluated with a Dunn-Šidák adjusted α value of 0.0064 to control for family-wise error rate. This sequential approach was used instead of two-way ANOVA because of combinations of variable types. To evaluate the effect of stressor intensity on resistance and resilience, we ran a series of mixed effects models (n = 6) with each sensitivity metric (i.e., resistance to wind, resistance to rain, and resilience) treated as the response, and a paired storm intensity metric treated as a fixed effect, with variable category within ecosystem type as a random intercept effect to account for variation in the way groups among ecosystems respond to stressor intensity. Interpretations were restricted to coefficients fit to specific variable/system combinations. We fit three models with resilience as the response, one for each of the storm metrics (i.e., wind speed, total rain, and maximum rain). Models with resistance as the response variable included one where resistance to wind was paired with maximum wind speed, and two with resistance to rain as the response and either total rain or maximum rain as the predictor. We split the models in this fashion to aid in interpretability of results.

We performed all analyses using the statistical program R version 4.0.3 (). We fit mixed effects models using the lme4 package () and reported results from conditional F tests with Satterthwaite degrees of freedom using lmerTest package for fixed effects (). To assess model assumptions, we evaluated randomized quantile residuals using the DHARMa library and found no violations (). We calculated variance explained using the r.squaredGLMM function in the MuMIn library and reported as marginal R2, variance explained by fixed effects only, and conditional R2 values, variance explained by whole models ().