Development of a New TRIPLEX-Insect Model for Simulating the Effect of Spruce Budworm on Forest Carbon Dynamics

The spruce budworm (SBW) defoliates and kills conifer trees, consequently affecting carbon (C) exchanges between the land and atmosphere. Here, we developed a new TRIPLEX-Insect sub-model to quantify the impacts of insect outbreaks on forest C fluxes. We modeled annual defoliation (AD), cumulative defoliation (CD), and tree mortality. The model was validated against observed and published data at the stand level in the North Shore region of Québec and Cape Breton Island in Nova Scotia, Canada. The results suggest that TRIPLEX-Insect performs very well in capturing tree mortality following SBW outbreaks and slightly underestimates current annual volume increment (CAI). In both mature and immature forests, the simulation model suggests a larger reduction in gross primary productivity (GPP) than in autotrophic respiration (Ra) at the same defoliation level when tree mortality was low. After an SBW outbreak, the growth release of surviving trees contributes to the recovery of annual net ecosystem productivity (NEP) based on forest age if mortality is not excessive. Overall, the TRIPLEX-Insect model is capable of simulating C dynamics of balsam fir following SBW disturbances and can be used as an efficient tool in forest insect management.

Defoliation varies with SBW population density and can be measured in terms of annual defoliation (AD) incurred.SBW mainly consumes foliage produced in the current year.Given that needles typically remain on the tree for five to six years, it can take many years for the SBW to kill a tree, especially at moderate or low defoliation levels (i.e., when not all foliage produced within a given year is consumed) [13].The last SBW outbreak started in the 1970s and 1980s in eastern North America, continued for almost two decades and was considered extremely severe (e.g., 80%-100% annual defoliation) [14].Tree defoliation and ultimately mortality vary greatly between stands during an outbreak as well as between outbreaks themselves [15].Older stands are more vulnerable than younger, more vigorous stands [16].Generally, C dynamics and insect disturbances lead to temporary changes in environmental conditions.However, insect outbreaks can cause long-term changes to ecosystem structure and function [17][18][19], which can result in long-term impacts in C dynamics.We therefore need to better understand how insect outbreaks interact with stand structure (e.g., how defoliation severity and duration interact with stand age) to improve our estimates of forest C losses and our understanding of regional and global C cycles.
Modeling is a complementary tool to traditional observation and experimental approaches that can provide long-term and large spatial-scale perspectives on the effects of forest insect epidemics as well as evaluate interactions among multiple factors [20].In recent years, there are some studies that have been conducted on C dynamics following insect outbreaks [21,22].For example, Albani et al. [23] used an ecosystem demography model in conjunction with a stochastic model to predict the impact of hemlock woolly adelgid (Adelge tsugae, Annand) on C dynamics in the eastern United States.Medvigy et al. [24] also simulated a linear decrease in net ecosystem productivity (NEP) with increasing defoliation intensity by the gyspsy moth (Lymantria dispar L.).In addition, Edburg et al. [25] simulated coupled C and nitrogen (N) dynamics following mountain pine beetle (Dendroctonus ponderosae (Hopkins), Coleoptera: Curculionidae, Scolytinae) outbreaks in the western United States of America.They found that the severity of the outbreak itself was an important factor which led to the initial decline in NEP, and required 80 to 100 years or more following the disturbance to recover.However, because CLM4 does not explicitly contain an age-class distribution model for trees, the impact of tree age on C dynamics following insect outbreaks is not known.In Canada, the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3) has been used to simulate and forecast ecosystem C dynamics impacted by SBW disturbances [10,26,27].Such studies found that a decrease in NEP may convert boreal forests from C sinks into C sources after SBW outbreaks.However, CBM-CFS3 is an empirically driven model that is not designed to estimate gross primary productivity (GPP) and respiration in forest ecosystems.It is therefore unable to explain in detail the flow of C dynamics under such conditions.
TRIPLEX, a process-based hybrid model [28,29], integrates the advantages of both empirical and mechanistic components.TRIPLEX 1.0 is able to simulate both short-and long-term forest growth and C and N dynamics and has been successfully calibrated and validated for different forest ages [28,30], tree species [31,32], and harvest disturbance types [33,34] in boreal ecosystems.The main objective of our study was to develop a new process-based model (TRIPLEX-Insect) to quantify forest mortality and C sequestration in response to SBW outbreaks at the stand level.Specifically, this study (1) provides and tests a novel approach to integrate the relationship between defoliation and mortality, and (2) confirms that our sensitivity experiments contribute to the evaluation and understanding of how defoliation affects the capacity of C sequestration in forest stands.

Study Site Description
The study area is located within the boreal forest of the North Shore region of Québec, Canada.The dominant species are black spruce (Picea mariana (Mill.)B.S.P.) and balsam fir, which are both host species of the SBW [16].In 2006, the first signs of a SBW outbreak were observed in the North Shore region of Québec, and subsequent tree mortality has been observed in recent years [6].Three sample plots of mainly balsam fir were established in 2006 (see Table 1) in the epicenter of the outbreak to study long-term forest response to the current SBW outbreak.We sampled AD on a 45-cm branch segment sampled in the mid-crown [35] of each tree (20 balsam fir per stand) in each plot in August of each year.We assessed AD in the laboratory using the Fettes method [36] for each year since 2006 [37].
We obtained monthly climate data, including air temperature, precipitation, and relative humidity from 1913-2016 from the CAUSAPSCAL (48.37 • N, 67.23 • W) meteorological station.Environment Canada has published the relevant climate data on their website (http://climate.weather.gc.ca/).We also used monthly data from 30-year averages (1913-1942) to replace earlier (1805-1912) climate data.Our study also used published data including plot locations and methods of annual assessments of defoliation in various regions of Cape Breton Island as well as density, basal area, volume, DBH and height per plot [38,39].

Model Structure
TRIPLEX 1.0 integrates three well-established models: 3-PG [40], TREEDYN3.0[41], and CENTURY4.0[42].TRIPLEX 1.0 can simulate forest growth and key C and N cycling processes which mainly include forest production, particularly GPP, net primary productivity (NPP), and NEP (see Figure 1) [28,30,31].Model simulations were conducted using a monthly time step, and monthly GPP was calculated from photosynthetically active radiation (PAR), leaf area index (LAI), forest age, mean air temperature, soil water, and the percentage of frost days per month.Boreal forest NPP was estimated using Equation (1): where Ra is autotrophic respiration, which is estimated using N, air temperature, and component C pools (e.g., wood, branches, foliage, coarse roots, and fine roots) [30].NEP is calculated as GPP subtracted by ecosystem respiration (ER), which is comprised of Ra and heterotrophic respiration (Rh).Rh is calculated by subtracting root respiration from soil respiration, expressed as the exponential function of temperature and Q 10 [30].Annual values are estimated by summing monthly results.
In TRIPLEX 1.0, soil C and N are based on the CENTURY soil decomposition sub-model [42].This sub-model provides realistic estimates of both C and N mineralization rates for Canadian boreal forest ecosystems [43].More detailed information on TRIPLEX 1.0 is described in Peng et al. [28] and Liu et al. [29].Although the insect sub-model is designed to resolve how C dynamics respond to SBW outbreaks, many potentially confounding factors, such as climate change, forest fires, and land-use changes, are not directly considered.[25]).Rectangles represent pools or state variables; ovals represent simulated processes; diamond-shaped boxes represent judgment process; dotted lines represent controls; and solid lines represent carbon (C), nitrogen (N), and water flow.

Tree Mortality Resulting from Defoliation
Observed annual defoliation (AD) can be considered as a measure of the percentage of leaf biomass lost to insects.Since this variable causes a decrease in LAI, it subsequently impacts forest ecosystem photosynthesis.Every year, SBW produces one generation that is divided into nine life stages comprised of eggs, larvae (which include stage one through six), pupa, and moth.Damage from the SBW mainly occurs from May to July when the larvae are in stage three through six [44][45][46].Temperature impacts the beginning of SBW feeding.We estimated that defoliation commenced in May when the mean temperature is generally greater than 5 °C in the study area [44,47,48].Thus, AD was distributed over a period of two to three months by setting different weights to AD in our model (see Table 2) to account for differences in SBW phenology [49].[25]).Rectangles represent pools or state variables; ovals represent simulated processes; diamond-shaped boxes represent judgment process; dotted lines represent controls; and solid lines represent carbon (C), nitrogen (N), and water flow.

Tree Mortality Resulting from Defoliation
Observed annual defoliation (AD) can be considered as a measure of the percentage of leaf biomass lost to insects.Since this variable causes a decrease in LAI, it subsequently impacts forest ecosystem photosynthesis.Every year, SBW produces one generation that is divided into nine life stages comprised of eggs, larvae (which include stage one through six), pupa, and moth.Damage from the SBW mainly occurs from May to July when the larvae are in stage three through six [44][45][46].Temperature impacts the beginning of SBW feeding.We estimated that defoliation commenced in May when the mean temperature is generally greater than 5 • C in the study area [44,47,48].Thus, AD was distributed over a period of two to three months by setting different weights to AD in our model (see Table 2) to account for differences in SBW phenology [49].Defoliation caused mortality (DM; %/year) is defined in our model and described by Equation ( 2): where K = 6.6774 is a parameter derived from empirical data taken from literature [50].The defoliation stress (DS) index was calculated by cumulative defoliation (CD) and the weight of current annual defoliation (WAD) using Equation ( 3): If annual current defoliation levels for the last five years were 0.25, 0.25, 0.5, 0.75, and 1, then CD would be 0.25, 0.5, 1, 1.75, and 2.75 [51].The WAD estimates will equal 3 × AD, 1.07 × AD, 1.6 × AD, and 0.52 when annual defoliation levels are 1 (AD: 0%-5%), 2 (AD: 6%-30%), 3 (AD: 31%-60%), and 4 (AD: 61%-100%), respectively.Additionally, both tree species (S) and ages (A) influence mortality during SBW outbreaks [52].For balsam fir forests, we used observed basal area data from 2010-2016 for the North Shore region of Québec to compare with model simulation results from the model (see Table A1) and then calibrated S = 0.4 and A as Equation ( 4): There are two mortality rates in the original TRIPLEX 1.0 model: natural mortality (NM; %/year) and competition mortality (CM; %/year).The CM rate is estimated based only on canopy competition for sunlight [28] and be given an initial value as 1.2%/year.In addition, we assume that DM rates greater than 3%/year means all dead trees are caused by severe defoliation, where the rates of NM and CM are zero (Table 3).Note that in this study we did not add new trees when defoliation caused mortality and dead trees were transferred to the litter pool in the next iteration.SBW preferentially eats flower cones, resulting in little to no seed production during an outbreak [53].

Biomass Loss Resulting from Mortality and Defoliation
The TRIPLEX-Insect model can simulate forest biomass loss by monthly time steps through the two following processes: (1) Mortality calculated in the Insect sub-model is sent to TRIPLEX 1.0 for NPP sub-model simulation.The DM is used to recalculate the survival and dead forest biomass, including branches, stems, and roots.Dead biomass becomes litter, which is then used in the simulation of soil C and N dynamics; (2) Leaf biomass loss is calculated based on AD (not mortality).When SBW-induced defoliation occurs, foliage is destroyed but not all foliage is consumed by insects; thus, we used a waste constant (W = 0.4) taken from Regniere et al. (1991) [49] to estimate supernumerary foliage litter.In other words, 40% of total C from defoliation is transferred directly to litter as litterfall, and 60% is distributed to insect lifecycles (e.g., insect biomass and insect respiration).

Compensatory Mechanisms of Foliage Loss
Typically, the physiological response to defoliation activates multiple compensatory mechanisms that increase total crown photosynthesis in evergreen species [54].In this study, we considered three compensatory methods when defoliation occurred, by (1) an increase in the photosynthetic efficiency of remaining leaves when December LAI was less than 20% of January LAI in the same year [54].For example, we recalculated alpha carbon (Cα) using Equation ( 5): (2) Defoliators reduced leaf biomass, resulting in more consumption of nonstructural carbohydrate (NSC), which was transferred to the growth of new leaves [55].Here, we assumed and estimated this following Equation ( 6): where NL is increment biomass of new leaves converted by NSC from other aspects of trees (i.e., stems, coarse roots and fine roots).F is the index of NL that is based on different CD conditions.For example, defoliation increases NSC consumption and reduces NSC production [56], we assumed NSC is exhausted when CD ≥ 4.5, so F = 0 in this case.The index of annual defoliation (IAD) was calculated by annual defoliation, and the biomass of new leaves scales the conversion from biomass of stems (85%), coarse roots (10%) and fine roots (5%).(3) A change in the nutritional distribution of foliage and roots [54,57], for example, we increased the leaf fraction (EtaF) from NPP using Equation ( 7), while we decreased the fine root fraction (EtaFR) from NPP using Equation ( 8): where EtaS is stem fraction from NPP and EtaCR is the coarse roots fraction from NPP.We selected data that included mean diameter at breast height (DBH), mean tree height, mean volume, and mean basal area from primarily balsam fir stands (greater than 85% balsam fir) in Cape Breton Island to calibrate the compensatory setting in our model (see Table A2).

Sensitivity Experiments
Using the TRIPLEX-Insect model, we conducted sensitivity experiments in order to test the effects of defoliation on C dynamics under different outbreaks durations, different initial forest age, and different defoliation severity in defoliated forests of balsam fir.Six different periods were estimated in which the duration of the outbreaks (10 or 20 years) and started stand ages (41, 101 or 161 years) varied given that younger stands are less vulnerable than older stands.Defoliation took 10 years in Short-Young (from 41 to 50 years tree age), Short-Mature (from 101 to 110 years stand age), and Short-Old (from 161 to 170 years stand age), while defoliation took 20 years in Long-Young (from 41 to 60 years stand age), Long-Mature (from 101 to 120 years stand age), and Long-Old (from 161 to 180 years stand age).We used five different defoliation scenarios representing a gradient from light to severe defoliation (no defoliation (control), 20% AD, 40% AD, 60% AD, and 80% AD) under the same AD severity level for the entire period of the outbreak.The TRIPLEX-Insect model was initialized using tree seedlings from the year of germination in all cases.

Model Validation
Although there are no direct data to verify C flux estimations made by the TRIPLEX-Insect model from 2006 to 2016 in the North Shore region of Québec, Canada, we surveyed 16 records of tree mortality in the QC1 (yearly record from 2014 to 2016), QC2 (yearly record from 2008 to 2016), and QC3 (yearly record from 2011 to 2014).We compared simulated mortality to field measurements to evaluate the accuracy of the TRIPLEX-Insect model.The model predictions fitted quite well the observed data (Figure 2A).In addition, we verified the compensatory mechanism in our model by using the current annual volume increment (CAI) data published by Ostaff and MacLean [35,36] from Cape Breton Island, Nova Scotia, Canada.Overall, TRIPLEX-Insect seems to underestimate the observed CAI by 10%-20% (Figure 2B). the same AD severity level for the entire period of the outbreak.The TRIPLEX-Insect model was initialized using tree seedlings from the year of germination in all cases.

Model Validation
Although there are no direct data to verify C flux estimations made by the TRIPLEX-Insect model from 2006 to 2016 in the North Shore region of Québec, Canada, we surveyed 16 records of tree mortality in the QC1 (yearly record from 2014 to 2016), QC2 (yearly record from 2008 to 2016), and QC3 (yearly record from 2011 to 2014).We compared simulated mortality to field measurements to evaluate the accuracy of the TRIPLEX-Insect model.The model predictions fitted quite well the observed data (Figure 2A).In addition, we verified the compensatory mechanism in our model by using the current annual volume increment (CAI) data published by Ostaff and MacLean [35,36] from Cape Breton Island, Nova Scotia, Canada.Overall, TRIPLEX-Insect seems to underestimate the observed CAI by 10%-20% (Figure 2B).

Sensitivity Analysis
Simulated outputs of GPP and Ra exhibited significant sensitivity to the level of defoliation.Five or 10 years after the beginning of an outbreak, we found a greater increase in GPP than Ra with 20% AD, while a greater loss in GPP compared to Ra with 60% AD and 80% AD (Table 4).Changes in GPP compared to no defoliation (control) increased from 1.32% to 4.25% while changes in Ra ranged from 0.74% to −1.64% as 40%AD levels after 5 years of an outbreak starting (Table 4).In addition, in all scenarios, with 60% and 80% AD, total GPP and Ra loss reached 100% because no trees survived.Indeed, medium and severe defoliation (40%-80% AD) led to the mortality of all trees when the SBW outbreak was programmed to last for 20 years.Therefore, we ended the simulation when all trees were dead.
The sensitivity experiments simulated negative NEP following defoliation for all scenarios evaluated (Figure 3A-F).Although we observed a pattern of rapid decline for all defoliation levels, more severe AD (from 20% to 80%) led to a greater decrease in annual NEP in balsam fir forests (Figure 3A-F).

Sensitivity Analysis
Simulated outputs of GPP and Ra exhibited significant sensitivity to the level of defoliation.Five or 10 years after the beginning of an outbreak, we found a greater increase in GPP than Ra with 20% AD, while a greater loss in GPP compared to Ra with 60% AD and 80% AD (Table 4).Changes in GPP compared to no defoliation (control) increased from 1.32% to 4.25% while changes in Ra ranged from 0.74% to −1.64% as 40%AD levels after 5 years of an outbreak starting (Table 4).In addition, in all scenarios, with 60% and 80% AD, total GPP and Ra loss reached 100% because no trees survived.Indeed, medium and severe defoliation (40%-80% AD) led to the mortality of all trees when the SBW outbreak was programmed to last for 20 years.Therefore, we ended the simulation when all trees were dead.
The sensitivity experiments simulated negative NEP following defoliation for all scenarios evaluated (Figure 3A-F).Although we observed a pattern of rapid decline for all defoliation levels, more severe AD (from 20% to 80%) led to a greater decrease in annual NEP in balsam fir forests (Figure 3A-F).We found that the duration of the outbreak and stand age also influence forest recovery.Following the end of the outbreak, younger (ages range from 40 to 60 years) stands (Figure 3A,D) recovered C sinks (NEP > 0) more rapidly and achieved higher levels of NEP compared to scenarios with older stands.Shorter duration of outbreaks resulted in a greater tree survival (20% AD and 40% AD in Figure 3A-C).For outbreaks lasting 20 years, only scenarios with 20% AD were able to recover (Figure 3D-F) because the other scenarios resulted in total tree mortality.

Effects of Defoliation on Carbon Fluxes and Stocks
The lowest GPP and NPP values from our model simulations occurred in 2014 for QC1 and QC3, and in 2015 for QC2 in the North Shore region of Québec, Canada.In comparing simulated results of insect disturbances to conditions of no defoliation (control), Figure 4A-C show that annual GPP had We found that the duration of the outbreak and stand age also influence forest recovery.Following the end of the outbreak, younger (ages range from 40 to 60 years) stands (Figure 3A,D) recovered C sinks (NEP > 0) more rapidly and achieved higher levels of NEP compared to scenarios with older stands.Shorter duration of outbreaks resulted in a greater tree survival (20% AD and 40% AD in Figure 3A-C).For outbreaks lasting 20 years, only scenarios with 20% AD were able to recover (Figure 3D-F) because the other scenarios resulted in total tree mortality.

Effects of Defoliation on Carbon Fluxes and Stocks
The lowest GPP and NPP values from our model simulations occurred in 2014 for QC1 and QC3, and in 2015 for QC2 in the North Shore region of Québec, Canada.In comparing simulated results of insect disturbances to conditions of no defoliation (control), Figure 4A-C show that annual GPP had decreased 98.23% (in 2014), 99.41% (in 2015), and 43.1% (in 2014), respectively, in the QC1, QC2, and QC3 plots.On the other hand, annual NPP decreased by 99.4% (Figure 4D), 99.99% (Figure 4E), and 44.44% (Figure 4F), which was slightly higher compared to GPP losses.NEP was sensitive to changes in the severity of AD and CD (Figure 4G-I).Two of the stands (QC1 and QC2) showed signs of weak recovery whereas QC3 did not.The simulation of NEP for all stands changed from positive to negative following the commencement of defoliation, which resulted in all plots becoming significant C sources after only a few years of SBW-induced defoliation.
Forests 2018, 9, x FOR PEER REVIEW 9 of 17 decreased 98.23% (in 2014), 99.41% (in 2015), and 43.1% (in 2014), respectively, in the QC1, QC2, and QC3 plots.On the other hand, annual NPP decreased by 99.4% (Figure 4D), 99.99% (Figure 4E), and 44.44% (Figure 4F), which was slightly higher compared to GPP losses.NEP was sensitive to changes in the severity of AD and CD (Figure 4G-I).Two of the stands (QC1 and QC2) showed signs of weak recovery whereas QC3 did not.The simulation of NEP for all stands changed from positive to negative following the commencement of defoliation, which resulted in all plots becoming significant C sources after only a few years of SBW-induced defoliation.When compared with control, SBW outbreaks simulations resulted in above (AGB) and below (BGB) ground decline in biomass for all plots (Figure 5A,B).In contrast, total litterfall (TL) stocks showed increasing trends for all simulations 2006 to 2016 (Figure 5D).The differences in litter fall were 28.09%, 252.56% and 7.97% higher than in scenarios with no defoliation for QC1, QC2 and QC3, respectively, in 2016.In addition, soil carbon (SC) stocks indicated almost no change for all plots compared to no defoliation (control) (Figure 5C).When compared with control, SBW outbreaks simulations resulted in above (AGB) and below (BGB) ground decline in biomass for all plots (Figure 5A,B).In contrast, total litterfall (TL) stocks showed increasing trends for all simulations 2006 to 2016 (Figure 5D).The differences in litter fall were 28.09%, 252.56% and 7.97% higher than in scenarios with no defoliation for QC1, QC2 and QC3, respectively, in 2016.In addition, soil carbon (SC) stocks indicated almost no change for all plots compared to no defoliation (control) (Figure 5C).

Prediction of Mortality and NEP with Different Defoliation Scenarios
Simulated results showed cumulative mortality was 24.4%, 88.1% and 24.6%, respectively, for QC1, QC2 and QC3 from 2006 to 2016.We used 30 years (1981-2010) of averages of monthly climate data and three different defoliation scenarios (no defoliation (control), 20% AD per year, and 45% AD per year) to forecast mortality from 2017 to 2025 in the North Shore region of Québec, Canada.In our predicted results, balsam fir mortality should reach 100% if continuous moderate defoliation (45% AD per year) occurs at QC1 until 2021, at QC2 until 2018, and at QC3 until 2022.This means that the risk for balsam fir stands will be significant if AD cannot be lowered.
In our predictions, both QC1 and QC2 showed increasing trends in NEP with no defoliation (control) and with 20% AD from 2017 to 2025 (Figure 6A,B).However, we also found that there was still a slight C source (NEP < 0) at both QC1 and QC2 between 2017-2025.Given its younger ages, a more rapid recovery was observed at QC3 compared to both QC1 and QC2 for all defoliation scenarios.Recovery in NEP reached 1.91 Mg C ha -1 year -1 (no defoliation; control) and 2.22 Mg C ha - 1 year -1 (20% AD) in 2020.Predicted showed C sink (NEP > 0) occurred again and kept with no defoliation (control) at QC3 from 2017 to 2025 (Figure 6C); however, QC3 quickly converted from a C sink (NEP > 0) to a source (NEP < 0) again with 20% AD.

Prediction of Mortality and NEP with Different Defoliation Scenarios
Simulated results showed cumulative mortality was 24.4%, 88.1% and 24.6%, respectively, for QC1, QC2 and QC3 from 2006 to 2016.We used 30 years (1981-2010) of averages of monthly climate data and three different defoliation scenarios (no defoliation (control), 20% AD per year, and 45% AD per year) to forecast mortality from 2017 to 2025 in the North Shore region of Québec, Canada.In our predicted results, balsam fir mortality should reach 100% if continuous moderate defoliation (45% AD per year) occurs at QC1 until 2021, at QC2 until 2018, and at QC3 until 2022.This means that the risk for balsam fir stands will be significant if AD cannot be lowered.
In our predictions, both QC1 and QC2 showed increasing trends in NEP with no defoliation (control) and with 20% AD from 2017 to 2025 (Figure 6A,B).However, we also found that there was still a slight C source (NEP < 0) at both QC1 and QC2 between 2017-2025.Given its younger ages, a more rapid recovery was observed at QC3 compared to both QC1 and QC2 for all defoliation scenarios.Recovery in NEP reached 1.91 Mg C ha −1 year −1 (no defoliation; control) and 2.22 Mg C ha −1 year −1 (20% AD) in 2020.Predicted showed C sink (NEP > 0) occurred again and kept with no defoliation (control) at QC3 from 2017 to 2025 (Figure 6C); however, QC3 quickly converted from a C sink (NEP > 0) to a source (NEP < 0) again with 20% AD.

Impact of Defoliation Intensity and Duration
Defoliated forests are predicted to lose C because photosynthesis is interrupted by the removal of foliage while respiration and decomposition continue [58,59].The model simulations showed that more severe defoliation and a longer outbreak duration would result in a higher mortality level, lower recovery capacity, and thus more C loss in balsam fir forest ecosystems.Gray [60] reported that more severe outbreaks and a longer duration would result young and rapid carbon-accumulating trees replacing old and slow carbon-accumulating trees.In our study, annual GPP and NPP recovered significantly only at QC3 in 2015 and 2016 (Figure 4C,F), and this was perhaps due to the fact that QC3 had both the younger age and lower CD than both QC1 and QC2.
In our sensitivity experiments, we found that when defoliation severity increases, both GPP and Ra decrease.However, we did not find any correlation between GPP or Ra dynamics and the various forest ages investigated under the same defoliation intensity (see Table 4).These results suggest that when defoliation-driven tree mortality is not excessively high, GPP's response to defoliation severity is greater than Ra's in both mature and immature balsam fir forest stands.The possible reason for this difference is that leaves are the only organ of photosynthesis rather than respiration, so defoliation should have a greater and direct impact on GPP than Ra.When AD is low, it tends to act as a thinning mechanism, such that NEP exceeds that found in non-defoliated stands.Virgin and MacLean [61] also found that subsequent forest structure depends on defoliation severity, although SBW outbreaks caused a "patchier" response compared to forestry (thinning) operations.In our simulations, C fluxes (i.e., GPP, NPP and NEP) had a rapid decrease when DM occurred in 2013 (at QC1 and QC3) and in 2010 (at QC2).However, the C fluxes were higher under outbreak conditions compared to the no disturbance scenario before DM occurred (Figure 4).
Medvigy et al. [24] reported that forest ecosystems respond strongly to C sequestration when defoliation duration varies from 5 to 15 years, but they exhibit a relatively weak response when defoliation duration varies by more than 15 years.In our study, simulated and predicted results also found that defoliation duration influences C losses and recovery in balsam fir forest stands.Longterm defoliation duration results in higher mortality and a lower recovery capacity in forests (Figure 3D,E,F).However, when the defoliation duration is less than 10 years, tree mortality is low because CD remains below lethal levels with the exception of the most severe defoliation cases.
Many studies have shown that insect damage could convert forests from C sinks to C sources [10,20,62], and this is consistent with our simulation and sensitivity experimental results.However, we observed a rapid recovery trend in annual NEP if defoliation is reduced (red line after 2015 in Figure 6A-C).In other words, the C source period may not protract following insect outbreaks if

Impact of Defoliation Intensity and Duration
Defoliated forests are predicted to lose C because photosynthesis is interrupted by the removal of foliage while respiration and decomposition continue [58,59].The model simulations showed that more severe defoliation and a longer outbreak duration would result in a higher mortality level, lower recovery capacity, and thus more C loss in balsam fir forest ecosystems.Gray [60] reported that more severe outbreaks and a longer duration would result young and rapid carbon-accumulating trees replacing old and slow carbon-accumulating trees.In our study, annual GPP and NPP recovered significantly only at QC3 in 2015 and 2016 (Figure 4C,F), and this was perhaps due to the fact that QC3 had both the younger age and lower CD than both QC1 and QC2.
In our sensitivity experiments, we found that when defoliation severity increases, both GPP and Ra decrease.However, we did not find any correlation between GPP or Ra dynamics and the various forest ages investigated under the same defoliation intensity (see Table 4).These results suggest that when defoliation-driven tree mortality is not excessively high, GPP's response to defoliation severity is greater than Ra's in both mature and immature balsam fir forest stands.The possible reason for this difference is that leaves are the only organ of photosynthesis rather than respiration, so defoliation should have a greater and direct impact on GPP than Ra.
When AD is low, it tends to act as a thinning mechanism, such that NEP exceeds that found in non-defoliated stands.Virgin and MacLean [61] also found that subsequent forest structure depends on defoliation severity, although SBW outbreaks caused a "patchier" response compared to forestry (thinning) operations.In our simulations, C fluxes (i.e., GPP, NPP and NEP) had a rapid decrease when DM occurred in 2013 (at QC1 and QC3) and in 2010 (at QC2).However, the C fluxes were higher under outbreak conditions compared to the no disturbance scenario before DM occurred (Figure 4).
Medvigy et al. [24] reported that forest ecosystems respond strongly to C sequestration when defoliation duration varies from 5 to 15 years, but they exhibit a relatively weak response when defoliation duration varies by more than 15 years.In our study, simulated and predicted results also found that defoliation duration influences C losses and recovery in balsam fir forest stands.Long-term defoliation duration results in higher mortality and a lower recovery capacity in forests (Figure 3D-F).However, when the defoliation duration is less than 10 years, tree mortality is low because CD remains below lethal levels with the exception of the most severe defoliation cases.
Many studies have shown that insect damage could convert forests from C sinks to C sources [10,20,62], and this is consistent with our simulation and sensitivity experimental results.However, we observed a rapid recovery trend in annual NEP if defoliation is reduced (red line after 2015 in Figure 6A-C).In other words, the C source period may not protract following insect outbreaks if defoliation intensity is not excessively high and the duration is relatively short.Albani et al. [23] predicted similar impacts of the hemlock wooly adelgid on eastern United States forests, i.e., an 8% decrease in NEP during the period 2000-2040 and a 12% increase without disturbance during the period 2040-2100.

Model Performance, Limitations, and Future Research
To improve the ability of TRIPLEX to model impacts of insect disturbances, we developed the new TRIPLEX-Insect model which is able to simulate and predict C dynamics and tree mortality caused by defoliation.This ability to test scenarios will be useful in planning SBW management.The TRIPLEX-Insect model is inherited from the original model (TRIPLEX 1.0) and expands its capabilities to simulate key variables related to SBW disturbances, such as changes in mean DBH, mean height, total stand volume, biomass, and C dynamics.The extensive application of this model, however, will first require further research on the five following aspects: (1) To accurately model the influence of defoliation on regional forest C dynamics, a comprehensive understanding of the spatial and temporal development of insect outbreaks is required [50,63].In our study, we only calibrated and validated the model at the stand level, yet spatial processes that drive outbreaks to cross thresholds and become catastrophic disturbances have not been completely simulated and are difficult to predict using the current version of the TRIPLEX-Insect model.In the future, we intend to link our model to remote sensing analyses to scale up the model from a stand level to a landscape level for regional-level applications.
(2) Another key element that will need to be included in our model is the influence of species composition.In the current application, we tested the model in balsam fir stands whereas the SBW has been observed to have a reduced effect in mixed-host hardwood stands [64][65][66].To more fully understand forest level C losses, we will also need to include other host species, such as white and black spruce.Even if these species are less vulnerable to mortality caused by the SBW [16], it has recently been shown that associational susceptibility occurs when spruce trees are mixed with fir trees, and as outbreak severity and duration length increase [37].
(3) Climate can also influence SBW outbreak characteristics, such as duration and intensity.For instance, seasonal temperatures and precipitation impact SBW life cycles (e.g., feeding and survival rates [49,67,68]) and population dynamics (e.g., dispersal area and develop rate [49,69,70]).To more accurately predict ecosystem response to SBW disturbances in the future, forthcoming modeling research should investigate the effects of future climate change on defoliation dynamics.
(4) Although our model is not capable of considering all SBW dynamics which impact ecosystem C budgets (e.g., insect respiration processes and decomposition of dead larvae), our method of simulating defoliation directly as an input parameter will be easily applicable to similar pest species from which defoliation data are regularly collected.Furthermore, since the TRIPLEX-Insect model can run with only an AD input, it should be possible to modify weights of monthly defoliation distribution and foliage waste for other defoliator species.
(5) TRIPLEX did not initially consider the concurrent simulation of tree mortality and regeneration; thus, the accuracy of the long-term simulations on C dynamics is limited.In this study, we only predicted short-term C dynamics (from 2017 to 2020) at sites in the North Shore region of Québec, Canada.This approach of only modeling C dynamics during the outbreak when establishment of new balsam fir was chosen so that we could avoid uncertainty, linked to the post-disturbance establishment of new trees, when testing our model.This could also cause some uncertainty related to model simulations due to simple assumptions concerning C pool totals and regeneration deficiencies.

Conclusions
There is currently a lack of accurate C budget quantification methods in conjunction with insect disturbances.Accordingly, we successfully developed the TRIPLEX-insect model, which is a new process-based sub-model that takes into account AD, CD, and tree mortality inputs to ascertain C dynamic response to balsam fir forests under SBW disturbances.This study found that the capacity of C sequestration is highly correlated with both defoliation intensity and duration in forest ecosystems.Furthermore, our sensitivity experiments also showed that tree age was a key factor which determined the capacity of recovery in the boreal forest following SBW outbreaks.Additionally, we found that GPP was more sensitive to defoliation intensity compared to Ra during outbreaks for both mature and immature forests.The simulated results showed that C fluxes (i.e., GPP, NPP and NEP) would increase to higher levels under conditions of insect disturbance than no disturbance (control) before DM occurred.The predicted results also indicated that the C source period would not be protracted if DM is not excessively high.Overall, our results demonstrated that the TRIPLEX-Insect model can be a useful tool for assisting forest insect management initiatives and provides more accurate C budget assessments of balsam fir forests under spruce budworm disturbances.

Figure 1 .
Figure 1.Framework and flow chart of TRIPLEX-Insect model (modified from Peng et al. (2002)[25]).Rectangles represent pools or state variables; ovals represent simulated processes; diamond-shaped boxes represent judgment process; dotted lines represent controls; and solid lines represent carbon (C), nitrogen (N), and water flow.

Figure 1 .
Figure 1.Framework and flow chart of TRIPLEX-Insect model (modified from Peng et al. (2002) [25]).Rectangles represent pools or state variables; ovals represent simulated processes; diamond-shaped boxes represent judgment process; dotted lines represent controls; and solid lines represent carbon (C), nitrogen (N), and water flow.
4 NM Normal mortality (%/year) 0.6 b 0 b CM Competition mortality (%/year) 1.2 b 0 b a Calibrated based on data from the Cape Breton Island study sites.b Denotes an assumed value.

Figure 4 .
Figure 4. Simulated annual GPP, NPP, and NEP under conditions of insect disturbance (red line) and no disturbance (control; green dotted line) by the TRIPLEX-Insect model using observed annual defoliation scenarios (red bar plot) from the QC1 (A,D,G,J), QC2 (B,E,H,K), and QC3 (C,F,I,L) plots from 2006 to 2016.

Figure 4 .
Figure 4. Simulated annual GPP, NPP, and NEP under conditions of insect disturbance (red line) and no disturbance (control; green dotted line) by the TRIPLEX-Insect model using observed annual defoliation scenarios (red bar plot) from the QC1 (A,D,G,J), QC2 (B,E,H,K), and QC3 (C,F,I,L) plots from 2006 to 2016.

Forests 2018, 9 , 17 Figure 6 .
Figure 6.Simulated NEP from 2010 to 2016 (on the left side of the vertical dotted line) and predicted NEP under three different defoliation levels (0% AD, 20% AD and 45% AD) from 2017 to 2025 (on the right side of the vertical gray dotted line) for the QC1 (A), QC2 (B), and QC3 (C) sites.Horizontal dotted line represent results for disturbance (control).

Figure 6 .
Figure 6.Simulated NEP from 2010 to 2016 (on the left side of the vertical dotted line) and predicted NEP under three different defoliation levels (0% AD, 20% AD and 45% AD) from 2017 to 2025 (on the right side of the vertical gray dotted line) for the QC1 (A), QC2 (B), and QC3 (C) sites.Horizontal dotted line represent results for disturbance (control).

Table 1 .
Stand location, soil type, and observed defoliation of the balsam fir study sites in the North Shore region of Québec, Canada.

Table 2 .
Monthly distribution of annual defoliation (AD) in May, June, and July based on the mean temperature in May.

Table 3 .
Key parameter values used in simulations of the TRIPLEX-Insect model.Parameters including Etas, EtaCR, EtaF, and EtaFR during endemic periods are based on the default values in TRIPLEX 1.0.For alpha carbon (Cα), EtaF, and EtaFR, severe defoliation was defined when LAI in December was less than 20% of LAI in January during the same year.For NM and CM, severe defoliation means DM greater than 0.03.

Table 4 .
Simulated total changes of GPP and Ra at 5 years and 10 years after the beginning of an outbreak compared to no defoliation (control) in the balsam fir forest.AD denotes annual defoliation.

Table 4 .
Simulated total changes of GPP and Ra at 5 years and 10 years after the beginning of an outbreak compared to no defoliation (control) in the balsam fir forest.AD denotes annual defoliation.