The response of ecological systems to environmental disturbances has been a focus of ecology for decades (Holling, 1973).
Understanding the consequences of disturbance underpins the theoretical frameworks of successional dynamics (Prach & Walker, 2011), interspecific coexistence (Chesson, 2000), and consumer-resource interactions (Silliman et al., 2013). Today, there is particular interest in the 'resilience' of systems to disturbances-primarily disturbances that have a negative impact-as the frequency and magnitude of extreme environmental pulse events increases with climate change (Jentsch et al., 2007). Specifically, natural resource managers could benefit from understanding how systems will respond to disturbance, how long recovery will take, and what management actions could improve resilience.
Resilience is a nebulous concept in ecology, with many definitions applying to different aspects of a system's response to a perturbation (Capdevila et al., 2020). Here we focus on the definition Holling (1996) termed 'engineering resilience': the magnitude of the initial impact of a pulsed disturbance on a system, and the time it takes for the system to recover to its pre-disturbance state after the disturbance ceases. Pulse disturbances are discrete events that abruptly change the state of a system (e.g. abundance or biomass), or ecological parameters such as mortality or fecundity (Jentsch & White, 2019;Yang et al., 2008), with the event stopping after some time. Examples of pulse disturbances in different systems include windstorms or wildfires (e.g. Buma & Wessman, 2011), atmospheric rivers flushing estuaries (Cheng et al., 2016) and mortality from red tides and marine heatwaves (Laurel & Rogers, 2020;Summerson & Peterson, 1990).
A useful way to characterize pulse responses of different systems in comparable terms is to separate the response into (a) the disturbance impact, the initial decline during and immediately following the disturbance event (this is often referred to as 'resistance'; Capdevila et al., 2020), and (b) the subsequent return trajectory, or recovery time, to the original state (Figure 1, Ingrisch & Bahn, 2018).
This framework presumes that the original state of the system was a stable attractor, that the pulse disturbance did not move the system into the attracting basin for a different stable state, and that environmental conditions do not preclude a return to the original state (though that framework could be adjusted to accommodate some of those scenarios; Yeung & Richardson, 2018). This bivariate framework is useful for conceptualizing differences in engineering resilience among disparate systems, but it does not provide a quantitative, dynamic explanation of why different systems would exhibit different responses, or what the shape and timing of the trajectory would be. Here we address this question of what factors determine the population dynamic response to a pulsed disturbance for the specific case of age-structured populations of a single species.
While much of the interest in resilience centers on higher levels of biological organization such as communities or ecosystems (e.g. Hillebrand & Kunze, 2020), those systems are ultimately comprised of multiple interacting populations and species. Furthermore, conservation and management actions are often focused on individual populations, such as in fisheries or endangered species. Moreover, age-structured population models capture the fundamental time scales associated with transient dynamics following disturbances (Capdevila et al., 2020), as set by mortality and maturity schedules, in a way that simpler unstructured models used to investigate resilience (e.g. logistic) cannot (Botsford et al., 2019;Hastings et al., 2018). As such, examining population dynamic responses to pulse disturbances can produce direct insights for conservation and management.
Much of our existing understanding of population responses to disturbance comes from the analysis of models of linear, agestructured populations (e.g. Leslie matrix models). Here we extend this thinking to nonlinear age-structured models, as others have recommended (Capdevila et al., 2020;Ezard et al., 2010;Stott et al., 2011). This is particularly important because linear models are more appropriate for describing the dynamics of populations at low density, where density-dependent vital rates would be unimportant, but those models will ultimately grow without bound as recovery proceeds-making them useful only over short-time horizons. By contrast, nonlinear models can capture the dynamics of populations near their non-zero steady states. The quantitative framework proposed by Ingrisch and Bahn (2018) implicitly assumes a system with a steady state, so our analysis of the factors leading to interspecific differences in resilience using nonlinear models fills an important gap. Further, we develop our analysis with a renewal equation approach, and we show how that provides an analytical method for characterizing the shape of the trajectory of the population response to different types of pulsed disturbances.
increasing or decreasing levels of disturbance) would have different effects on population resilience: delayed responses for species in which biomass is concentrated in older age classes, and for disturbances that become progressively more severe.
age-structured population model, marine heatwave, pulse disturbance, recovery time, resilience, resistance The recovery of a population from a pulse disturbance comes either from reproduction in the local population, or immigration from other locations. We focus here on the case where adult movement is relatively restricted, so recovery depends on new juveniles entering the affected population. This applies to many types of species with a juvenile stage with high dispersal potential, such as winddispersed seeds, aquatic insects, or the planktonic larvae of marine fish and invertebrates. The new juveniles can be produced locally by the reproductive adults in the disturbed population, or they can arrive from other locations. Often the relative contribution of local vs. external reproduction is unknown, so we focused on the two extreme cases marking the end-points of a continuum: (1) closed populations (i.e. no dispersal) and ( 2) open populations with recruitment occuring at a constant level from external sources. We consider both cases because the degree of demographic openness can shift which demographic rates have the greatest influence on dynamics (Yau et al., 2014).
We developed an analytical model of population responses to pulse disturbances, which allows us both to understand general principles and to predict the responses of populations whose life history characteristics have been quantified. We structure our analysis as an investigation of demographic features that we expect to affect both the impact and recovery time following disturbance: the degree of demographic openness, the range of age classes affected by the disturbance, and whether the population is harvested. Additionally, we consider the case of longer 'pulses' that could vary in intensity over time and bridge the gap into 'press' disturbances (Capdevila et al., 2020). These analyses provide a broad understanding of the dynamics of this type of resilience in age-structured populations with density-dependence.
A renewal equation is a useful way to represent age structured population dynamics, particularly when a major source of variability in the population is fluctuations in the production or survival of offspring each year (Botsford et al., 2019). Because the recovery from disturbance depends on new offspring ('recruits'), this representation focuses on recruits, although the dynamics produced are identical to that of models based on Leslie matrices (Botsford et al., 2019). Instead of keeping track of the vector of abundance at each age a at time t, n a,t , a renewal equation simply represents the current abundance of a population at time t in terms of past recruitment to the population (i.e. the entry of individuals into the first age class, a = 0) at time t -a for each age-a cohort, multiplied by the fraction of that recruit cohort that survives to age a (up to maximum age A). Thus, if recruitment for an age-a cohort is R t-a and survivorship to age a is a , the expression for total abundance at time t, N t , illustrates this equivalence: Recruitment, R t , is calculated as the total number of eggs (or live offspring) produced in time step t, E t , multiplied by the probability of surviving the juvenile period, t , before entering the first age class of the censused adult population (e.g. this would correspond to survival during the pelagic larval period for marine fishes or the time to fledging in birds):
(1)
F I G U R E 1 Schematic illustrating the components of response of population biomass to a pulsed disturbance (grey bar indicating several years of lower juvenile survival). The two variables used to quantify resilience, impact and recovery time (the latter defined as the time to return to 95% of the original state) are illustrated.
We modelled dynamics for both a demographically open population and a demographically closed population. Equations ( 1) and ( 2) describe the dynamics of both, but in the former, E t represents offspring arriving from some outside source (presumed to be constant), and in the latter, E t depends solely on reproduction within the population (which would change if the population is reduced by a disturbance). Presumably most real populations fall in between those two extremes, but quantifying connectivity precisely is challenging, particularly for species with very small propagules. Therefore we use these two extreme cases as bounds for potential model outcomes (as in White et al., 2013;Yau et al., 2014).
In the open population model, an undisturbed population with constant recruitment R comes to a stable equilibrium (i.e. a steady state to which the system would return if perturbed). The closed model, however, would not come to a stable equilibrium unless some process is density-dependent. Therefore, we assumed that juveniles experience Beverton-Holt-style within-cohort density-dependent survival prior to entering the census population:
where is the survival at low density and sets the maximum density.
To keep results similar across different simulated species, we parameterized such that the population would persist deterministically so long as E t was >20% of the value in an unperturbed population. Note that both the open and closed population models are nonlinear with stable, non-zero equilibria.
We initially considered disturbances that affected only the survival of juveniles, so effectively reductions in t . This type of disturbance is observed, for example, in the responses of some marine fish and seabird populations to marine heat wave conditions (e.g. Laurel & Rogers, 2020;Piatt et al., 2020), the response of salmon populations to drought conditions that reduce survival of outmigrating smolts (Notch et al., 2020), or the response of prairie bird populations that experience reduced fledgling survival during drought (Yackel Adams et al., 2006).
We investigated 16 species of demersal nearshore fishes, most of them in genus Sebastes (rockfishes; Table S1). These species represent a wide range of life histories, in terms of natural mortality rates (and thus average lifespan), growth rates and ages at first maturity (Table S1). In addition, as a case study, we investigated the dynamics of Pacific cod Gadus macrocephalus in the Gulf of Alaska.
The age-structured model described by Equation (1) describes abundance as the state variable, so the right-hand side of Equation (1) contains two terms: the number of recruits at time t - a and the age-dependent function a giving the proportion of those recruits that will survive to age a. If we were interested in a different state variable, such as a biomass, then a different agedependent function would be needed. Additionally, reproduction is an age-dependent process. These various age-dependent functions are important determinants of the dynamics, and we will show in the results that the disturbance response differs if one quantifies it in terms of changes in biomass rather than abundance. Therefore we describe these functions, also termed influence functions (Botsford et al., 2019) here.
The first age-dependent relationship is the probability of survival to age a, a . This depends on the age-specific instantaneous natural mortality rate, M a,t , and in analyses that include harvest, the instantaneous harvest rate, F. It is calculated as the cumulative probability of surviving from each age to the next, up to age a ( 0 is by definition equal to 1 because that is the recruit age class):
where a gives the proportion of the incoming age class a that is exposed to harvest. Note that the mortality rate M a,t could also vary with time and age, such as if a pulse disturbance increased the mortality rate of certain age classes for a period of time. In those cases we define M a,t = M a 1 + t such that t is the time-varying proportional increase in mortality.
We also calculate biomass-at-age, assuming that biomass is related to length. Length-at-age is described by a von Bertalanffy function (applicable to many species with indeterminant growth) with asymptotic maximum length L ∞ , growth rate k, and age-atlength-0 equal to 0. Biomass-at-age is then a function of length-at-age, with allometric parameters u and : b a = uL a v . Finally, fecundity, expressed here as egg production (for the demographically closed model) is approximately proportional to biomass (as in many species in which clutch size increases with maternal size):
where a is the proportion of age-class a that is reproductively mature, and is the number of eggs produced per unit biomass.
The relationships implied by Equations ( 4)-( 6) of abundance, biomass, fecundity, and fishery yield as functions of age are depicted in Figure 2, using black rockfish Sebastes melanops as an example, with different values of the natural mortality rate, M, to illustrate the sensitivity of the shape of those curves to that life history trait.
It is convenient to refer to these relationships as influence functions because they quantify the relative influence of different age
classes on different demographic processes (Botsford et al., 2019).
Conveniently, many possible changes to model assumptions, such as a mortality or fecundity rate that eventually declines with age, would simply change the shape of the influence functions but not alter our general conclusions.
In the case of an open population where E t is a specified constant, Equation (1) has the useful interpretation of being a weighted moving average of recruitment, R t , in which the weights are the values of the survival-to-age function, a . Mathematically, this is the convolution of the two functions R t and a . Note from Equation (1) that in convolution, the ages of the weightings increase in the direction opposite to time; in other words, the recruitment time series and the weighting function are moving in opposite directions in the time domain (see Supplemental Information for an illustration of this calculation). The benefit of this insight is that the solution to the open population model can be calculated from the convolution of the recruitment function R t (as modified by disturbance) and any of the age-dependent functions in Figure 2, without direct simulation.
Moreover, this calculation also applies when the pulse disturbance itself has different patterns over time; for example, a gradual reduction in juvenile survival ( t ) over multiple years as conditions worsen, rather than a discrete period of anomalous low survival.
We began with an analytical examination of the model to determine what factors determined the impact and recovery time following a pulsed disturbance (as in Figure 1), as well as shaped the population's recovery trajectory. That analysis revealed in addition to the intensity and duration of the disturbance, the most important factor was the natural mortality rate of the species (see Section 3). Based on those results we then simulated the response of populations to pulse disturbances ranging from 10% to 70% reductions in early life survival (parameter t ) for 5-20 years (values chosen for illustrative purposes). We made those calculations for the 16 California species, to illustrate how the response depends predictably on life history variables. We also compared responses in abundance to those in biomass to illustrate how the different patterns of those metrics as a function of age (Figure 2) affect the population response, and we compared open and closed population dynamics. We initially considered only disturbances affecting the survival of offspring, and thus the recruitment of juveniles into the population but also extend that analysis to disturbances affecting all age classes. We also considered the effects of varying levels of harvest, F.
Finally, because actual populations likely experience a mixture of disturbance effects to both the early life stages and adult life stages, we considered as a case study the type of disturbance experienced by Pacific cod during the 2014-2016 marine heat wave (Barbeaux et al., 2021;Di Lorenzo & Mantua, 2016). During that event of prolonged, anomalously warm, low-productivity water, the adult cod mortality rate increased by an estimated 68% (Barbeaux S1), with yield given for the value of the harvest rate associated with depletion to 33% of the unfished biomass (the most recent estimate for that stock). Distributions are shown for the actual estimated natural mortality rate for blue rockfish (black curve; M = 0.18 year -1 ) as well as values of M 50% less (blue curves) and 50% greater (red curves) to illustrate the effect of that parameter. et al., 2021), and the reduction in spawning habitat with a suitable temperature range reduced survival of early life stages of cod to 58% of its pre-heatwave level (Laurel & Rogers, 2020). We did not attempt to recreate the exact population trajectory of the cod population, which depends idiosyncratically on the history of harvesting, socio-economic factors, and stochastic recruitment events preceding the disturbance (and which was already characterized by Barbeaux et al., 2021). Rather, we examined the relative contributions of those realistic levels of juvenile and adult disturbance to the dynamics of a population that was at a stable equilibrium and age distribution prior to the disturbance. We made separate simulations with both disturbances and each one individually to compare their effects. We assumed the population had harvest rate F such that the population's reproductive output (E t ) was approximately 26% of its unfished value (as estimated for that stock in 2014; Barbeaux et al., 2021), and in this analysis we examine the trajectory of both abundance and fishery yield following disturbance. We assumed that the population was demographically closed and used life history parameters drawn from the Barbeaux et al. (2021) stock assessment.
All analysis performed in Matlab R2020b (9.9.0.1524771).
Code available at https://doi.org/10.5281/zenodo.7199549 (White et al., 2022). This study did not require ethical approval.
We began by examining how total abundance, N t , would respond to a disturbance that reduces the survival of new offspring, t by some proportion for a duration years. Prior to the pulse, equilibrium abundance, N * , will be at a constant level determined by the value of constant recruitment R t = R, specified for the open population model (when there is no harvest):
Once the disturbance begins, abundance of recruits will decrease by factor 1δ each year (so the population will have a geometrically decreasing trajectory), and the decrease associated with each annual cohort will propagate through the age structure of the population as determined by a . If we define the time immediately prior to the onset of the disturbance as t = 0 (so that t = 1 is the first year in which recruitment is reduced), then the abundance at the end of the final year of the disturbance, t = , will be The first term on the right-hand side of this equation contains the age classes that recruited prior to the pulse disturbance, and the second term on the right-hand side are those age classes that occurred after, hence were affected. This solution clarifies the three factors affecting the initial impact on population abundance. Not surprisingly, the impact is greater if the disturbance or duration are greater. But the value of natural mortality also plays a role: if survival to age a declines faster with age (i.e. if natural mortality M is greater) then the young vulnerable age classes make up a greater proportion of the population at equilibrium, so reductions to total abundance will be greater (essentially, the age classes represented in the second term on the right-hand side of Equation ( 8) account for a greater proportion of the population age structure).
Once the disturbance ceases, recruitment will return to the normal value of R. The expression for abundance at some time after that point, + Δ , will then have an additional term reflecting the new age classes arriving after the disturbance:
As the time after disturbance increases, the third term on the righthand side of Equation ( 9) will increase and the second term (representing the cohorts affected by disturbance) will decrease as the population convergences on the steady state given by Equation ( 7).
The abundances in that second term will decrease by a factor e -M with each progressive year, so the recovery trajectory will have a shape proportional to 1 -e -M . If we define recovery time as the time required to reach 95% of the initial population size, then the recovery time will be M -1 ln Impact ∕ (1 -0.95) , where Impact is expressed as the ratio of N to pre-impact abundance N 0 .
These dynamics are illustrated in the simulations in Figure 3, using population parameters for black rockfish (Table S1) as an example. The initial decline in abundance during the pulse disturbance is a geometric decay that lasts for the duration of the disturbance, followed by an inverse geometric recovery (1 - e -M ) to the initial state. The initial impact is greater for greater reductions in recruitment ( ) (Figure 3a), longer pulses of disturbance ( ) (Figure 3b), and populations with higher natural mortality rates (M; and thus shorter average lifespans; Figure 3c). The recovery time after disturbance is longer for greater initial impacts, if populations have the same demographic parameters (Figure 3a,b). However, the recovery time is faster for populations with higher mortality rates (shorter lifespans) despite the greater initial impact on abundance from the same relative reduction in recruitment (Figure 3c).
The patterns of impact and recovery time for biomass are similar (Figure 3d-f), but slower, with the maximum impact not realized until after the disturbance ceases, reflecting the delay in an affected recruit cohort becoming a large proportion of population biomass (Figure 2b). That delay also produces a more rounded and 'spread out' trajectory of decline and recovery, essentially because the maximum biomass of any given cohort occurs over a range of ages (Figure 2b), whereas the maximum abundance of a cohort occurs at age 1 and declines exponentially from there (Figure 2a).
Interpreting the open population model as the convolution of the pulse disturbance to recruitment and either the abundance-or biomass-at-age influence functions (Figure 2) sheds additional light on the expected shape and duration of the population response to disturbance. In general, the convolution of two peaked functions will (7)
yield a function with a peak that is approximately as wide as sum of the widths of the two input peaks, and a centroid (or center of mass) that is the sum of the centroids of the two input functions (Bracewell, 1965). This explains why the impact of a constant 'square wave' disturbance (green curves in Figure 4) produces a slower impact and recovery time for biomass than for abundance, because the influence function for biomass is both broader and has a centroid Expanding beyond these examples, we calculated impacts and recovery times for all 16 fish species in our dataset for several different levels of disturbance to recruitment during a 5-year pulse event.
To compare species, we present these results in bivariate plots. For a given level of disturbance, the species' responses fall along a nearly straight line, with the order along the line determined by natural mortality rate (M): lower M corresponds to low impact but long recovery time, and higher M corresponds to higher initial impact but shorter recovery time (Figure 5a), as in the blue rockfish example (Figure 3c).
The solution in Equation ( 9) and depicted in Figure 5a Curves in each panel represent different intensities of disturbance (a, d; 10%, 30% or 50% reductions in recruit survival), durations of disturbance (b, e; 5, 10 or 15 years) or values of the natural mortality rate, M (0.5, 1, and 1.5 times the baseline value of 0.18 year -1 ). The light blue curves are the baseline scenario and identical across panels. Simulations use life history parameters for black rockfish, Sebastes melanops (Table S1), with no harvesting. where the first term represents survivorship when the age-a cohort was younger, prior to disturbance, and the second term is survivorship during the disturbance. The entire solution would be cumbersome to write out, but it is straightforward to note that additional organisms lost in any given age class a in 1 year would be (assuming M is constant with age)
The overall effect of the disturbance on the population will be similar to the recruit-only case (but not identical), with greater impact for longer or more intense disturbances, and a recovery time that is correspondingly faster with higher M. However, once most or all of the age classes are affected by disturbance, there is no longer an effect of M on the initial impact, because M no longer determines the relative proportion of individuals in the affected cohorts. This is confirmed by examining simulations with a proportional increase in M applied to all age classes for the 16 fish species (Figure 5b): there remains a clear negative relationship between M and recovery time but little effect on impact. In fact for the highest values of M, the impact was actually slightly lower than for lower M (Figure 5b); this is because the relative difference between e -(M) and e -(M+ M) (i.e. the reduction in survival due to disturbance) decreases as M increases.
In a closed population, the R terms in the solution in Equation ( 9 Nonetheless, the closed model retains the strong effect of M on the relative impact and recovery time, with a pattern nearly identical to the open population scenario (Figure 5c).
If a population is harvested, then F in Equation ( 4) is nonzero, and the main effect on the population is a reduction in survival to old truncates the age structure but does not affect recruitment, so the effect of higher harvest is simply the same as higher natural mortality rates (M): the impact of a given disturbance is greater but the recovery time is faster, because the cohorts affected by disturbance make up a greater proportion of the age structure. This is illustrated in Figure 6 using blue rockfish as an example. By contrast, in a closed population, the effect of harvest on reproductive output is important, so scenarios with higher F have both greater impacts (as in the open population case) and slower recovery times (Figure 6).
We considered the relative effects of realistic levels of disturbance to either the youngest age class or all classes, or both, using a closed population model of Pacific cod. If the disturbance affects only the first year class, the response in terms of biomass (Figure 7a, red curve) would lag the disturbance (as depicted in the example in Figure 3d) due to the time required for the affected cohorts to accumulate biomass. By contrast, a disturbance only affecting adult mortality had an immediate effect on biomass (Figure 7a, blue curve), and was similar in magnitude to the impact felt by simulations with both types of disturbance together (Figure 7a, black curve). The latter suggests that the overall population impact of higher adult mortality is much more substantial than the reduction in recruitment, for the levels of disturbance we simulated.
The relative effects of the two disturbance types were also distinct in their effect on fishery yield (Figure 7b). Because cod do not reach harvestable size until they are approximately 4 years old, the impact of a disturbance on recruits only did not affect yield until after the 3-year disturbance had ceased, and reached its maximum impact 4 years later. Similarly, the yield when both disturbances occurred (Figure 7b, black curve) continued to decline for 2 years after biomass had already begun to recover (Figure 7a).
Our analysis provides new quantitative insights into what factors would lead to interspecific differences in two key components of resilience: magnitude of initial impact (also known as the resistance to impact) and recovery time (Capdevila et al., 2020;Ingrisch & Bahn, 2018). The baseline (pre-disturbance) natural mortality rate is an important determinant: species with a higher mortality rate are shorter-lived, so disturbances affecting only one or a few age classes have a greater proportional effect on total abundance.
However, those shorter-lived species also have faster transient dynamics during the recovery period, and are predicted to return to their pre-disturbance state faster. These insights about resilience dynamics allow us to predict how different species would respond to pulse-disturbance events. They are also consistent with observations from disturbed communities comprised of species with a range of longevities, such as benthic invertebrates in dredged seafloors (e.g. Rijnsdorp et al., 2018), and with analyses of structured population models of terrestrial plants, insects and vertebrates (Morris et al., 2008).
Many recent examples of pulsed disturbances affecting only a portion of the population age structure (e.g. only young-of-the-year) et al., 2011). Therefore, we also investigated how harvesting affected the population response to such disturbances. In general, increased harvest approximates increased natural mortality, so more heavily harvested populations would experience a greater reduction from a disturbance to a few cohorts. However, when most of the age classes are affected by disturbance, there is no longer an effect of either natural or harvest mortality on the level of reduction. In a closed system, higher harvest rates also lengthen recovery times because of lower per-capita reproductive output. This suggests that reducing harvest rates or protecting portions of the population in reserves (prior to disturbances) would be beneficial for resilience. Additionally, our simulations emulating the marine heatwave disturbance to Pacific cod revealed that harvested populations, the impact and recovery of harvest yield will differ from that of population biomass (assuming harvest remains constant during and after the disturbance). Specifically, the additional time lag required for new offspring to grow into the harvested portion of the adult population will likely prolong the impact on yield, even beyond the time at which overall biomass has begun to recover.
The natural mortality rate is the key demographic rate explaining variability among species in their response to disturbance. However, our general analysis does omit nuances that are important in particular applications. For example, in systems with tightly coupled consumer-resource or competitive species interactions, the nature of those interactions will also affect recovery trajectories (e.g. Commander & White, 2020). Alternatively, a system with multiple attractors (e.g. Allee effects) would not always recover to their original state if the disturbance pushed the system into a different attractor basin (Carpenter & Brock, 2006). Additionally, we examined disturbances that produced some specified reduction in survival for certain age classes. The actual link between an anomaly in some environmental variable and organisms' physiological tolerance or sensitivity to that change, that is, the reaction norm, is a separate, ecophysiological question. The nature of those reaction norms and the importance of the natural mortality rate in determining the response to disturbance could also be understood in the context of demographic buffering. The demographic buffering hypothesis posits that because fluctuations in demographic parameters tend to reduce the long-term population growth rate (Lewontin & Cohen, 1967), the parameters with greater elasticity on the overall growth rate will be under greater selection and should be less variable (Hilde et al., 2020;Pfister, 1998).
For the marine fish species we examined, variability in adult survival over time is much less than variability in larval or juvenile survival.
Because our analysis shows that resilience is most sensitive to the adult mortality rate, this could be taken as support for the demographic buffering concept, although we caution that our analysis is focused on ecological rather than evolutionary time scales.
Efforts to enumerate the occurrence of pulse disturbances typically define a pulse as a certain level of departure from normal conditions F I G U R E 7 Simulated disturbances like those experienced by Pacific cod in the 2014-2016 marine heatwave. Curves reflect simulations with a 42% reduction in recruitment (red curves), a 68% increase in adult mortality (blue curves), or both (black curves) during the three-year window indicated by grey shading. The response variable was (a) total population biomass or (b) fishery yield, both relative to their pre-disturbance levels. Simulations assumed closed demographics and used life history parameters for Pacific cod (Table S1). Frauenthal, 1986) but to our knowledge this is the first application of the idea to disturbance dynamics.
There is a long history of studies investigating population resilience in the context of returning to a stable age distribution after a perturbation (Ezard et al., 2010;Stott et al., 2011). Most such analyses have been conducted using linear models without density dependence. That approach has the advantage of opening up a wealth of analytical solutions and approximations (Botsford et al., 2019;Caswell, 2001) but in practice limits the applicability of the analysis to systems that have been driven to low abundance where densitydependence would be negligible. Additionally, in such models the recovery is to a state that is geometrically increasing (though with a stable age or stage distribution). Thus, it is possible that the estimated recovery time would be longer than the time it would take a population to grow in abundance enough for density dependence to become non-negligible. In other words, such models are best used to describe only short-term dynamics. By taking a nonlinear approach, we are able to better approximate pre-disturbance conditions at a steady state (and consider the scenario in which a very abundant population is disturbed but perhaps not to a density so low as to have linear dynamics) and capture the longer term convergence on a steady state rather than a constant growth rate. These are also the type of dynamics envisioned by the impact/recovery framework propsed by Ingrisch and Bahn (2018); in fact quantifying the 'impact' as we have done is not possible using a linear model. Additionally, while we did not explicitly consider the role of the 'strength' of density dependence in our analysis; we can implicitly learn about that by comparing our open and closed scenarios. Bearing in mind those differences between linear and nonlinear dynamics, analysis of linear matrix models (like those described by Ezard et al., 2010 andStott et al., 2011). has provided a number of metrics to characterize the transient recovery trajectory. Analogous calculations are not possible for nonlinear models with a non-zero steady state, such as ours. Nonetheless it is instructive to compare some insights between the two approaches. In linear matrix models, the logarithm of the damping ratio approximates the instantaneous rate of convergence of the population on the stable age distribution (Stott et al., 2011). The damping ratio is the ratio of the magnitudes of the first and second eigenvalues of the projection matrix. In general, the second eigenvalue will have greater magnitude (slower convergence) when the age distribution is truncated and reproduction is concentrated into a narrower band of age classes, such as when an iteroparous population is harvested (Botsford et al., 2019). This is consisent with our results for closed populations: the recovery time slowed as harvest increased.
A second well-known metric from linear structured models is
Cohen's distance (Cohen, 1979), which estimates the relative 'distance' of a population from its stable age/stage distribution (i.e. the time required for convergence). In a sense, that metric combines information on both the impact and the recovery time because the convergence time depends both on the severity of impact and on the dynamics of recovery. By separating those two aspects of resilience we were able to examine the factors affecting both, which can differ (e.g. the natural mortality rate may or may not affect the level of impact, depending on the range of age classes affected by disturbance).
On the other hand, an advantage of using linear models is the ability to perform sensitivity analyses on the eigendecomposition of the projection matrix to identify the demographic rates contributing most to resilience (Caswell, 2007;Morris et al., 2008). Nonetheless, we were able to use our analysis of the renewal equation to identify the natural mortality rate as a major determinant of both impact and recovery time. Our simulations then bore out that finding, with results from 16 species with a wide range of growth and maturity patterns falling out largely on an axis of variability determined by the mortality rate.
One concern in any examination of resilience and recovery time is that ecological systems are likely quite rarely at a true stable equilibrium and are rather in some stage of transient recovery towards an attractor, if not in a limit cycle or fluctuating chaotically (Hastings et al., 2018;Rapacciuolo et al., 2019). Thus, it is not necessarily clear when 'recovery' should be assessed in a real system, even if it has dynamics similar to those modelled here. Additionally, in the presence of ongoing stochastic variation, average population growth rates are slowed, which would further delay the recovery times modelled here (Lewontin & Cohen, 1967;Tuljapurkar, 1990).
Nonetheless, it is instructive as a first exercise to use deterministic analyses such as ours to understand the mechanisms governing dynamics in a simplified system. An extension of this analysis would be to ask how populations respond to repeated disturbances or environmental variability, either randomly distributed (e.g. white, red or pink noise; Vasseur & Yodzis, 2004) or with characteristic frequency content (e.g. El Niño-Southern Oscillation; Schmidt et al., 2018). In that case, the response should depend strongly on the generation time, because populations tend to amplify disturbances that occur at time scales similar to the mean age of reproduction (Bjørnstad et al., 2004;Botsford et al., 2019). That effect