Simulated impacts of insect defoliation on forest carbon dynamics

Many temperate and boreal forests are subject to insect epidemics. In the eastern US, over 41 million meters squared of tree basal area are thought to be at risk of gypsy moth defoliation. However, the decadal-to-century scale implications of defoliation events for ecosystem carbon dynamics are not well understood. In this study, the effects of defoliation intensity, periodicity and spatial pattern on the carbon cycle are investigated in a set of idealized model simulations. A mechanistic terrestrial biosphere model, ecosystem demography model 2, is driven with observations from a xeric oak–pine forest located in the New Jersey Pine Barrens. Simulations indicate that net ecosystem productivity (equal to photosynthesis minus respiration) decreases linearly with increasing defoliation intensity. However, because of interactions between defoliation and drought effects, aboveground biomass exhibits a nonlinear decrease with increasing defoliation intensity. The ecosystem responds strongly with both reduced productivity and biomass loss when defoliation periodicity varies from 5 to 15 yr, but exhibits a relatively weak response when defoliation periodicity varies from 15 to 60 yr. Simulations of spatially heterogeneous defoliation resulted in markedly smaller carbon stocks than simulations with spatially homogeneous defoliation. These results show that gypsy moth defoliation has a large effect on oak–pine forest biomass dynamics, functioning and its capacity to act as a carbon sink.


Introduction
Recent studies using atmospheric CO 2 observations and inverse modeling have inferred a terrestrial carbon sink of several petagrams of carbon (Pg C) per year over the past Content from this work may be used under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. few decades (Le Quéré et al 2009, Beaulieu et al 2012, and forest inventory-based analyses have identified a net sink in boreal and temperate forests of 1.2 Pg C yr −1 since 1990 (Pan et al 2011). However, there is substantial regional and temporal variability in this C sink. In temperate ecosystems, this variability is at least in part due to disturbance, including insect epidemics (Kurz et al 2008a, 2008b, Amiro et al 2010, Pfeifer et al 2011, Hicke et al 2012. As one example, gypsy moth (Lymantria dispar L.) outbreaks are a risk throughout the northern hemisphere, including Canada, Europe and Japan  (Johnson et al 2005). Over 41 million m 2 of basal area in the northeastern US are thought to be at risk of attack by gypsy moth, and thus far gypsy moths have invaded about 23% of the forested area that provides suitable habitat for them throughout the US (Hicke et al 2012). Field observations and remote sensing are increasing our ability to detect these attacks, and to assess the implications for ecosystem structure and functioning (Liebhold et al 1992, Barron and Patterson 2008, de Beurs and Townsend 2008, Schäfer et al 2010. Such observations are critical for the development and testing of predictive, mechanistic numerical models of the terrestrial biosphere. At one measurement site where a gypsy moth defoliation event was well-observed, numerical experiments that did not include gypsy moth defoliation failed to realistically simulate the observations of CO 2 fluxes and tree growth and mortality rates (Miao et al 2011, Scheller et al 2011. Here, we seek to build on our understanding of forest responses to gypsy moth defoliation. Because there has been relatively little research on the impacts of gypsy moth defoliation on the carbon cycle, we work in an idealized setting, and focus on understanding the chains of events, processes and mechanisms, rather than attempting to predict the specific future trajectory of a forest. Consequently, we ignore potentially confounding factors including climate change, fires and land-use change. The test-bed for our analysis is the Silas Little Experimental Forest (39.9 • N, 74.6 • W), chosen because gypsy moth defoliation was observed in detail at this site (Schäfer et al 2010. We describe how a terrestrial biosphere model was modified to enable it to represent defoliation events and then evaluate the model against eddy-flux and biometry measurements. We use the model to address the following questions. (1) What is the response of the carbon cycle to defoliation events of varying intensity? (2) How does the carbon cycle respond to the periodicity of defoliation events? (3) Is the spatial representation of defoliation important?

Methods
Our analysis was underpinned by recent observations at the Silas Little Experimental Forest (SLEF), located in Pemberton Township in the Pine Barrens of southern New Jersey, USA (39.9 • N, 74.6 • W). Detailed descriptions of the site already exist (Skowronski et al 2007, Schäfer et al 2010. The site's mean annual temperature is 11.5 • C, and its mean annual precipitation is 1123 mm. The soil is derived from Cohansey sands, and has low nutrient content and cation exchange capacity (Schäfer 2011 In 2007, the site was completely defoliated by gypsy moth from the end of May until mid-July , 2012, Schäfer et al 2010, Schäfer 2011. Following the cessation of herbivory, leaf area recovered to about 50% of the summer maxima that had occurred in 2005 and 2006. This defoliation event was probably unusually intense because the pines (then comprising about 27% of the live basal area) were defoliated and the plastic flagging on some of the trees was eaten. In the few years following defoliation, approximately 1/3 of the oaks (on a basal area basis) experienced mortality (table 1), although the different oak species varied in their ability to recover from the attack (Schäfer 2011). The pines exhibited little mortality following defoliation (table 1), and their recovery may have been facilitated by epicormic budding. Other areas of the New Jersey Pine Barrens were defoliated in 2007. At pine-oak and pine-scrub oak stands, oaks but not pines were defoliated . Other studies have also indicated that oak is preferred over pine by the gypsy moth (Campbell and Sloan 1977, Davidson et al 2001, Barron and Patterson 2008.
Our modeling was carried out with the ecosystem demography 2 model (ED2) (Medvigy et al 2009(Medvigy et al , 2010. The model uses a system of size-and age-structured partial differential equations to closely approximate the first moment behavior of a corresponding individual-based stochastic gap model to realistically represent fine-scale heterogeneity in ecosystem structure within each grid cell. ED2 differs from most other large-scale terrestrial models by formally scaling up physiological processes through vegetation dynamics to ecosystem scales, while simultaneously modeling fires, natural disturbances, land use and the dynamics of recovering lands (Moorcroft et al 2001, Hurtt et al 2002, Albani et al 2006. Several important modifications were made to ED2 for this study. First, although ED2 has previously been shown to accurately simulate eddy-flux measurements and tree growth and mortality measurements throughout the northeastern US  (Medvigy et al 2009, Medvigy andMoorcroft 2012), the model was never specifically evaluated for xeric sites such as SLEF, and did not include plant functional types (PFTs) for xeric habitats. Unsurprisingly, preliminary analysis indicated that ED2 did not capture net ecosystem productivity at SLEF as accurately as for mesic sites. We addressed this by introducing two new PFTs, a xeric hardwood and a xeric conifer. Whittaker and Woodwell (1968) derived allometric relationships for the species characteristic of the xeric northern coastal plain, and we use their Quercus coccinea allometry for our xeric hardwood and their Pinus rigida allometry for our xeric conifer. With the additional exceptions of the parameters noted in table 2, the parameterizations for the xeric hardwood and xeric conifer were set to be identical to ED2's existing parameterizations for mesic mid-successional hardwoods and mesic northern pines, respectively. A second new aspect developed in this work was the characterization of a defoliation event in ED2, done through ED2's already-existing capability of representing sub-grid-scale heterogeneity related to disturbance (Albani et al 2006, Medvigy et al 2009. For example, suppose that the landscape initially consists of N sub-grid-scale tiles. To simulate a defoliation event, each of the N tiles is split into two daughter tiles. One of the tiles is set to be of size f relative to its parent tile, and the other daughter tile is set to be of relative size 1 − f . The tile of relative size f will undergo complete defoliation of all trees (see below), regardless of tree size, age, or previous disturbance history. The tile of relative 1 − f remains completely undisturbed throughout the defoliation event. We subsequently refer to f as the 'defoliation intensity'. This procedure has the effect of doubling the number of tiles at each defoliation event, and this can lead to a computationally unwieldy number of tiles. The model therefore periodically reviews all tiles and merges tiles with similar forest composition and structure.
The temporal evolution of defoliation events was derived by using observations from SLEF. Although only measurements of maximum seasonal leaf area index (LAI) are available, sap flow and LiDAR observation, ongoing since 2005, can be used to infer the seasonal trajectory of LAI (Clark et al 2012). One possibility would be to prescribe the inferred 2007 LAI time series (figure 1) for any simulated year in which SLEF experiences a gypsy moth outbreak. However, in long-term simulations spanning years and decades, there will be changes in the simulated structure and composition of SLEF, as well as the maximum potential LAI. To deal with the possibility of changes in maximum potential LAI, we defined the relative LAI as the ratio of the 2007 LAI to Of central interest here is the model's parameterization of mortality. The central unit of ED2 is the 'cohort': the density of trees of a given plant functional type, size (with separate leaf, fine root, sapwood and heartwood biomass pools), and disturbance history. The 'forest' in the model typically consists of a collection of cohorts that dynamically evolves over time. Mortality, applied on a monthly time step, acts by decrementing cohort number density. The mortality rate depends on a cohort's carbon balance. When a cohort is in positive overall carbon balance, mortality proceeds at a small, background rate. However, mortality rates can become large if the cohort is in negative or nearly negative carbon balance. In the model, the gypsy moth does not directly kill trees. Instead, defoliation is modeled by reducing the size of an affected cohort's leaf biomass pool. If the cohort consequently falls into an unfavorable carbon balance, it will then experience large mortality rates. For further details on mortality in ED2, see Albani et al (2006).
In a preliminary analysis, we found that this parameterization enabled the model to give satisfactory simulations of net ecosystem productivity, tree diameter growth and tree mortality for the 2005-9 period at SLEF, which included the All simulations were carried out for a single model grid cell. The grid cell is intended to represent the approximate footprint of the SLEF eddy-flux tower (1 km 2 ), but the size of the grid cell does not affect model results. The meteorological forcing, initial stand structure and composition, and soil texture are assumed to be the same for the entire grid We carried out 3 sets of simulations to investigate the effects of defoliation intensity, defoliation periodicity and spatial pattern of defoliation. (i) Intensity experiments: the fraction of the grid cell that was defoliated varied among simulations, and ranged from 0% to 60% (see table 3). For simplicity, we assume that defoliated trees become completely defoliated. Also, because hardwood trees have been most severely affected by historical gypsy moth outbreaks (Davidson et al 2001, Barron andPatterson 2008), we assume that only the hardwoods are defoliated. In each of these experiments, defoliation was prescribed to occur every 5 yr. This is consistent with observational records, which indicate a 5 yr periodicity in 1 km × 1 km quadrats in xeric oak-pine forests of the northeastern US (Johnson et al 2006). (ii) Periodicity experiments: defoliation was prescribed at 40% intensity in each simulation, comparable to historical defoliation events (Johnson et al 2006), but the periodicity of defoliation varied from 5 to 60 yr. We again assume that all defoliated trees are completely defoliated, and that only the hardwoods are defoliated. (iii) Spatial pattern: we investigated if a partial defoliation of all trees in the grid cell would differ from complete defoliation of some trees and leaving others completely intact. Depending on the size of the grid cell, both situations can be realistic. At the level of the forest stand, outbreaks are spatially heterogeneous and often out-of-sync (Johnson et al 2006), while on large spatial scales outbreaks are more synchronous (Peltonen et al 2002). This comparison is also of interest because many terrestrial biosphere models run at coarse resolutions (e.g., 100-200 km) and do not represent spatial heterogeneity very well. A potential way for these models to represent defoliation would be a partial defoliation of all trees in the coarse grid cell. This aggregation is also apparent when attempting to use spectral reflectance and derived data products (NDVI, GPP, etc) from satellite observations with low spatial resolution, such as MODIS. We therefore created an alternate defoliation prescription consisting of a weighted average of the 2005 (60%) and 2007 (40%) LAI profiles, and applied this to all trees in the grid cell with a 5 yr defoliation periodicity.
Simulations are named according to their intensity and periodicity; e.g., 40% defoliation with 5 yr periodicity would be named I40-P5. All simulations use the default heterogeneous representation of defoliation unless explicitly suffixed by '-HOM' to denote homogenous representation of foliage. This notation is summarized in table 3.

Intensity of defoliation
In all simulations, long-term net ecosystem productivity (NEP; equal to photosynthesis minus respiration) takes on positive values for well over a century, indicating that SLEF continues to act as a carbon sink regardless of defoliation ( figure 2(a)). The transient behavior is similar in the different simulations and is characterized by rapidly increasing NEP over the first 10-20 yr, stabilizing, and finally a slow decline beginning by about year 60. There is a clear ordination of the NEP time series corresponding to defoliation intensity, with I0 generally having the largest NEP and I60-P5 the smallest NEP ( figure 2(a)). Simulation-averaged NEP decreases approximately linearly with disturbance intensity at a rate of about 7 g C m −2 yr −1 for every 10% of defoliation. Note that the figure shows NEP averaged in 6 yr bins to enhance the legibility of the long-term trends. NEP does take on negative values in some years, especially those years with large defoliation intensity.
Our results can be compared to Edburg et al (2011), who simulated forest recovery from a mountain pine beetle outbreak. Whereas the simulations presented here prescribe periodic defoliation, their simulations prescribed a single mortality event near the simulation start time. They found that 50% mortality would result in an NEP anomaly of about −50 g C m −2 y −1 in the decade following the outbreak. This anomaly is quite similar to the approximately −70 g C m −2 y −1 anomaly between our I50-P5 and I0 despite the aforementioned differences in simulation design.
We find a nonlinear relationship between aboveground live biomass (AGB) and defoliation intensity ( figure 2(b)). The 200 yr time series from I0 and I10-P5 are generally closely matched, but then there is a wide gap between I10-P5 and I20-P5, and finally the simulations with the largest defoliation intensities (I40-P5, I50-P5, I60-P5) are generally clustered together. In the model, changes in AGB are determined as the balances between gains (tree growth and recruitment) and losses (mortality). We examined the time series of growth, recruitment, and mortality, and found that the mortality term (figure 2(c)) had the largest nonlinear response to defoliation intensity. Simulation-averaged mortality was nearly equal in I0 and I10-P5, increased as intensity increased from 10% to 40%, and then increased at a slower rate as intensity increased from 40% to 60%. This is consistent with the differences in live AGB among simulations ( figure 2(b)). Interestingly, Baker (1941) also showed a nonlinear increase in mortality with defoliation intensity, however in the cases studied there, the response was not as gradual as the modeled response in this investigation (figure 2(c)).
We also identified interactions between defoliation intensity and drought. Mortality rates during years 30-89 are shown in figure 2(d). In all simulations with defoliation, there is a spike in mortality every 5th year, in accord with the disturbance frequency. However, I0 has spikes in mortality in every 6th year, corresponding to years with meteorological forcing based on 2010 observations. That summer was New Jersey's hottest and 8th driest since records began in 1895. 5 Because I0 had the smallest late summer evapotranspiration per unit leaf area (figure 2(e)), it had warmer leaves than the other simulations. The higher foliar temperatures led to increased water stress, and the water stress ultimately led to increased mortality. This is consistent with a previous study that investigated the interactive effect of drought and defoliation (Bréda and Badeau 2008).
Insect defoliation has a strong effect on simulated ecosystem composition. We expected that the xeric conifers would eventually die out in all simulations because they require fire for germination, and our simulations did not include fire. Without insect disturbance, conifers indeed become locally extinct by around year 85 of the simulation (figure 2(f)). However, because the defoliation selectively impacts the oaks, it gives a competitive advantage to conifers, and the persistence of conifers is much stronger in the simulations with defoliation. Conifer biomass actually increases as a proportion of the total for the first few decades in simulations with defoliation intensities of 30% or larger (figure 2(f)). Some historical records have shown similar behavior with differential species responses to climatic and biotic events, thus changing species composition over time (Albertson and Weaver 1945). Nevertheless, the simulations probably overstate the actual ability of pines to capitalize on oak defoliation because pines are more susceptible to windthrow and lightning when the canopies of the larger oaks are reduced, and this is not accounted for in the simulations. Furthermore, defoliation of pines has been observed to occur during unusually extreme outbreaks in the New Jersey Pine Barrens, and this was not accounted for by the model.

Periodicity of defoliation
There was a strong nonlinear response of the carbon cycle to defoliation periodicity. When averaged over the 200 yr simulation, both gross primary productivity (GPP) and total ecosystem respiration (R tot ) responded positively to increasing periodicity, but the stronger response of GPP drove increases in NEP ( figure 3(a)). NEP increases by over 40% as defoliation periodicity increases from 5 to 15 yr, but additional increases in periodicity have a comparatively small effect on NEP. This model prediction suggests that there is a compounding effect of defoliation that is mainly operative for periodicities shorter than 15 yr, and is similar to previous findings that individual oak trees in New England require about 10 yr to recover from a single heavy defoliation event (Campbell and Sloan 1977). However, at nutrient-poor sites, the timescale to recovery may be longer if there is a lag before canopy nitrogen content returns to pre-defoliation levels. Such nutrient dynamics were not implemented in our simulations.
The end-of-simulation live AGB values also increased by over 40% as defoliation periodicity increased from 5 to 15 yr, and end-of-simulation live AGB was insensitive to additional increases in periodicity ( figure 3(b)). However, increasing defoliation periodicity from 15 to 60 yr was associated with increases in the live AGB in the largest trees, and decreases in the live AGB of intermediate size classes. Furthermore, increasing defoliation periodicity was also correlated with decreases in conifer fraction; the I40-P5 end-of-simulation conifer fraction was 6% (compare figure 2(f)), while in I40-P60, conifers became locally extinct by about year 115.

Spatial structure of disturbance
Comparing homogenous (i.e. suffix HOM in the simulation, figure 4) versus heterogeneous defoliation events whereby a different fraction of foliage is consumed showed a large spatial and temporal variability. The 'I40-P5-HOM minus I40-P5' difference in annual NEP is presented in figure 4(a), and is a substantial fraction of total NEP in I40-P5 (see figure 2(a)). A spectral analysis found that most of the power was in the 5 yr and 6 yr cycles, which are associated with the defoliation and meteorological forcing, respectively.  However, I40-P5-HOM consistently has larger NEP than I40-P5, especially at the beginning of the simulation. Both runs have NEP values approaching zero toward the end of the simulations as the ecosystems approach their respective equilibria.
Throughout most of the simulated period, I40-P5-HOM has 20-30% more live AGB than I40-P5 ( figure 4(b)). This arises because I40-P5-HOM has both greater GPP than I40-P5 (by 7% on the simulation average) and less mortality (by 44% on the simulation average). At the level of the entire canopy, the difference in GPP arises because of the nonlinear, saturating relationship between GPP and LAI. However, I40-P5 also has a distinct mortality peak every 5 yr following each defoliation event that is absent in I40-P5-HOM. In I40-P5-HOM, all oaks carry out some photosynthesis with >60% leaves throughout the defoliation event compensating for the lost foliage (VanderKlein and Reich 1999). After defoliation, the lost leaf biomass can in part be made up from the year's photosynthate and from reallocation. In contrast, 40% of the oaks in I40-P5 cannot do any photosynthesis whatsoever during defoliation. After defoliation, these oaks have not accrued photosynthate to draw on, reallocation from roots must fill a larger gap, and many oaks fall into negative carbon balance and may experience mortality. Moreover, the absence of carbon reserves may weaken oaks and allow for fungal attack. In the field, many oaks that died after the 2007 defoliation event at SLEF were observed to be ringed by basiocarps in fall 2009.

Conclusions
This study investigates a spectrum of ways in which periodic gypsy moth outbreaks can affect forest biomass dynamics, functioning and composition on annual to century time scales. We address this issue by analyzing idealized simulations from a structured terrestrial biosphere model capable of representing the small-scale heterogeneity arising from disturbance, and we focus on SLEF as a case study because of its strong observational record. Few previous mechanistic models have successfully simulated changes in carbon cycling during these disturbances, and so one of our objectives here was to establish a baseline that can be evaluated in future studies.
We found that NEP decreased approximately linearly with increasing defoliation intensity, but live AGB and ecosystem composition exhibited strong nonlinearities. In the idealized scenario without defoliation, the ecosystem was more sensitive to drought and experienced more drought-induced mortality than in the more realistic scenarios that included defoliation. We also found a nonlinear response to disturbance periodicity for both NEP and live AGB, with ecosystem sensitivity generally being much larger for 5-15 yr periodicities than for 15-60 yr periodicities.
A major implication of this study is that it is essential to correctly specify the spatial and temporal patterns of defoliation events in order to accurately simulate the corresponding carbon dynamics. This is relevant to both historical and future outbreaks, and becomes particularly important when considering large spatial scales because regionally-averaged defoliation amounts can be much less than local defoliation amounts in patches. Partial defoliation of all stands was simulated to have a much weaker effect on site carbon budgets than a complete defoliation of some stands. In the case where each tree 'pays' for defoliation with the same fraction of its leaves, net carbon uptake is ∼20% larger than the case where some trees 'pay' nothing and other trees are completely defoliated. This result highlights a challenge faced by aggregated terrestrial biosphere models that do not represent biotic heterogeneity resulting from disturbances, and thus overestimate NEP at larger spatial scales.
This study included idealized simulations to glean some insights into the century-scale impacts of frequency and intensity of gypsy moths defoliation on the carbon cycle. Actual ecosystem responses will depend on a range of complex factors including future drought frequency and severity, changes in the frequency of defoliation, and interactions with other disturbances, such as fire and windthrow. There is also uncertainty in the applicability of some of the details of our modeled defoliation events. For example, we prescribed the same temporal evolution of all defoliation events (figure 1). This was because we only had measurements at a sufficiently high temporal resolution for a single oak-pine forest defoliation event (SLEF in 2007). For simplicity in the model formulation, we also assumed that all oak trees in a defoliated stand share the same defoliation rate, while the pines are not defoliated at all. In future work, our model can be generalized to allow trees of size classes to be defoliated at different times, or to allow for the possibility of pine defoliation. Finally, our results were driven by information from one particularly well-measured site. In ongoing work, we are carrying out additional simulations to assess the regional-scale impacts of gypsy moth defoliation on the forests throughout the eastern US.