SIMULATION OF POST-HURRICANE IMPACT ON INVASIVE SPECIES WITH BIOLOGICAL CONTROL MANAGEMENT

. Understanding the eﬀects of hurricanes and other large storms on ecological communities and the post-event recovery in these communities can guide management and ecosystem restoration. This is particularly impor-tant for communities impacted by invasive species, as the hurricane may aﬀect control eﬀorts. Here we consider the eﬀect of a hurricane on tree communities in southern Florida that has been invaded by Melaleuca quinquevervia (melaleuca), an invasive Australian tree. Biological control agents were introduced starting in the 1990s and are reducing melaleuca in habitats where they are established. We used size-structured matrix modeling as a tool to project the continued possible additional eﬀects of a hurricane on a pure stand of melaleuca that already had some level of biological control. The model results indicate that biological control could suppress or eliminate melaleuca within decades. A hurricane that does severe damage to the stand may accelerate the trend elimination of melaleuca with both strong and moderate biological control. However, if the biological control is weak, the stand is resilient to all but extremely severe hurricane damage. Although only a pure melaleuca stand was simulated in this study, other plants, such as natives, are likely to accelerate the decline of melaleuca due to competition. Our model provides a new tool to simulate post-hurricanes eﬀect on invasive species and highlights the essential role that biological control has played on invasive species management.

1. Introduction. The damage to tree communities periodically caused by hurricanes has played a crucial role in structuring the plant communities of southern Florida [10,14,24]. Understanding the effects of hurricanes and the process of posthurricane recovery in these communities, including both native and invasive species, can guide ecosystem restoration and subsequent management strategies. Studies frequently support the idea that large-scale disturbances accelerate the spread of invasive species [11,23], as more light reaches the ground layer following a hurricane. In addition, hurricanes may set up a positive invasive feedback loop where increased litter input to soil improves soil nutrient availability and moisture, subsequently favoring invasive plants [23,28,31]. Snitzer et al. [23] observed positive responses by invasive species to light the first year after a hurricane [23]. In contrast, other studies have shown that invasive forest species exhibit the same range of ecological roles as native forest species and compete with native species for the same kind of regeneration opportunities after a severe hurricane [1,12], consequently might not favor invasive species over native species to any great degree.
Most studies of hurricane effects on invasive plants have utilized examples where biological control is absent [3,16]. However, the use of biological control, using natural enemies to target and control invading plant species, is now integrated with several integrated plant management strategies [2,22,27]. In general, herbicidal and biological control methods are different in the sense that the former method usually requires repeated treatment of the invasive plants [1,8,13] while the latter is generally capable of being self-sustained on the targeted invasive plant species and may require only a single treatment [30]. Since biological control is qualitatively different from herbicide treatment, the recovery of invasive species after a hurricane may also be different.
Melaleuca quinqueneriva (Cav) S.T. Blake (Myrtaceae) (hereafter melaleuca) is a large, perennial tree imported to the United States from Australia beginning in the early 1900's [32]. During the 20th century, through natural expansion and several planting campaigns, melaleuca became one of the worse invaders of natural areas, including The Everglades, a World Heritage Site [32]. In addition to largescale herbicide and physical removal efforts, a suite of host-specific, phytophagous insects has been released to control melaleuca since 1996, are well established and have been contributing to the biological control of this invasive tree [5,7,17,25,26]. Herein, we aimed to understand the recovery of melaleuca following a hurricane event through modeling simulations.
We developed a size-structured matrix model that allows us to take into consideration the differential effects of the hurricane on various size classes of trees and simulate how invasive species respond to a hurricane. Hurricane disturbance has different effects on trees of different size, so stand age must be taken into account [10]. In particular, we assumed sensitivity to wind damage for canopy trees increases with the size of trees [23]. Therefore, in older stands, large trees, comprising a high fraction of basal area, are most likely to fall, shifting a stand towards younger trees, while younger stages will develop faster in open stands created by the hurricane.

Materials and methods.
2.1. Matrix model of the combined effects of biological control and hurricane. Here we describe a size-structured matrix simulation model incorporating the underlying mechanisms of stand dynamics for a pure melaleuca stand. The reason for using the matrix model rather than, say, an agent-based model, is that it is more convenient for quickly setting up new scenarios and parameterizations of biocontrol and hurricanes and evaluating the consequences than agent-based models that we know of. Matrix models are well suited to stands of a single species, as they can represent age-or size-structure in a mathematically simple way. An alternative approach, individual-or agent-based modeling, is more appropriate for mixed-species stands. The matrix model was an adequate representation here, as it was applied only to a pure melaleuca stand, for which unpublished empirical data over sixteen-year period are available to evaluate model accuracy. We performed the simulations in a small plot (0.01 ha), so that all individuals in the plot are affected by competition from larger trees in the simulated plot. In this sense, the model resembles the gap-phase models of forest succession [4].
A state transition matrix A, with Markov chain properties, is used in an equation that updates the stand structure each year; where, where Xi is the number of tree stems in size class i and where A multiplies the vector of tree number per diameter size class. For illustrative purposes, the matrix below represents 8 size classes, of diameters (0-5, 5-10, 10-15,...35-40 cm); X 1 (t + 1) X 2 (t + 2) X 3 (t + 3) X 4 (t + 4) X 5 (t + 5) X 6 (t + 6) X 7 (t + 7) However, the matrix model used in our simulations has a finer resolution of 32 size classes, so each of the size classes in the matrix model is subdivided into four size subclasses. Each time step of the equation is assumed to advance the growth process of the trees by 2 years. The parameters in the equation have the following meanings: a ii is the probability that the tree will still be in the same size class in the next time step (it has neither died or moved into a different size class during the year; a i+1,i is the probability that the tree will grow from size class i to size class i+1 in a time step. a ij (j≥j reproduction ) is the number of seedlings that become established that are contributed by each individual in size class i in time step, where j reproduction is the size class at which reproduction starts.
The matrix elements incorporate survival, growth and reproduction, which are assumed in the model to have the functional forms below.
Survival. Survival per time step for trees in size class i is expressed as where m 0 is a mortality rate and where d i is the mid-point of the size class. The exponential decay of mortality represents the assumption that natural mortality rate declines with size. The 1 ⁄2 power dependence on the diameter in the exponent is an assumption that intends, with the fitted parameter, k2, to allow mortality per time step to decline realistically with size.
Growth. To represent the probability of growing from one size class, i, to the next, i + 1, during a time step, we use The denominator reflects the assumption that asymmetric competition among trees slows the growth rate in a hierarchical fashion, with smaller trees being affected by all larger ones, but not by trees in the same or smaller size class. The model therefore assumes that trees of a given size class, i, are negatively affected by all trees that are larger, with the effect of larger trees being proportional to their diameter at breast height (DBH) raised to the 3 ⁄2 power. This represents some combination of shading and competition for space, water, and nutrients. Again, this is an assumption, but it provides reasonable stand simulation. We also make the assumption, in the numerator, that the rate of growth of diameter slows with age as a tree approaches maximum size, which is represented by the exponential decline. That is, the growth rate coefficient from one size class to the next to decline with size, although the fitted parameter, k 1 , is relatively small, making this decline gradual.
Fertility. We assume for simplicity that, starting at a certain size, j reproduction , the fertility of each size class has a value is f. Then the non-zero matrix elements are given by, where the forms taken by a ii and a i+1,i mean that mortality is assumed to occur first, followed by growth. The assumption of the same reproduction for all trees greater than 13 cm DBH (diameter at the breast height) is an assumption. The five parameters outlined above are used to calibrate it to data; m 0 , k 1 , k 2 , q, and f. We have fit these parameters to data from a nearly pure stand of melaleuca (M. B. Rayamajhi, unpublished) prior to any effects of biological control. Although we do not know the exact stand age, we assumed that it was several decades old. We estimated parameters by minimizing the sum of squares of the differences between the model values of tree size distribution after 200 years of simulation and the unpublished empirical data for a nearly pure melaleuca stand.
Biological control effects. Field studies indicate that biological control can reduce both growth rates and rates of seedling production, so these effects are incorporated in the model when biological control is started. The following replacements are made at the start of biological control; and where c p,i and c f,i are, respectively, effects of biological control on growth rates and reproductive rates of the melaleuca trees. Three different levels of biocontrol, incorporated in values of c p,i and c f,i were compared. The first corresponds to the level that is currently observed in the field, while the second is more moderate biocontrol and the third is weak biocontrol. Hurricane effects. When there is a hurricane, its effects are variable as they depend on the intensity, direction, etc., of the hurricane, as well as other factors. It is assumed here that the hurricane reduces the number of stems, but it also has effects on growth rate, resulting from the thinning of tree density in various size classes. In the year of the hurricane, the numbers of trees in size classes are changed as follows; Two different sets of parameters for hurr i were assumed, representing a moderate and a strong hurricane.
All parameter values are shown in Table 1. The parameters were based on the expectation that the trees follow certain life cycle patterns, such as increasing survival probability with size. The parameters m 0 , k 1 , k 2 , q, and f were fit by minimizing least squares around data on size distribution of melaleuca preceding biocontrol (M. B. Rayamajhi, unpublished data). The resultant parameters, with slight adjustments reasonably fit empirical data collected from two study sites on which the size distributions of melaleuca trees were followed over a period of years from 1997 through 2013 (M. B. Rayamajhi, unpublished data) on a 10 x 10 m plot of pure melaleuca in southern Florida. Sensitivity of the results in Figures 1 and 2 to the above five parameters was performed using simulations with +/-10 % of the baseline values in Table 1. The general nature of the results was similar to Figures  1 and 2 in all cases, although both the slope and the asymptotic value of the total basal area were sensitive to m 0 , k 2 . That is expected, as those parameters relate directly to tree mortality. Estimating these parameters more accurately should get more attention in the future.

2.2.
Simulations using a size-structured matrix model. The purpose of the modeling is to provide comparisons of different hypothetical combinations of biological control and hurricane damage. The matrix model of 32 size classes was used to project growth of a pure melaleuca stand impacted first by biological control and then by a hurricane, where both hurricane severity and given levels of effect of biological control could be adjusted. Simulations of three levels of biological control (weak, moderate and strong) were made, where the last level corresponds to what is believed to be the effect of the biological control agents currently being applied to melaleuca (parameters listed in Table 1). Two levels of hurricane severity were applied: moderate and strong. If a hurricane was imposed, it was assumed to occur 16 years after the start of biological control. For the moderate hurricane, it was assumed that the largest three size classes were reduced by approximately one-half, the middle two size classes by 30%, and only minimal damage occurred to the smallest size classes. The strong hurricane eliminated the two largest size classes, reduced the three next smaller size classes by 80% and reduced the smaller size classes by about one half ( Table 1). The strong hurricane was applied only in the case of biological control, as it did not have a qualitatively different effect than the moderate hurricane for the other two cases. The simulations were performed over a long enough time for the stand to reach steady state before biological control was imposed. The basal area of the stand reached a steady state by about year 200, starting from an initial 100 seedlings. The growth to a mature stand was relatively smooth, with some small oscillations due to the delay of about 25 years between seedling establishment and reproduction by that cohort (Figure 1).

3.
Results. Simulation scenario with strong biological control. In the model runs, biological control is assumed to be fully saturated at simulation year 200. Strong biological control reduced growth rate by about 85% and reproduction by about 50%. No direct effects on mortality were assumed, although the decline in reproduction could be assumed to include increased seedling mortality. Slowing growth rate exposed a cohort to the higher mortality rates of small size classes for a longer time period. Simulated size distributions after the imposition of biocontrol are shown for four, eight, 20, and 40 years, in Figures 2a, b, c, and d, respectively, after the start of biological control, showing biological control having a significant effect. In particular, the strong negative effect on reproduction creates a bimodal distribution of stems by depressing the first two size classes in these early years. This level of biological control, if sustained would lead to the elimination of the melaleuca. Because only a pure melaleuca stand was simulated here, other plants, such as natives, were not simulated, but their competition would likely accelerate the decline of melaleuca. Total basal area was computed in the simulation and shows a steady decline after the onset of biological control (Figure 3a).
Simulation scenario with strong biological control and moderate hurricane.A moderate hurricane was imposed on year 216 of the simulation; that is, 16 years after application of strong biocontrol. It increases the rate of decline of melaleuca for several years, although it does not have an obvious additive effect on the long-term decline (Figure 3b). A strong hurricane (not shown here) has larger, but qualitatively similar results. Simulation scenario with moderate biological control. It is interesting to study whether lower intensities of biological control would still have the effect of eliminating melaleuca. The assumed moderate biological control was a 40% reduction in growth rate and a 30% reduction in reproduction. Simulated size distributions are found for four, eight, 20, and 40 years, respectively, in Figures 4a, b, c, and d, after the start of biological control. Note that the numbers   of stems in the smallest size class do not decline very rapidly, as the effect of the biological control on reproduction is much less in this case than in the case of strong biological control. Even with this moderate level of biological control, melaleuca is on the path to eventual extinction in the modeled stand. The total basal area declines continuously, although not as rapidly as with the strong biological control (Figure 5a). Simulation scenario with moderate hurricane.A moderate hurricane imposed 16 years after application of moderate biological control accelerates the decline of   melaleuca for a couple of decades (Figure 5b). Again, as strong hurricane (not shown) causes a faster rate of decline, but does not have a great overall effect. Simulation scenario with weak biological control. A reduction of biological control to 20% reduction in growth rate and 10% reduction in reproduction still has a pronounced effect on the melaleuca. Here we show only the change in basal area, first without a hurricane, as shown in Figure 6a. The melaleuca declines rapidly at first, but then levels off at a lower, but stable level (Figure 6a). The reason for the stabilization at a lower level is that the matrix model incorporates density dependence, or compensation, in the growth rate due to shading and other competitive effects. Reduction in the coefficients of growth and reproduction cause a thinning of the canopy that reduces shading and allows faster growth and thus greater survival of the trees in smaller size classes, which compensates to some extent for the reductions due to biological control. The resilience of the melaleuca stand at its new reduced size can be further illustrated by imposing a hurricane. The melaleuca recovers well within a few decades from a moderate hurricane (Figure 6b). However, a strong hurricane reduces the melaleuca stand for nearly two centuries, showing the limits of resilience (Figure 6c). The biological control contributes to the delay in recovery, as can be seen from the effects of a hurricane on the stand when there is no biological control (Figure 6d). Although the effect on basal area is severe, recovery is faster than in the case where there is weak biological control.  4. Discussion. A size-structure matrix model was used to approximate the combined effects of biological control followed by a hurricane on a pure stand of the invasive melaleuca tree in southern Florida. The main result is that the current level of biological control, denoted as strong biological control here, if maintained, should virtually eliminate the melaleuca over several decades. A more moderate level of biological control also appears to be capable of eliminating the melaleuca, though over a longer time period. In both cases, a moderate hurricane that followed the imposition of biological control by 16 years steepened the decline for several years, but did not appreciably qualitatively affect the long term trend. The effects of an assumed weak biological control were different. In that case the density dependence in the model produced enough of a compensatory effect that the weak biological control could only reduce the melaleuca to a lower level. The reduced stand was resilient to a moderate hurricane, but a severe hurricane on top of the weak biological control suppressed the melaleuca for more than two centuries. The limitations of the size-structured model must be acknowledged. It used assumptions on size-specific growth, survival, and fecundity, as well as on the densitydependent mechanisms of shading and competition that affected these processes. These assumptions were made based on the general understanding of these processes, but there was little specific knowledge, so parameters were chosen that gave a reasonable fit to the main patterns of the existing data (M.B. Rayamajhi, unpublished data) to describe the existing empirical data. The model also assumed that the effectiveness of the biological control agents were not affected by the hurricane. This seems reasonable, as surveys of areas heavily impacted by Hurricane Irma still showed a that biological control agents were present in high densities. Therefore, we consider the results of the models to be first cut suggestions of the succession following a hurricane, though more detailed information is needed to provide reliable forecasts.
Although these limitations in the model must be kept in mind, the results may still be useful. The matrix model showed the effect of biological control on the steady state structure of a pure melaleuca stand, and it also showed that even weak biological control greatly reduced recovery of a stand to its pre-hurricane structure following a strong hurricane. These predictions will be interesting to test against future observations. Melaleuca dominance in the simulated habitats without biological control is due to its faster growth and higher reproduction rates compared to native trees [18,19], because it lacks natural enemies in southern Florida. Biological control reduced both growth and reproduction, previous work has used an individual based forest model, JABOWA, to show melaleuca largely declined in basal area [29], though its stem density remained high, because reproductive output, even with limitation was still relatively high. Consequently, the dramatic decline of melaleuca stimulated recovery of native species. When native species recolonized the habitat, they suppressed melaleuca further, since melaleuca is less competitive after long-term damage inflicted by biological control. For this reason, biological control suppression of melaleuca on plots in which native species occur may be much stronger than on pure melaleuca stands.
In general, major storm events are thought to favor invasive species by creating favorable disturbed conditions and dispersing propagules. For instance, Hurricane Katrina in Florida has caused high mortality of the local species and enhanced the invasion of invasive species [9]. These generalizations stress the spread of seeds; however, they do not take into consideration effects of top-down controls, which can play a suppressive role when present. Here, we simulated biological control's effect on melaleuca based on long-term empirical measurements, that biological control agents inflict strong defoliation effects on melaleuca, rendering them more vulnerable than native vegetation to breakage from hurricane-force winds [14]. Additionally, biological control vastly reduces melaleuca seed production [6], seedling establishments [21] and tree growth [20]. Simulations showed that hurricanes may have negative effects on melaleuca recovery in the presence of effective biological control. Introducing biological control in our simulations caused rapid reductions of density and basal area of melaleuca, supporting the hypothesis that intense wind disturbance, such as hurricane or other high windstorm may actually have a net negative effect on invasive species recovery when effective biological control agents are present [15]. This can now be tested in areas of southern Florida that were affected by Hurricane Irma in 2017. An empirical study is under way to estimate the damage of Hurricane Irma at sites occupied by melaleuca and relate the relative damage to the impact that biological control agents Oxyops vitiosa Pascoe and Boreioglycaspis melaleucae (Moore) and Lophodiplosis trifida Gagné have had on stands. Data from the empirical study will be used in a computer simulation model to project the long term impact of the hurricane damage on the effectiveness of biological control.

5.
Conclusions. Our model simulations support the hypothesis that biological control at its present level will suppress the melaleuca. Hurricanes can steepen the decline of melaleuca on which biological control is imposed but might not change the general trend of decline too much. Biological control will and greatly delay recovery following a severe hurricane to pre-hurricane total basal area.