Modelling long-term impacts of mountain pine beetle outbreaks on merchantable biomass , ecosystem carbon , albedo , and radiative forcing

The ongoing major outbreak of mountain pine beetle (MPB) in forests of western North America has led to considerable research efforts. Yet many questions remain unaddressed regarding its long-term impacts, especially when accounting for the range of possible responses from the non-target vegetation (i.e., deciduous trees and lower-canopy shrubs and grasses). We used the Integrated BIosphere Simulator (IBIS) process-based ecosystem model along with the recently incorporated Marauding Insect Module (MIM) to quantify, over 240 years, the impacts of various MPB outbreak regimes on lodgepole pine 5 merchantable biomass, ecosystem carbon, surface albedo, and the net radiative forcing on global climate caused by the changes in ecosystem carbon and albedo. We performed simulations for three locations in British Columbia, Canada, having different climatic conditions, and four scenarios of various coexisting vegetation types with variable growth release responses. The impacts of MPB outbreaks on merchantable biomass (decrease) and surface albedo (increase) were similar across the 12 combinations of locations and vegetation coexistence scenarios. The impacts on ecosystem carbon and radiative forcing, however, 10 varied substantially in magnitude and sign depending upon the presence and response of the non-target vegetation, particularly for the two locations not subjected to growing-season soil moisture stress; this variability represents the main finding from our study. Despite major uncertainty in the value of the resulting radiative forcing, a simple analysis also suggested that the MPB outbreak in British Columbia will have a smaller impact on global temperature over the coming decades and centuries than a single month of global anthropogenic CO2 emissions from fossil fuel combustion and cement production. Moreover, we 15 found that: (1) outbreak severity (i.e., per-event mortality) had a stronger effect than outbreak return interval on the variables studied, (2) MPB-induced changes in carbon dynamics had a stronger effect than concurrent changes in albedo on net radiative forcing, and (3) the physical presence of MPB-killed dead standing trees was potentially beneficial to tree regrowth. Given that the variability of pre-outbreak vegetation characteristics can lead to very different regeneration pathways, the four vegetation coexistence scenarios we simulated probably only sampled the range of possible responses. 20

Abstract.The ongoing major outbreak of mountain pine beetle (MPB) in forests of western North America has led to considerable research efforts.However, many questions remain unaddressed regarding its long-term impacts, especially when accounting for the range of possible responses from the non-target vegetation (i.e., deciduous trees and lower-canopy shrubs and grasses).We used the Integrated BIosphere Simulator (IBIS) process-based ecosystem model along with the recently incorporated Marauding Insect Module (MIM) to quantify, over 240 years, the impacts of various MPB outbreak regimes on lodgepole pine merchantable biomass, ecosystem carbon, surface albedo, and the net radiative forcing on global climate caused by the changes in ecosystem carbon and albedo.We performed simulations for three locations in British Columbia, Canada, with different climatic conditions, and four scenarios of various coexisting vegetation types with variable growth release responses.The impacts of MPB outbreaks on merchantable biomass (decrease) and surface albedo (increase) were similar across the 12 combinations of locations and vegetation coexistence scenarios.The impacts on ecosystem carbon and radiative forcing, however, varied substantially in magnitude and sign, depending upon the presence and response of the non-target vegetation, particularly for the two locations not subjected to growing-season soil moisture stress; this variability represents the main finding from our study.Despite major uncertainty in the value of the resulting radiative forcing, a simple analysis also suggested that the MPB outbreak in British Columbia will have a smaller impact on global temperature over the coming decades and centuries than a single month of global anthropogenic CO 2 emissions from fossil fuel combustion and cement production.Moreover, we found that (1) outbreak severity (i.e., per-event mortality) had a stronger effect than outbreak return interval on the variables studied, (2) MPBinduced changes in carbon dynamics had a stronger effect than concurrent changes in albedo on net radiative forcing, and (3) the physical presence of MPB-killed dead standing trees was potentially beneficial to tree regrowth.Given that the variability of pre-outbreak vegetation characteristics can lead to very different regeneration pathways, the four vegetation coexistence scenarios we simulated probably only sampled the range of possible responses.

Introduction
The mountain pine beetle (MPB; Dendroctonus ponderosae Hopkins) is an insect native to forests of western North America, from northern Mexico to British Columbia, Canada (Safranyik and Carroll, 2006).Outbreaks of this bark beetle are characterized by high stand-level mortality of the target species, primarily lodgepole pine (Pinus contorta var.latifolia), but also other pines and, occasionally, other genera (NRCan, 2012).The MPB outbreak that started at the end of the previous century has reached an unprecedented level of documented severity, particularly in British Columbia where 18.1 Mha of forests have been affected (British Columbia, 2012b) and more than half of the merchantable pine volume has been killed (Walton, 2013).
Forests generally appear to recover well following MPB outbreaks (Axelson et al., 2009;Kashian et al., 2011;Hansen, 2014;Alfaro et al., 2015), which have recurred in western North America for thousands of years (Brunelle et al., 2008).However, forest managers face the decision of whether to proceed with salvage logging of MPB-killed dead standing trees (DSTs) and how best to do it (Griesbauer and Green, 2006;Bowler et al., 2012;Amoroso et al., 2013;Hawkins et al., 2013;Mathys et al., 2013;Landry and Ramankutty, 2015).MPB impacts also go beyond timber losses by modifying ecosystem carbon storage, thereby possibly affecting the ongoing climate change.The recent MPB outbreak has been estimated to decrease ecosystem carbon storage (cumulative values) by 270 Tg C between 2000 and 2020 in British Columbia (Kurz et al., 2008), by 580 Tg C between 1999 and 2050 in British Columbia (Arora et al., 2016), and by 15-26 Tg C between 2000 and2009 in the western United States (Ghimire et al., 2015), the last value increasing to 19-35 Tg C when including mortality from other bark beetles.
Recent reviews have identified various lingering knowledge gaps limiting the understanding of ecological and climatic effects caused by outbreaks of MPB and other forest insects (Liu et al., 2011;Seidl et al., 2011;Hicke et al., 2012a;Landry and Ramankutty, 2015).First, the mortality of many large trees often causes a growth release of the surviving non-target species and smaller host trees generally avoided by the MPB, which can alter the competition balance among the plant types present (Romme et al., 1986;Heath and Alfaro, 1990;Stone and Wolfe, 1996;Axelson et al., 2009;Amoroso et al., 2013;Hawkins et al., 2013;Hansen, 2014;Alfaro et al., 2015;Campbell and Antos, 2015).This growth release likely explains why field-based studies using the eddy covariance method have found the forest carbon balance to be more resilient than expected during an MPB outbreak or for close to a decade afterwards (Bowler et al., 2012;Brown et al., 2012;Reed et al., 2014).Therefore, modelling studies should allow for the possibility of a compensatory response from the surviving vegetation, including lower-canopy shrubs and grasses.Second, there is a need for more studies assessing the range of responses to different outbreak mortality levels and return intervals under the same background conditions, because comparisons performed across forest types and climates can be misleading.Third, the recurrence of MPB outbreaks calls for a long-term perspective going beyond a single mortality event.Fourth, the impact of the physical presence of MPB-killed DSTs on local exchanges of energy and water as well as the influence of these modified exchanges on carbon cycling have hitherto not been studied.Fifth, the climatic effects of MPB outbreaks are not limited to the carbon cycle, as the post-outbreak fall of DST needles and stems increases the reflection of incoming solar radiation, especially over seasonally snow-covered forests (Bright et al., 2013;Vanderhoof et al., 2014).In the only study to date aiming to quantify the net global climatic impact of MPB outbreaks, this albedo-induced cooling was estimated to be stronger than the warming from reduced ecosystem carbon storage (O'Halloran et al., 2012).
The main objective of our study was to use a modelling approach to evaluate the impacts of MPB outbreaks on four variables relevant to the forestry sector, land-atmosphere exchanges of carbon and energy, and global climate change, while explicitly addressing the five knowledge gaps identified above.Although more uncertain than empirical studies, process-based modelling approaches can provide a longerterm perspective on the impacts of MPB outbreaks and help assess interactions among several factors.Our purpose was not to forecast stand-level forest attributes (e.g., species-level basal area) but rather to contrast responses for very different scenarios about the presence and response of "non-target vegetation", which consisted of deciduous trees and lowercanopy shrubs and grasses that are never targeted by the MPB.Similarly, we did not account for all the factors affecting MPB population dynamics because we imposed idealized outbreak regimes, seeking here to provide initial insights on how impacts varied as a function of outbreak severity and return interval (e.g., Dietze and Matthes, 2014).

Overview
Our approach involved a set of different scenarios of coexistence between the MPB-targeted trees and non-target vegetation types; for each of these scenarios, we compared, over 240 years and for three locations in British Columbia, the impacts of various MPB outbreak regimes.In each instance, we included the explicit representation of interactions between MPB-killed DSTs and the carbon, energy, and water cycles.

Modelling the effects of MPB mortality
We used the recently developed Marauding Insect Module (MIM) incorporated within the Integrated BIosphere Simulator (IBIS) to simulate the effects of insect outbreaks.Here, we provide an overview of IBIS-MIM and refer readers to Landry et al. (2016) for more details.
The IBIS global ecosystem model was originally developed to estimate, within a single and consistent modelling framework, the land-atmosphere exchanges of carbon, energy, water, and momentum required by climate models, while simulating how vegetation phenology and spatial distribution respond to climate (Foley et al., 1996;Kucharik et al., 2000).IBIS represents coexisting upper (trees) and lower (shrubs and grasses) vegetation canopies as well as various soil and snow layers.The simulated exchanges of carbon, energy, and water depend upon the state of both canopies and the soil, including snow when present.Vegetation diversity is represented through various plant functional types (PFTs) characterized by different climatic constraints and parameters related to physiology, carbon dynamics, and energy exchanges.Photosynthesis and autotrophic respiration are typically computed on an hourly time step as a function of input climatic conditions.Competition balance and vegetation changes are determined at the end of each year based on the annual carbon balance of each PFT, except for the leaf phenology of deciduous PFTs, which is updated daily.For each PFT that can exist in the grid cell based on prevailing climatic conditions, leaf area index (LAI) cannot become lower than a very small, but non-zero, value; if a PFT undergoes 100 % mortality in a grid cell (e.g., as occurred with our Peak regime; see Sect.2.3), this "seed" LAI can therefore initiate regeneration.IBIS does not simulate establishment of many individuals for the same PFT in a grid cell.Annual litterfall is divided into daily transfers to soil, where carbon decomposition is modelled as a function of microbial biomass, soil temperature, and moisture.IBIS results compare relatively well with empirical data over large regions and several field sites, including in Canada (Foley et al., 1996;Delire and Foley, 1999;Kucharik et al., 2000;Lenters et al., 2000;El Maayar et al., 2001, 2002;Kucharik et al., 2006).
MIM was designed to simulate the effects of insect outbreaks within process-based ecosystem models similar to IBIS (Landry et al., 2016).MIM prescribes, at a daily time step, the direct insect-caused vegetation damage (i.e., defoliation and/or mortality), an approach that is similar to the "pathogen and insect pathways" from Dietze and Matthes (2014).The resulting impacts on vegetation dynamics and land-atmosphere exchanges of carbon, energy, and water are estimated by the supporting ecosystem model as a function of the post-outbreak state of the vegetation.MIM currently represents the effects of vegetation damage caused by outbreaks of three insect functional types (IFTs): broadleaf defoliators, needleleaf defoliators, and bark beetles.The bark beetle IFT used here was parameterized based on MPB-caused mortality of lodgepole pine.When a MPB outbreak occurs, mortality is assumed to begin on 1 August and increases linearly over 50 days (Landry et al., 2016) until reaching the userprescribed annual mortality level for the specific year (see Section 2.3).Killed trees become DSTs that interact with the exchanges of energy and water (e.g., absorbing solar radiation) but do not transpire or perform photosynthesis.DST carbon is gradually transferred to litter based on a predefined schedule: at year end for fine roots, over the 3 years following the year of mortality for needles, and, after a 5-year delay period, over 20 years for stem and coarse roots (Landry et al., 2016).IBIS then subdivides these annual amounts into daily transfers to soil.IBIS-MIM results after a simulated MPB outbreak generally agreed with previous studies (Landry et al., 2016).The only bias identified in that previous assessment of IBIS-MIM consisted of a lower snow depth/amount following MPB mortality vs. the no-outbreak control, contrary to the conclusion of most -but not all -previous studies.This bias likely resulted from overestimation by IBIS of heat storage within tree stems (Pollard and Thompson, 1995;El Maayar et al., 2001), including DSTs.IBIS-MIM might therefore underestimate the length of the snow cover season in MPB-attacked stands, thereby underestimating the consequent increases in springtime albedo and reflected solar radiation.This possible bias seems unlikely to be serious for the current study, because (1) areal snow coverage, which matters more for albedo than snow depth/amount, was the same for the outbreak and control results during most of the snow cover season; and (2) the generally earlier and faster snowmelt caused by MPB outbreaks (Pugh and Small, 2012;Mikkelson et al., 2013) is consistent with IBIS-MIM results.

Simulation design
IBIS requires input data for soil texture and climatic variables related to temperature, humidity (including precipitation and cloud cover), and wind speed.For all variables, including the ones provided in Table 1, we used the same input data as Landry et al. (2016): the mean 1961-1990 atmospheric CO 2 level, gridded 1961-1990 monthly mean data for climate (New et al., 1999), and survey data from versions 2.1 and 2.2 of the Soil Landscapes of Canada for soil (http://sis.agr.gc.ca/cansis/nsdb/slc/index.html).Note that contrary to gap models, IBIS does not have an intrinsic spatial resolution; however, as the computation of radiative forcing requires a specific area (see Appendix A), we used a nominal area of 1 ha here.Although we did not assess the effect of climate change, seeking to first understand the ecosystem responses within a stable climate regime, we considered the effect of varying climatic conditions by studying three locations in British Columbia, henceforth designated as northern, central, and southern (Fig. 1).These three locations have witnessed substantial MPB mortality since 2000 (Walton, 2013) under different climatic conditions (Table 1).The southern location is warmer than the central and northern locations.These last two locations have similar mean annual temperature, but the northern location has warmer summers and colder winters.Annual precipitation is similar in all locations, but summer rainfall is much lower in the southern location and results in drier conditions during the growing season.
We divided the simulations into four groups of different vegetation coexistence scenarios (Table 2).Five IBIS PFTs can coexist at the three locations considered here: the needleleaf evergreen (NE) trees targeted by MPB, broadleaf deciduous (BD) trees (e.g., trembling aspen; Populus tremuloides Michx.), and evergreen shrubs, cold-deciduous shrubs, and C 3 grasses in the lower canopy.The NEonly scenario allowed only NE trees and thus did not account for the possible response of the non-target vegetation.The NE-LC scenario allowed for the coexistence of NE trees and lowercanopy PFTs, but excluded BD trees.Note that in IBIS, coexisting PFTs compete for the same incoming solar radiation and soil water instead of being segregated into inde-Table 2. The four different scenarios simulated for the coexistence of plant functional types (PFTs).NE is needleleaf evergreen tree (i.e., the target PFT); LC is lower canopy (i.e., the sum of evergreen shrubs, cold-deciduous shrubs, and C 3 grasses); BD is broadleaf deciduous tree.

Scenario
PFTs allowed NEonly NE NE-LCcons NE and LC, with constant LC biomass from the first outbreak onwards NE-LC NE and LC AllPFT NE, LC, and BD pendent sub-grid tiles as in many similar models, so that tree PFTs actually shade the underlying lower canopy.The NE-LCcons scenario was similar to NE-LC except that the total biomass of the lower canopy was kept constant from the year of the first MPB outbreak onwards.Consequently, the lower canopy could increase its net primary productivity (NPP, in kg C m −2 yr −1 ) following MPB mortality (e.g., due to higher light penetration), but the additional productivity was transferred to litterfall so that the lower canopy was unable to grow and expand, thereby preventing further increases in productivity.The purpose of this constraint was to account for the effect of vegetation growth limitations not included in IBIS, such as nutrient availability.Finally, the AllPFT scenario allowed the five PFTs to freely compete throughout all years.Previous modelling studies quantifying MPB impacts on C eco or net ecosystem productivity (NEP, in kg C m −2 yr −1 ) (Kurz et al., 2008;Edburg et al., 2011;Ghimire et al., 2015;Arora et al., 2016) have resorted to NEonly-type approaches that did not account for the possible growth release of the non-target vegetation.We performed 12 sets (three locations and four coexistence scenarios) of five independent simulations; note that all simulations at a given location had the exact same weather.For a given set, the five independent simulations branched from the same 400-year spin up (which was sufficient for carbon pools to stabilize; see Landry et al. (2016) for more details) and consisted of (1) a no-outbreak control run; (2) a single 100 % MPB mortality event occurring during year 1 after the spin up, used to assess the theoretical maximum impacts and designated the Peak regime; (3) a regime of periodic single-year MPB outbreaks, also starting in year 1, with a mortality level of 16.6 % at return intervals of 20 years, designated the Small regime; (4) similar to Small, but with a single-year mortality level of 33.3 % at return intervals of 40 years, designated the Medium regime; and (5) similar to Small and Medium, but with a single-year mortality level of 50.0 % at return intervals of 60 years, designated the Large regime.We simulated a spatially homogeneous distribution of MPB-killed trees, based on the observation that MPB mortality was spatially more regular than the underlying distribution of trees in a 0.25 ha plot (Kashian et al., 2011) and to avoid the complications of introducing sub-grid spatial heterogeneity in IBIS.All simulations were continued for 240 years after the spin up, for a total of 640 years.Note that over these last 240 years, the mean mortality was 0.83 % yr −1 for the three periodic regimes (e.g., 16.6 % mortality every 20 years for Small), thereby allowing the effect of outbreak severity vs. return interval to be compared for the same mean mortality, but was only 0.416 % yr −1 for the Peak regime (i.e., a single 100 % mortality event over 240 years).We simulated single-year mortality events instead of manyyear outbreaks for two reasons.First, we wanted to focus on long-term results and a previous study with a model similar to IBIS found that, for the same level of total mortality occurring over 1, 5, or 15 years, differences in NEP became very small 25 years after the outbreaks (Edburg et al., 2011).Second, simulating a multiyear outbreak would raise the issues of its length and precise unfolding over consecutive years, introducing other complicating factors to our experimental setup that already considers the combination of three locations, four coexistence scenarios, and four outbreak regimes.

Response variables
We studied four variables: lodgepole pine merchantable biomass (B merch , in kg C m −2 ), ecosystem carbon (C eco , in kg C m −2 ), surface albedo (α, unitless), and radiative forcing (RF, in W m −2 ).B merch is highly relevant for the forestry sector, as it indicates the amount of lodgepole pine having commercial value.Note that, in the AllPFT scenario, BD trees (e.g., trembling aspen) could also have some commercial value, but we intentionally limited B merch to lodgepole pine due to the major commercial importance of this species in British Columbia.C eco goes beyond timber and accounts for all the carbon contained in live and dead pools including the soil, so that changes in C eco directly affect atmospheric CO 2 levels.Changes in α affect energy exchanges, with increases in α implying a cooling influence on the global climate.Finally, RF is used to assess the net impact of different perturbations on the global mean atmospheric surface temperature, without performing simulations with climate models (Myhre et al., 2013).In this study, the net RF (indicating warming if RF > 0, cooling if RF < 0) following MPB outbreaks was the sum of the radiative forcing from changes in atmospheric CO 2 and α.Unlike the other variables, RF is defined as a change between two different states; hence, RF results cannot be provided on a relative change (%) basis but must be absolute values for a given outbreak area.We stress that even if MPB impacts were simulated for a nominal area of 1 ha, the RF results we report are, by definition, for the net effect on global climate.We explain in Appendix A how we used IBIS output to compute B merch and RF.

Simplified analysis of maximum global climatic impact
Our simulation design was too simple to allow for a precise estimate of the global climatic impact from the MPB outbreak in British Columbia, but we used our RF results to bound the maximum value of this net warming or cooling impact.To do so, we identified among all our Peak simulations the two instances that led to the most extreme (positive and negative) annual RF values.We then recomputed these RF trajectories for an area of 18.1 Mha, which is the total area affected by the MPB outbreak (British Columbia, 2012b).Finally, we determined the value of a single pulse of actual (positive RF) or avoided (negative RF) fossil fuel CO 2 emissions that would invariably have, throughout the 240 years, a stronger radiative forcing than the MPB-caused RF (see Appendix A).This approach likely overestimated the maximum annual impact of MPB on the global climate for three reasons.First, not all area affected by MPB suffered 100 % mortality as prescribed in Peak.Second, the singleyear Peak mortality led to stronger extreme RF values compared to a gradual increase and decrease of the outbreak over many years.Third, the amount and composition of non-target vegetation are highly variable among MPB-attacked stands (Axelson et al., 2009;Amoroso et al., 2013;Hawkins et al., 2013;Pelz and Smith, 2013;Alfaro et al., 2015;Campbell and Antos, 2015) and this variability appears consequential for RF (see Sect. 3); hence the RF response from the same vegetation coexistence scenario was unlikely representative of the mean response over the entire area affected.

Comparison to previous studies
IBIS-MIM results following a MPB outbreak have already been found to agree with results from 38 field-, satellite-, and model-based studies for 28 variables related to vegetation dynamics and the exchanges of carbon, energy, and water, over daily to multiannual timescales (Landry et al., 2016).That previous assessment also provided time-sincedisturbance NPP results for different PFTs under the NE-LCcons scenario.Here we performed an additional assessment of IBIS-MIM, comparing its results for different variables related to the carbon cycle and for α to previous studies (see Tables S1-S2 in the Supplement).We restricted these comparisons to empirical studies that included control stands, except for the modelling study of Arora et al. (2016) that used a model similar to IBIS-MIM under a NEonlytype setting and provided results for one grid cell close to our central location.Given the differences in locations, cumulative mortality levels, and temporal patterns of mortality, IBIS-MIM agreed reasonably well with previous stud- ).The columns correspond to the three locations (Fig. 1) and the rows to the four vegetation coexistence scenarios (Table 2).Control values differed among locations and vegetation coexistence scenarios, so the same relative change (in %) across the 12 panels does not correspond to the same absolute change.
ies, providing a measure of confidence in the realism of the results shown below.

Transient results
MPB-caused reductions in B merch were similar across the three locations and four vegetation coexistence scenarios (Fig. 2).For the Peak regime, the single 100 % mortality event removed all B merch , after which 20-50 years were needed before NE trees became big enough to have any commercial value.The slower recovery of B merch in AllPFT compared to other scenarios resulted from the growth release of BD trees, which were able to grow strongly for a few decades but were very poor long-term competitors at all locations.For the three periodic regimes, the recurring MPB outbreaks prevented B merch from recovering to the control value as in Peak.We found some evidence of a biophysical interaction between DSTs and regrowing NE trees for the Peak regime in the NEonly, NE-LCcons, and NE-LC scenarios: after the 100 % mortality event, NPP of NE trees increased rapidly but, after about 20-25 years, decreased no-ticeably for ∼10 years before resuming again (similarly to the results shown in Fig. 2 from Landry et al., 2016).We believe this response resulted from the interception of radiation by DSTs, which warmed the upper canopy and initially allowed regrowing NE trees to perform photosynthesis more efficiently at a higher temperature.This warming effect decreased as DSTs fell, causing productivity of the regrowing NE trees to decline even though they were intercepting more light.In the AllPFT scenario, the slower regrowth of NE trees caused them to miss the warming effect due to the presence of DSTs.
The impacts of MPB on C eco (Fig. 3) were much more complex than on B merch .In the NEonly scenario, the three periodic regimes led to gradual declines in C eco that showed indications of possible stabilization towards the end of the simulations, whereas for the Peak regime, C eco was still recovering after 240 years.The results were qualitatively the same in the NE-LCcons scenario, albeit with much smaller C eco because the lower-canopy growth release partially compensated for the death of NE trees.At the southern location, the ).The columns correspond to the three locations (Fig. 1) and the rows to the four vegetation coexistence scenarios (Table 2); the y axis scale differs across the four rows.Control values differed among locations and vegetation coexistence scenarios, so the same relative change (in %) across the 12 panels does not correspond to the same absolute change.
growing-season soil moisture stress probably explains why the growth release of non-target PFTs was only marginally stronger in the unconstrained NE-LC and AllPFT scenarios vs. NE-LCcons.Conversely, MPB outbreaks substantially increased total NPP at the northern and central locations for NE-LC and AllPFT, by inducing a strong growth release of the non-target vegetation and fostering the increased coexistence of PFTs occupying different ecological niches (upper vs. lower canopies, and evergreen needleleaf vs. deciduous broadleaf strategies) compared to undisturbed forests dominated by lodgepole pine.Here, the higher productivity of the non-target vegetation exceeded the productivity losses and gradual decomposition of killed NE trees; hence, after a delay of a few years to a few decades, C eco switched from negative to positive (Fig. 3g, h, j, and k).MPB outbreaks increased α for all locations and vegetation coexistence scenarios (Fig. 4).Results were very similar across locations and scenarios, except for smaller α increases in AllPFT (panels j, k, and l) for the Peak regime due to the absorption of radiation by BD trees (even when leafless dur-ing winter, as IBIS accounts for the snow-masking effect of stems) following their growth release.
RF varied substantially across locations and vegetation coexistence scenarios (Fig. 5; also see Fig. S1 in the Supplement for the CO 2 -based and α-based components of the total RF for Peak).For NEonly and NE-LCcons, and for the other scenarios at the southern location only, RF was almost always positive, indicating a warming effect of MPB outbreaks on the global climate.The small negative RF values observed for the Peak regime in three instances (panels d, e, and l) came from a combination of still slightly positive α, almost complete return of C eco to the control values (see Fig. 3), and long time lags in atmosphere-ocean CO 2 exchanges following vegetation regrowth (see Landry and Matthews (2016) for such time lags after the simulation of fire disturbances in a coupled climate-carbon model; here, these time lags were accounted for by the convolution approach used to compute RF as explained in Appendix A).At the northern and central locations, MPB outbreaks under the NE-LC and AllPFT scenarios led instead to a net cooling effect on the global www.biogeosciences.net/13/5277/2016/Biogeosciences, 13, 5277-5295, 2016 ).The columns correspond to the three locations (Fig. 1) and the rows to the four vegetation coexistence scenarios (Table 2).Control values differed among locations and vegetation coexistence scenarios, so the same relative change (in %) across the 12 panels does not correspond to the same absolute change.
climate, even though RF values were initially positive for at least 4 years (panels g, h, j, and k).

Mean effect
Figure 6 shows the mean time-averaged effect of MPB outbreaks on the four variables over the 240 years following the first outbreak.In the NEonly scenario, the results were almost equal in all three locations for a given MPB regime despite the different climatic conditions.When other PFTs coexisted with the target NE trees, however, the influence of climate became noticeable, especially in the NE-LC and AllPFT scenarios.Figure 6 also reveals that a higher perevent mortality generally caused stronger absolute effects than a shorter return interval.Indeed, for a given location and vegetation coexistence scenario, the departure from zero for the periodic outbreak regimes was generally in the following order: Large > Medium > Small.Moreover, the single 100 % mortality event under the Peak regime had a mean effect comparable to, and in many instances greater than, the mean effect under the three periodic regimes, despite a 240-year mean mortality of 0.416 % yr −1 for Peak vs. 0.83 % yr −1 for the periodic regimes.Figure 6 also suggests that RF was more closely linked to C eco than to α, which is supported by the transient results where RF basically mirrored C eco (compare Figs. 3  and 5; also see Fig. S1).Even though the α-caused cooling effect offset a fraction of the CO 2 -based warming when C eco was negative or added to the CO 2 -based cooling when C eco was positive, the sign of RF was primarily related to C eco changes.Spearman's rank correlation coefficient between C eco and RF across the 48 outbreak simulations presented in Fig. 6 was −0.99, indicating that greater decreases (increases) in C eco were almost invariably associated with greater increases (decreases) in RF, thereby leading to a greater warming (cooling) effect.On the contrary, Spearman's rank correlation coefficient between α and RF was only +0.21, indicating a weak association that was opposite to the actual effect of increased α on RF. , for 1 ha outbreaks) compared with the no-outbreak control (first outbreak occurred in 1).The columns correspond to the three locations (Fig. 1) and the rows to the four vegetation coexistence scenarios (Table 2); the y axis scale differs across the four rows.

Comparison of MPB to anthropogenic CO 2 emissions
The highest (positive) yearly RF value came from the NEonly scenario at the central location (Fig. 5b), whereas the lowest (negative) yearly RF value came from the NE-LC scenario at the northern location (Fig. 5g).We therefore recomputed these two RF trajectories over an area of 18.1 Mha to obtain the bounding positive and negative responses, respectively.As shown in Fig. 7, these bounding RF responses invariably had smaller absolute impacts than a single pulse of +0.83 Pg C of fossil fuel CO 2 for the warming case or a single pulse of −0.80 Pg C (avoided emissions) for the cooling case.

Influence of the non-target vegetation
Our main finding is that the non-target vegetation has a major influence on ecosystem carbon storage and net climatic impact following MPB outbreaks.The substantial variability in C eco (Fig. 3) and RF (Fig. 5) responses occurred despite almost identical recovery of the MPB-targeted dominant tree species across the four vegetation coexistence scenarios and three locations (Fig. 2).Previous studies not accounting for growth release of the non-target vegetation, including shrubs and grasses, may therefore have overestimated the MPB-caused decreases of C eco .Strong compensatory responses following MPB or other bark beetle outbreaks have also been reported in previous studies that used controls or considered several mortality levels (Heath and Alfaro, 1990;Stone and Wolfe, 1996;Klutsch et al., 2009;Griffin et al., 2011;Amoroso et al., 2013;Hawkins et al., 2013), including cases of NPP or aboveground tree biomass being higher with than without insect outbreaks for some period of time (Romme et al., 1986;Seidl et al., 2008;Pfeifer et al., 2011;Hansen, 2014).Modelling studies on a centennial timescale also found that C eco can be higher when insect disturbances are simulated than without them (Seidl et al., 2008;Albani et al., 2010) Figure 6.Mean effect over 240 years of the different MPB outbreak regimes on lodgepole pine merchantable biomass (B merch ), ecosystem carbon (C eco ), surface albedo (α), and radiative forcing (RF; in pico-W m −2 , for 1 ha outbreaks) compared with the no-outbreak control, for the three locations (Fig. 1) and four vegetation coexistence scenarios (Table 2).Control values differed among locations and vegetation coexistence scenarios, so the same relative change (in %; panels a through i) does not correspond to the same absolute change.
tion can contribute to a fast recovery of NEP following MPB outbreaks (Bowler et al., 2012;Brown et al., 2012) and even to stable growing-season NEP despite ongoing increases in MPB mortality (Reed et al., 2014).Another field-based study found C eco to be almost equal in control stands and stands affected 25-30 years earlier by a ∼25 %-mortality MPB outbreak (Kashian et al., 2013).We do not believe that one of the vegetation coexistence scenarios we simulated is fundamentally more realistic than the others.Rather, we believe that these four scenarios sample the ensemble of possible responses to MPB outbreaks, because the amount and composition of non-target vegetation vary substantially among MPB-attacked stands, even over short distances, which leads to variable regeneration pathways (Axelson et al., 2009;Amoroso et al., 2013;Hawkins et al., 2013;Pelz and Smith, 2013;Alfaro et al., 2015;Campbell and Antos, 2015).Given this substantial variability, trying to bracket the range of possible vegetation responses, like we did here, appears safer for large-scale modelling studies than aiming to forecast one specific course of events.Contrary to our simulation protocol, however, it is unlikely that the composition of non-target vegetation would remain unchanged from one outbreak to the next.Moreover, most forests in western North America undergo recurrent standclearing fires or wood harvests, resetting stands on trajectories that, according to the evidence currently available, are not strongly affected by the previous occurrence of MPB outbreaks (Hicke et al., 2012b;Harvey et al., 2014).

Climatic impact
The only previous study estimating the net impact of MPB outbreaks on global climate found a negative RF throughout the 14-year period studied due to an α-based cooling that was invariably stronger than the C eco -based warming (O'Halloran et al., 2012).This outcome is at odds with IBIS-MIM results, because our net RF depended critically upon the sign of the C eco change and, even for instances of overall net cooling, our RF values were positive during the first four years at least (see transient results in panels g, h, j, and k of Fig. S1 for instances of α-based cooling being temporarily stronger than C eco -based warming under Peak).This discrepancy likely involves methodological differences between the two studies but might also come from a spa-tiotemporal mismatch that could have affected the RF results from O' Halloran et al. (2012): their α was representative of MPB-killed stands (coming from time-since-mortality comparisons at the stand level), whereas their C eco was based on the regional-level results of Kurz et al. (2008), which were for stands killed at different times averaged with unaffected stands.
Our estimate of the global climatic impact due to the MPB outbreak in British Columbia (Fig. 7), although very simple, seems appropriate to bound the range of possible values.For the warming case, the maximum decrease in C eco (based on Peak from Fig. 3b, over 18.1 Mha) was equal to 818 Tg C ∼ 50 years after mortality; for this same case, the decrease 21 years after mortality was 490 Tg C. By comparison, Kurz et al. (2008) 2016) simulated a 1999-2050 decrease of 580 Tg C with a vegetation coexistence scenario similar to NEonly in an IBIS-like model.For the cooling case, the lower-canopy growth release in the unconstrained NE-LC scenario was likely too strong, causing the increase in C eco to be overestimated.Consequently, the actual impact probably lies within these bounding responses, which have smaller absolute RF than a single pulse of +0.83 Pg C (warming case) or −0.80 Pg C (cooling case) of fossil fuel CO 2 .Pulses of such magnitude represent approximately 1 month of current global CO 2 emissions from fossil fuel combustion and cement production (Le Quéré et al., 2015).Even though these results suggest a marginal impact on global temperature, the current MPB outbreak in British Columbia could add or offset a sizable share of the warming due to greenhouse gas (GHG) emissions from the province alone or Canada as a whole.Total GHG emissions were estimated at 61.5 Mt CO 2 eq in British Columbia for 2012 (British Columbia, 2012a) and at 726 Mt CO 2 eq in Canada for 2013 (Environment Canada, 2014), which are roughly equivalent to 0.017 and 0.20 Pg C of fossil fuel CO 2 , respectively.Therefore, the upper bound on the maximum global climatic impact of the current MPB outbreak in British Columbia, for either warming (+0.83 Pg C) or cooling (−0.80 Pg C), is equivalent to GHG emissions over roughly 50 years for British Columbia and 4 years for Canada.A more adequate assessment would require a dedicated study going beyond the simplified analysis presented here, which could only provide an upper bound for either a net warming or net cooling effect.
It is important to remember that the RF concept does not apply to changes in local temperature.Since the α-based cooling is local whereas the CO 2 -based effect is global, one could expect that MPB outbreaks always decrease local temperature.However, this perspective neglects the post-MPB changes in sensible and latent heat fluxes that also modulate local temperature.In summer, high levels of MPB mortal- ity have been found to increase surface temperature by up to a few • C due to reduced evapotranspiration (Griffin et al., 2011;Bright et al., 2013;Maness et al., 2013), a response Landry et al. (2016) obtained in their detailed assessment of IBIS-MIM.Landry et al. (2016) also found that MPB outbreaks decreased surface temperature in winter, but they could not find any empirical observations on this variable.

Management and research implications
Management activities aiming to prevent or respond to MPB outbreaks must consider several factors, including economic impacts (e.g., Patriquin et al., 2005) and potential effects on fire behaviour (see the review from Hicke et al., 2012b).Although we did not account for these factors or model management activities explicitly, our study suggests the following.First, our results are in line with an emerging body of empirical literature pointing towards the resilience of carbon storage in MPB-affected forests due to the growth release of the surviving vegetation (Bowler et al., 2012;Brown et al., 2012;Reed et al., 2014).Similar to previous studies (Stone and Wolfe, 1996;Klutsch et al., 2009;Griffin et al., 2011;Bowler et al., 2012;Vanderhoof et al., 2014), we found a growth release of shrubs and grasses -and not only trees -that could contribute to ecosystem-level resilience in carbon storage.This benefit could be compromised if non-target vegetation is damaged during salvage logging operations, which might therefore be detrimental to carbon stewardship.
We found indications of the potential growth release of surviving target trees in (1) the fast NPP recovery of NE trees for the three periodic regimes and (2) the comparable B merch mean decrease for the three periodic regimes vs. Peak despite www.biogeosciences.net/13/5277/2016/Biogeosciences, 13, 5277-5295, 2016 higher average mortality rates (Fig. 6).These outcomes are consistent with the relative stability of aboveground wood NPP for mortality levels below ∼ 60 % in a forest-level manipulation experiment (Stuart-Haëntjens et al., 2015).Second, our Peak results indicate that high amounts of DSTs could facilitate the growth of the surviving trees by warming the surrounding air.Salvage logging would therefore prevent this beneficial impact (in cool growing seasons) and dampen the growth release of surviving trees.The warming effect of DSTs appears coherent with empirical evidence on vegetation-temperature interactions in northern latitudes (Liu et al., 2006;Lee et al., 2011) but is likely smaller than simulated here due to the overestimation by IBIS of heat storage within tree stems.The impact of DSTs on the amount and partitioning (i.e., direct vs. diffuse) of solar radiation absorbed by the surviving vegetation, which all modulate NPP, is also probably more complex than simulated by IBIS-MIM.Field studies comparing the vertical profiles of temperature, total solar radiation, and diffuse solar radiation of stands with and without DSTs would be useful to resolve such questions.
Third, since MPB outbreaks do not necessarily warm the global climate, outbreak-preventing activities like preemptive logging might not mitigate climate change.Assessing the net climatic impact of salvage logging is even more complicated, because the exercise needs to go beyond comparing RF for salvaged and unsalvaged stands and account for the landscape-level redistribution of harvesting activities as well as the differences in the production and fate of the salvagedderived wood products compared to the no-salvage baseline (Lemprière et al., 2013;Landry and Ramankutty, 2015).
Fourth, IBIS-MIM transient results suggest that MPB impacts vary substantially not only across space but also through time.While the spatial variability of vegetation responses (Griesbauer and Green, 2006;Axelson et al., 2009;Amoroso et al., 2013;Hawkins et al., 2013;Pelz and Smith, 2013;Alfaro et al., 2015;Campbell and Antos, 2015) underlines the need for studies in many stands, the temporal variability reported here calls for continued or periodic data gathering over decades at the same sites.A better appreciation of the long-term effects of MPB outbreaks should foster adequate management responses in affected forests within the insect's native range as well as in the Canadian boreal forest where the MPB has already established and may spread as the regional climate warms (Cullingham et al., 2011;Nealis and Cooke, 2014).

Study limitations
Eddy covariance measurements of ecosystem-level carbon exchanges following MPB outbreaks (Bowler et al., 2012;Brown et al., 2012;Reed et al., 2014) extend for less than a decade and lack formal controls, hence limiting our capacity to validate long-term IBIS-MIM output.Many dendrochronological studies have gathered data over longer time periods, sometimes with formal controls, but the results do not translate directly into changes in stand-level NPP.Romme et al. (1986) estimated the effect of MPB on aboveground tree growth, but only for six stands 10-20 years after an outbreak.A single study by Kashian et al. (2013) quantified the impacts of MPB on C eco , through allometric equations and soil samples obtained at 12 stands 25-30 years after mortality; however, this study dealt primarily with fire and the MPB results were qualified as "preliminary" by the authors.Other empirical and modelling studies accounting for possible growth releases support the plausibility of IBIS-MIM responses when simulating additional PFTs besides NE trees, yet our results must be considered tentative because process-based modelling of vegetation competition (Kucharik et al., 2006;Moorcroft, 2006;Purves and Pacala, 2008) and non-stand-replacing disturbances (Bond-Lamberty et al., 2015) remains challenging.
We could not capture the effect of nitrogen availability on post-MPB vegetation recovery (Edburg et al., 2011) because IBIS does not simulate nutrient cycling.Another limitation of IBIS is the representation of a single NE tree PFT and a single BD tree PFT, whereas tree species other than lodgepole pine and trembling aspen often coexist in MPBattacked forests (Griesbauer and Green, 2006;Axelson et al., 2009;Amoroso et al., 2013;Hawkins et al., 2013;Campbell and Antos, 2015).Accounting for these other species could increase the range of post-MPB responses for C eco and RF and also partly offset the B merch reductions simulated here, which included lodgepole pine only.Since IBIS does not represent different age cohorts within the same PFT, we could not formally assess the growth release of individual surviving NE trees.The post-MPB response of younger trees would likely differ from those of mature trees and MPB impacts on tree demographics could further complicate standlevel responses following an outbreak, for example by increasing total biomass despite reduced productivity because of a strong decrease in competition-related mortality (Pfeifer et al., 2011).IBIS parameters for the different PFTs were also not specifically based on data gathered from British Columbian forests.
Finally, the overestimated heat storage in tree stems could underestimate the MPB-caused α increase, thereby biasing RF towards a warming effect.Conversely, another potential bias could overestimate the cooling effect of α increase: while a seminal study concluded that α and atmospheric CO 2 act with the same "efficacy" on the global surface temperature (Hansen et al., 2005), other studies have found that α has a lower efficacy than CO 2 (Hansen et al., 1997;Davin et al., 2007).MPB outbreaks can also affect the global climate through other mechanisms that are poorly constrained and that we did not consider, for example changes in atmospheric chemistry (Arneth and Niinemets, 2010) or in the partitioning between latent and sensible heat fluxes (Ban-Weiss et al., 2011).

Conclusions
Despite major progress over the last decades, various knowledge gaps still limit the understanding of the consequences of MPB outbreaks.In this study, we used a climate-driven process-based ecosystem model to estimate the long-term impacts of prescribed MPB outbreaks.We found that the differences in vegetation coexistence scenario and location had little influence over MPB impacts on lodgepole pine merchantable biomass (B merch ; Fig. 2) and surface albedo (α; Fig. 4).On the contrary, accounting for the non-target vegetation invariably reduced losses of ecosystem carbon (C eco ) and, at the two locations not subjected to growing-season soil moisture stress, even led to post-outbreak C eco increases for the two vegetation coexistence scenarios with the strongest growth release (Fig. 3).Although MPB-induced increases in α always had a cooling influence, the net global warming or cooling effect of MPB outbreaks was determined by the stronger carbon-based responses (see Fig. 6, and compare Fig. 5 with Fig. 3).A simple analysis suggested that the MPB outbreak in British Columbia will have less influence on global temperature over the coming centuries than 1 month of global anthropogenic CO 2 emissions at the 2014 level (Fig. 7).
The management and research implications of our study are fourfold.First, salvage logging, particularly when performed as clear-cut harvesting, may be detrimental to carbon stewardship when surviving trees and other lower-canopy vegetation are abundant.Second, salvage logging could slow forest recovery if, following high MPB mortality, tree productivity is indeed increased due to the physical presence of dead standing trees, a hypothesis that should be subject to empirical studies.Third, MPB disturbances might not necessarily lead to global warming, so activities aiming to prevent or control outbreaks (e.g., preemptive logging) should not be heralded as climate mitigation strategies without more detailed analyses.Fourth, the substantial spatiotemporal variability in MPB-induced changes suggests a need to support field studies that encompass a wide range of stand conditions and are maintained over several decades.

Data availability
Simulation results (as NetCDF files; http://www.unidata.ucar.edu/software/netcdf/) are available upon request from the corresponding author.

Figure 1 .
Figure 1.Three locations studied; the province of British Columbia is shaded.

Figure 2 .
Figure2.Transient effect of the different MPB outbreak regimes on lodgepole pine merchantable biomass (B merch ) compared with the no-outbreak control run (first outbreak occurred in year 1).The columns correspond to the three locations (Fig.1) and the rows to the four vegetation coexistence scenarios (Table2).Control values differed among locations and vegetation coexistence scenarios, so the same relative change (in %) across the 12 panels does not correspond to the same absolute change.

Figure 3 .
Figure3.Transient effect of the different MPB outbreak regimes on ecosystem carbon (C eco ) compared with the no-outbreak control run (first outbreak occurred in year 1).The columns correspond to the three locations (Fig.1) and the rows to the four vegetation coexistence scenarios (Table2); the y axis scale differs across the four rows.Control values differed among locations and vegetation coexistence scenarios, so the same relative change (in %) across the 12 panels does not correspond to the same absolute change.

Figure 4 .
Figure4.Transient effect of the different MPB outbreak regimes on surface albedo (α) compared with the no-outbreak control run (first outbreak occurred in year 1).The columns correspond to the three locations (Fig.1) and the rows to the four vegetation coexistence scenarios (Table2).Control values differed among locations and vegetation coexistence scenarios, so the same relative change (in %) across the 12 panels does not correspond to the same absolute change.

Figure 5 .
Figure5.Transient effect of the different MPB outbreak regimes on radiative forcing (RF; in pico-W m −2 , for 1 ha outbreaks) compared with the no-outbreak control (first outbreak occurred in 1).The columns correspond to the three locations (Fig.1) and the rows to the four vegetation coexistence scenarios (Table2); the y axis scale differs across the four rows.
simulated a 2000-2020 decrease of 270 Tg C with an inventory-based model omitting the possible growth release of the non-target vegetation and projections of MPB-caused mortality almost 40 % higher than more recent estimates (Walton, 2013), whereas Arora et al. (

Figure 7 .
Figure7.Comparison of the strongest warming and cooling radiative forcing (RF) responses from the MPB Peak outbreaks with the RF from a pulse of fossil fuel (FF) CO 2 emissions (in milli-W m −2 ).The MPB RF were computed for an outbreak area of 18.1 Mha; the warming response came from the NEonly scenario at the central location and the cooling response from the NE-LC scenario at the northern location.

Table 1 .
Input climate data and soil texture for the three locations.