Implementation of a Marauding Insect Module (MIM, version 1.0) in the Integrated BIosphere Simulator (IBIS, version 2.6b4) Dynamic Vegetation–Land Surface Model

1Department of Geography, McGill University, Montréal, Canada 2Natural Resources Canada, Canadian Forest Service, Northern Forestry Centre, Edmonton, Canada 3Liu Institute for Global Issues and Institute for Resources, Environment, and Sustainability, University of British Columbia, Vancouver, Canada 4Earth and Environmental Sciences and Biology, Irving K. Barber School of Arts and Sciences, University of British Columbia, Kelowna, Canada 5Department of Geography, Planning and Environment, Concordia University, Montréal, Canada 6Currently at the Department of Geography, Planning and Environment, Concordia University, Montréal, Canada


Introduction
The damage to plants caused by insects, particularly during outbreaks defined by sudden and major changes in insect population, are pervasive in terrestrial ecosystems and affect 25 not only vegetation dynamics, but also carbon, nutrient, energy, and water exchanges, and 2 even atmospheric chemistry (Landsberg and Ohmart, 1989;Hunter, 2001;Lovett et al., 2002;Kurz et al., 2008;Amiro et al., 2010;Arneth and Niinemets, 2010;Clark et al., 2010Clark et al., , 2012Stinson et al., 2011;Bowler et al., 2012;Brown et al., 2012Brown et al., , 2014Edburg et al., 2012;Hicke et al., 2012;Yang, 2012;Bright et al., 2013;Maness et al., 2013;Mikkelson et al., 2013a;Pugh and Gordon, 2013;Metcalfe et al., 2014;Reed et al., 2014;Seidl et al., 2014;5 Turcotte et al., 2014;Vanderhoof et al., 2014;Landry and Parrott, 2016). Yet the simulation of insect-induced plant damage in climate models has lagged behind the simulation of fire, even though the two disturbance types were recognized as climate-related phenomena worthwhile of explicit representation in Dynamic Global Vegetation Models (DGVMs) more than 15 years ago (Fosberg et al., 1999). 10 Since the term "DGVM" is often used for interactive vegetation models that estimate only some of the exchanges of carbon, energy, water, and momentum with the atmosphere (Prentice et al., 2007;Quillet et al., 2010), we will refer here to the subset of DGVMs that compute all required land-atmosphere exchanges while accounting for dynamic vegetation as Dynamic Vegetation-Land Surface Models (DVLSMs) to prevent possible confusion 15 (e.g., many DGVMs do not compute the land-to-atmosphere fluxes of shortwave and longwave radiation). Insect damage has been represented in DVLSMs in a handful of cases. Based on the empirical relationships of McNaughton et al. (1989), the ORganizing Carbon and Hydrology in Dynamic EcosystEms (ORCHIDEE) DVLSM accounts for background leaf consumption by herbivores (not limited to insects), but the realism of the resulting impact 20 on simulated tree mortality has been questioned by the authors themselves (Krinner et al., 2005). The effects of prescribed mortality due to mountain pine beetle (MPB; Dendroctonus ponderosae Hopkins) outbreaks in western US on coupled carbon-nitrogen dynamics  and water and energy exchanges  have been studied in the Community Land Model (CLM) DVLSM. Medvigy et al. (2012) used 25 the Ecosystem Demography version 2 (ED2) DVLSM to simulate the impacts of defoliation by the gypsy moth (Lymantria dispar Linnaeus) on vegetation coexistence and carbon dynamics in the eastern US. Background herbivory or insect outbreaks have also been simulated in DGVMs and other climate-driven terrestrial models (Randerson et al., 1996;Seidl et al., 2008;Wolf et al., 2008;Albani et al., 2010;Schäfer et al., 2010;Chen et al., 2015) less comprehensive than DVLSMs. However, most previous studies lacked realism by representing insect damage as end-of-year instantaneous events (instead of simulating their unfolding over many weeks during the growing season) and/or by imposing the assumed consequences of insect activity (e.g., reduced total canopy conductance) rather than letting 5 the model estimate these changes as a function of the new, insect-modified state of the vegetation. Moreover, many previous studies considered a single insect species, limiting their potential for global-scale studies, and failed to provide sufficient detail on the simulation of insect damage to efficiently guide modellers wanting to add insect disturbances to other DVLSMs. 10 Here, we present the "Marauding Insect Module" (MIM) we developed to simulate insect activity in DVLSMs and address the shortcomings identified above. MIM simulates insect activity with an approximated intra-annual schedule, prescribes only the plant damage caused directly by insects, and contains templates to allow for the inclusion of different insect functional types (IFTs). The concept of IFTs allows simplification of the huge diver- 15 sity of insect species by grouping species that cause similar impacts (Cooke et al., 2007;Arneth and Niinemets, 2010), and has recently been applied under the name of "Pathogen and Insect Pathways" in a simple ecophysiological model (Dietze and Matthes, 2014). We then illustrate, using MIM coupled to an existing DVLSM, the effects of a simulated MPB outbreak on many variables related to vegetation dynamics and exchanges of carbon, en- 20 ergy, and water, over daily to centennial timescales, and compare the results obtained to previous studies.

Overview
MIM was developed to be embedded within a host DVLSM and simulate the effects of in- 25 sect activity on vegetation dynamics, and biogeochemical and biogeophysical exchanges.

4
The underlying philosophy of MIM is to prescribe only the direct damage to the vegetation caused by insect activity, letting the host DVLSM quantify the resulting consequences for the post-mortality competition among the different vegetation types and the exchanges of carbon, energy, water, and momentum, based on the new conditions in the grid cells affected. Prescribing insect activity is less sophisticated than its prognostic simulation, but 5 nevertheless allows relevant questions to be addressed concerning the climatic and ecological impacts of insect-caused plant damage. We designed MIM so that it could be implemented in other DVLSMs in addition to the Integrated BIosphere Simulator (IBIS) we used in the current study. Furthermore, the structure of MIM is sufficiently flexible to allow for the representation of different insect species based on the templates we developed for three 10 IFTs.

Integrated BIosphere Simulator (IBIS)
We provide here only a short description of IBIS and refer readers to Foley et al. (1996) and Kucharik et al. (2000) for more details. IBIS represents two vegetation canopies (trees in the upper canopy, shrubs and grasses in the lower canopy), multiple soil layers (six in this 15 study, down to a depth of 4 m), and three snow layers when needed; both canopies intercept water and snow. Exchanges of radiation (shortwave and longwave), latent and sensible heat fluxes, and momentum between the atmosphere and the surface depend upon the state of each canopy. Water exchanges with the atmosphere consist of evaporation from intercepted water and the soil surface (including snow), as well as plant transpiration that is calculated 20 consistently with photosynthesis and removes moisture from each soil layer according to an exponential root profile. Fluxes of heat and moisture between soil layers, with drainage at the bottom, are influenced by the soil texture class, which is provided as input data. A time step of 60 min is sufficient to update all fluxes and state variables in offline (i.e., not coupled to a climate model) simulations. 25 IBIS represents vegetation diversity through a limited set of plant functional types (PFTs) characterized by different climatic constraints and physiological parameters. Photosynthesis and autotrophic respiration are computed on the same time step as land surface physics 5 (i.e., 60 min in this study) as a function of incoming radiation, CO 2 and O 2 concentration, temperature, and soil moisture stress. Changes in vegetation structure, including the proportions of competing PFTs, are determined at the end of each year, except for the leaf phenology of deciduous PFTs that is updated daily. Competition among PFTs accounts for the two-strata structure of vegetation (i.e., trees capture light first, but grasses have pref-5 erential access to water as they have a higher proportion of their roots in the upper soil layers) and is based on the annual carbon balance of each PFT. Litterfall is estimated annually based on PFT-specific parameters for each biomass pool and partitioned into daily transfers to the soil over the following year. Carbon decomposition and transfers among the different soil pools, which are influenced by microbial biomass and soil temperature and 10 moisture, are computed daily.
IBIS is arguably the first DVLSM to have been fully coupled to an Atmospheric General Circulation Model (Foley et al., 1998). Previous studies have shown that IBIS results compare reasonably well with observations, both over large regions (Foley et al., 1996;Kucharik et al., 2000;Lenters et al., 2000) and for field sites around the world (Delire and Foley, 1999;15 El Maayar et al., 2001Kucharik et al., 2006). Model intercomparisons also demonstrated that IBIS results were similar to other DGVMs, except for a stronger CO 2 fertilization with version 2 of the model (Cramer et al., 2001;McGuire et al., 2001;Friedlingstein et al., 2006).
We downloaded source code for version 2.6b4 of IBIS from the Center for Sus-20 tainability and the Global Environment (SAGE) website (http://nelson.wisc.edu/sage/ data-and-models/lba/ibis.php) with the required input data for climate (modified from the Climate Research Unit dataset CRU CL version 1.0 (New et al., 1999) by SAGE researchers for compatibility with IBIS) and for soil texture (based on an International Geosphere-Biosphere Program (IGBP) global dataset). The climate input data consist of different vari- 25 ables related to temperature, humidity (including precipitation and cloud cover), and wind speed. These climate data, which were provided for each month at a spatial resolution of 0.5 • , are temporally downscaled by a random weather generator built into IBIS to simulate daily and hourly variability (see Kucharik et al., 2000, for more details). We modified the 6 IBIS code before performing the illustrative simulations in Canadian forests (see Sect. 3) as follows: 1. We replaced the IGBP global soil dataset with survey data from the Soil Landscapes of Canada, versions 2.1 and 2.2, provided by the Canadian Soil Information System (http://sis.agr.gc.ca/cansis/nsdb/slc/index.html). 5 2. We modified the soil spin-up procedure due to the long time needed to reach equilibrium in Canada, the new procedure now taking 400 years (see Appendix A).
3. We improved the leaf-to-canopy scaling procedure for photosynthesis and transpiration, by: (1) replacing a mathematical simplification with the exact expression; and (2) adjusting the code that was used for the scaling integral (see Appendix B). Although 10 the current study used a constant CO 2 concentration, it is worth noting that these changes reduced the strength of CO 2 fertilization in IBIS.
4. We slightly increased (from 2.5 to 2.7 years) the mean carbon residence time for the needle pool of the boreal needleleaf evergreen PFT, which resulted in a better spatial distribution of the PFTs that exist in Canada, as well as a better succession dynamics 15 among these PFTs when starting a simulation from bare ground.

5.
We fixed an error in the random weather generator code that had previously prevented consecutive wet days from ever occurring.
6. We modified various elements related to energy exchanges: (1) we updated the nearinfrared optical properties of the lower-canopy leaves, based on values from version 20 4.0 of CLM (Oleson et al., 2010); (2) based on empirical data (Wang and Zeng, 2008), we constrained the variation of snow albedo as a function of solar zenith angle; and (3) we decreased the visible and near-infrared snow albedo parameters (see Appendix C). Following these changes, IBIS results for land surface albedo over Canada better matched MODIS-based values, both with (Barlage et al., 2005) and without 25 (MOD43B3-derived Filled Land Surface Albedo Product) snow cover. 7 7. We added a subroutine to confirm that the full annual carbon cycle, including the effect of MIM, balanced to a numerical precision of at least 1 × 10 −5 kg C m −2 .

Marauding Insect Module (MIM)
MIM aims to represent the effect of insect activity, from both outbreaking and nonoutbreaking insect species, on the coexistence of different PFTs and the land-atmosphere 5 exchanges of carbon, energy, water, and momentum. MIM does not currently simulate insect population dynamics, hence user-prescribed damage levels on defoliation and mortality (both in %) are required each year for each grid cell. It is the user's responsibility to ensure that prescribed damage levels over multiple years or grid cells are appropriate and that, for defoliators, prescribed vegetation defoliation and mortality are consistent with each 10 other (e.g., a single 5 % defoliation event is very unlikely to result in 80 % mortality). For host DVLSMs that, like IBIS, do not represent many individuals of the same PFT, a 5 % defoliation event translates into 100 % of the trees losing 5 % of their leaf area; in other DVLSMs, this same 5 % defoliation event could be assigned differently, for example by removing 100 % of the leaf area from 5 % of the trees. For each year and grid cell, MIM then 15 implements all the required changes in vegetation characteristics. The only input data for MIM are the prescribed levels of annual insect-caused defoliation and mortality, and the only state variables of the host DVLSM directly modified by MIM are the biomass and litter pools (to conserve carbon, new variables tracking insects respiration and biomass must also be added to the host DVLSM; see below). In fact, for DVLSMs that, unlike IBIS, simulate PFT 20 mortality explicitly (e.g., as a function of carbohydrates reserves), MIM would not need input data on prescribed mortality in the case of defoliators. We designed MIM to operate with a daily time step to simulate the intra-annual unfolding of insect activity and the resulting impacts, without the undue complications that would have stemmed from a sub-daily time step. Nevertheless, MIM could be adjusted to work under a shorter or longer time step. 25 MIM can currently simulate the activity from three IFTs parameterized to represent major outbreaking insect species in forests of North America: 8 -IFT #1: based on the forest tent caterpillar (Malacosoma disstria Hübner) and the gypsy moth, can defoliate (daily damage) and kill (year-end damage) broadleaf deciduous (BD) trees.
The user can prescribe damage from different IFTs to occur concurrently within the same grid cell, but for simplicity a given PFT cannot currently be targeted by more than one IFT. 10 For each IFT, the daily damage (defoliation for IFTs #1 and #2, mortality for IFT #3) unfolds by the same amount each day over the pre-defined duration of insect activity, thereby reaching the user-prescribed value at the end of the annual period of insect activity. The daily damage level (damage, in %) for a specific day d is thus given by: 15 where input user is the user-prescribed damage level for the year (in %), duration IFT is the duration of insect activity during the year (in days), and start IFT is the specific day of the year when insect activity starts.
Since MIM does not model insect population dynamics, we used fixed parameters for the values of start IFT and duration IFT (see Table 1 for values and corresponding literature 20 sources), except for start IFT of IFT #1: in this case, the activity begins on the same day as the IBIS-simulated beginning of leaf onset for the target tree, in accordance with the degree of synchrony between these two events for broadleaf defoliators (Dukes et al., 2009;Foster et al., 2013). (The duration of leaf onset simulated by IBIS is much shorter than duration IFT for IFT #1, so there is no risk that this defoliator of deciduous trees will consume leaves 25 9 faster than their simulated onset.) In reality, the start and duration of annual insect activity depend upon the phenological development of insects, for example the ending of the annual dormancy period for diapausing insects. Similarly, the linear unfolding of insect activity (i.e., equal day-to-day damage over the entire duration; see Eq. 1) is a simplification that could be refined in future implementations of MIM; yet, it provides a reasonable approximation of 5 the intra-annual progression of damage caused by the IFTs considered (Régnière and You, 1991;Cook et al., 2008;Hubbard et al., 2013). For example, although the individual feeding rate for the fifth and sixth larval instars of the eastern spruce budworm is much higher than for younger instars, the decreasing population density throughout summer leads to an approximately linear progression of total defoliation (Régnière and You, 1991). 10 Each day, the relevant biomass pools (leaves for IFT #1, needles for IFT #2, and all biomass pools for IFT #3) are decreased as a function of damage(d). More precisely, in IBIS-MIM damage(d) is multiplied by the "equilibrium values" (without insect damage) of the relevant biomass pools, and the result is then subtracted from the current value (on day d) of the relevant biomass pools. This approach was implemented here, because IBIS com- 15 putes these "equilibrium values" at the end of the previous year, when updating vegetation structure and proportions of competing PFTs; in other DVLSMs, however, this specific element of MIM's implementation might need to be adjusted. Besides daily defoliation, IFTs #1 and #2 can kill trees (also according to user-prescribed damage levels); when this happens, mortality of the PFT targeted by IFT #1 or #2 occurs as a one-time event at the end of the 20 year. We explain below how MIM deals with trees killed during a given year, either through daily (IFT #3) or sudden (IFTs #1 and #2) simulated mortality. Note that PFTs entirely defoliated by IFTs #1 or #2 behave exactly as dead trees if no reflush is allowed (see below), even if these killed PFTs are not labelled as "dead" before the end of the year.
The carbon contained in leaves or needles consumed by IFTs #1 and #2 based on 25 damage(d) needs to be accounted for to obey the conservation laws that form the basis of DVLSMs. Consequently, MIM divides all the defoliated carbon among three pathways: respired (i.e., immediately returned to the atmosphere as CO 2 ), excreted as frass that is then treated as leaf/needle litterfall by IBIS, or stored in IFT biomass (see Table 1). This last variable will be very relevant if MIM is eventually expanded to simulate insect population dynamics; currently, the biomass of defoliator IFTs is simply exported out of the simulation domain at the end of each year, and IBIS net ecosystem carbon balance accounts for this export, as well as IFT respiration. At present, MIM does not quantify the stem carbon consumed by IFT #3 and the resulting IFT biomass; given the difference between the total 5 biomass of bark beetles and the biomass of the trees they killed, this should have very small impacts on the simulated carbon fluxes. Many tree species can produce a second flush of foliage after an early-season defoliation event (Jones et al., 2004;Schäfer et al., 2010). We therefore allowed for the possibility of reflush in MIM, as this phenomenon can substantially influence simulated land-atmosphere 10 exchanges and vegetation competition. The amount of reflush (in %) occurring during day d is given by: where total reflush is the total amount of leaf reflush (in % of the total leaf biomass lost to defoliation earlier in the same year), duration reflush is the duration of the reflush (in days), 15 and start reflush is the specific day of the year when reflush starts (see Table 1; please note that reflush starts after defoliation is completed). Each day, the leaf biomass pool of the defoliated PFT is then increased based on the value of reflush(d) and the total amount of defoliation before the reflush. Although duration reflush is currently determined by phenology algorithms from IBIS, approaches based on empirical data could be implemented instead. 20 The value of total reflush for the year can be set to zero to prevent unrealistic reflush when the defoliation level is very low, or when trees have already been weakened by previous defoliation events, or if mortality is also prescribed for the same year. When mortality is prescribed, MIM also needs to account for the carbon remaining in IFTkilled trees, both for mortality simulated as a sudden event at the end of the year (IFTs #1 25 and #2) and for daily mortality (IFT #3). We therefore added a new feature to IBIS, whereby a PFT killed by an IFT instantaneously becomes a dead standing tree (DST) conserving 11 the same carbon pools. DSTs interact with energy, water, and momentum exchanges in the same way as live PFTs (e.g., interception of precipitation and absorption of radiation), but do not transpire or contribute to canopy photosynthesis. The simplest approach to simulate the fate of DSTs would have been to transfer all their carbon to IBIS litter pools at the end of the year when mortality happens. However, this would cause unrealistically large and sudden 5 changes in litterfall and canopy structure, because insect-killed trees initially remain standing and fall gradually on the forest floor. Consequently, the carbon contained in DST pools is progressively transferred to the appropriate litter pools based on a prescribed schedule. MIM currently offers five possible schedules corresponding to the snagfall dynamics of different tree species: 10 -Case #1: BD tree PFT killed by IFT #1, fate of DST based on trembling aspen (Populus tremuloides Michx.) in eastern Canada (Angers et al., 2010).
-Case #2: BD tree PFT killed by IFT #1, fate of DST based on trembling aspen in western Canada (Hogg and Michaelian, 2015). The transfer of carbon from DST pools (i.e., fine roots, leaves/needles, and stems, the latter including coarse roots and branches) towards IBIS litter pools starts after a delay period and then occurs at a specific rate; note that the delay refers to the time of attack, even for Case #5. Table 2 gives the value of these parameters for the five cases currently a one-time event, at the end of the year of mortality (note that IBIS partitions all annual DST transfers into daily amounts over the following year). For deciduous PFTs (i.e., Cases #1 and #2), the transfer of DST leaves also occurs as a one-time event. On the other hand, the DST needles are transferred to litter pools over many years for evergreen PFTs (i.e., Cases #3-5). Finally, the DST stems are also transferred to litter pools over many years, 5 usually starting after a 5-year delay period (see Fig. 1). As with the IFT-related parameters, all these aspects of DST dynamics can easily be modified as a function of new data or to accommodate other tree species. Moreover, for large-scale studies, the IFT-and DSTrelated parameters could vary spatially to reflect within-species variation, instead of having uniform values as we have used here (e.g., needlefall for Case #5 could occur over more 10 than three years).

Simulation design
To illustrate the performance of IBIS-MIM, we conducted six simulations using the MPBinspired IFT (i.e., IFT #3 from Table 1) with DST dynamics based mostly on MPB-killed 15 lodgepole pine (i.e., Case #5 from Table 2). We performed an outbreak simulation and a control simulation in each of three different locations in British Columbia, Canada, henceforth designated as the northern, central, and southern grid cells (Table 3). These three locations, which we used as proxies to assess the influence of climate on the main outcomes, have suffered substantial MPB-caused mortality since 2000 (Walton, 2013). The 20 mean annual temperature was almost equal in the northern and central grid cells, but summer was warmer and winter was colder in the former; the southern grid cell was warmer throughout the year. Annual precipitation was very similar in the three grid cells, but summer rainfall was substantially lower in the southern grid cell, leading to lower soil water content during the growing season. 25 13 All simulations started with the new 400-year spin-up procedure and were performed under a constant climate. In each grid cell, we prescribed a single 100 % mortality event happening in year 401 (i.e., in year 1 following the spin-up period) and continued the simulation up to year 1000. This single, 100 % mortality event does not aim to represent actual MPB outbreaks, but was implemented for the sake of simplicity, to increase the signal-to-5 noise ratio of the results, and to test the theoretical upper limit of impacts. We used the same climate data and weather generation for the outbreak simulation and the no-mortality control simulation performed in a given grid cell. In addition to yearly results throughout the entire simulation, we saved daily (monthly) results during 10 (200) years after the mortality event. We excluded the boreal BD tree PFT from simulations due to the generally low den-10 sity of such trees within MPB-attacked stands in British Columbia (Hawkins et al., 2012). Consequently, competition took place among four different IBIS PFTs: boreal NE trees (i.e., the target PFT), evergreen shrubs, cold-deciduous shrubs, and C 3 grasses. Figure 2 shows the effect of the single MPB outbreak on net primary productivity (NPP) in 15 the three grid cells. In all cases, simulated NPP of the target NE trees decreased to zero while NPP of the lower canopy substantially increased following the 100 % mortality event;

Responses over different timescales
the productivity of the different PFTs then gradually resumed towards the pre-outbreak levels (Fig. 2a, c, and e). However, the growth release of the lower canopy was much stronger in the northern and central grid cells than in the southern grid cell, where conditions were 20 drier during the growing season. Such positive impacts on lower canopy have often been reported following outbreaks from MPB and other bark beetles (Stone and Wolfe, 1996;Griffin et al., 2011;Simard et al., 2011;Bowler et al., 2012;Brown et al., 2012;Vanderhoof et al., 2014).
In the northern and central grid cells, the lower canopy growth release exceeded the pro- 25 ductivity losses coming from the death of NE trees, so that total post-outbreak NPP soon exceeded NPP in the control runs ( Fig. 2b and d). The increase in ∆NPP in the central grid cell from year ∼ 600 onwards came from the impact of the outbreak on the competition balance among PFTs: although NPP seemed relatively stable at the end of the spin up (i.e., years 300-400 in Fig. 2c), lower canopy NPP decreased markedly between years 600 and 750 in the control simulation, whereas the MPB outbreak released the lower canopy and postponed this decline, which seems 'bound to happen' in the long term. Empirical Belovsky and Slade, 2000) and modelling (Mattson and Addy, 1975;5 Seidl et al., 2008;Albani et al., 2010;Pfeifer et al., 2011;Hansen, 2014) studies of insect damage have previously shown that total productivity, biomass, or carbon storage can be higher in disturbed than in undisturbed ecosystems. As was the case in IBIS-MIM, the mechanisms identified in these previous studies involved responses from non-target vegetation, i.e., other species or non-attacked age classes of the target species. In the southern 10 grid cell, on the other hand, total NPP was reduced for about 75 years and then increased marginally for a few decades, before returning to the control level (Fig. 2f).
The previous results also exhibited an interesting feature: in all grid cells, the recovery of the NE trees was initially very rapid, but was then reversed after ∼ 20-25 years before resuming again (Fig. 2a, c, and e). Although additional simulations would be required to con-15 firm our hypothesis, we believe that this "dip" came from indirect biogeophysical interactions between recovering NE trees and decaying DSTs in the relatively cold climate considered here. After MPB mortality, the interception of radiation (shortwave and longwave) by DSTs warmed the surrounding air, allowing photosynthesis in the recovering NE trees to occur faster at a higher temperature than if DSTs had been absent. As DSTs gradually fell, NE 20 trees captured more light but had a lower needle temperature, resulting in lower NPP. Such strong photosynthesis-temperature responses have been found to play a major role when simulating future vegetation dynamics (Sitch et al., 2008;Medvigy et al., 2010) and carbon cycle-climate feedbacks (Matthews et al., 2005). Figure 3 shows the impact of the outbreak on four variables (two related to carbon cy- 25 cling, one to energy exchanges, and one to water cycling) over different timescales (yearly, monthly, and daily). The changes in net ecosystem productivity (NEP; Fig. 3a) were driven mostly by NPP, including the increases in total NPP ∼ 5 years post-mortality in the northern and central grid cells. Changes in heterotrophic respiration (R h ) were generally smaller, but contributed to the NEP local minimum around year 25 (particularly visible in the central and southern grid cells) and progressively offset the NPP increase in the northern and central grid cells, so that ∆NEP became negligible after roughly a century. The total amount of aboveground litter (Fig. 3b) slightly decreased for a few years after the mortality event, because the total litterfall from DSTs in the outbreak simulations was initially lower than from 5 live trees in the control simulations. After a few years, however, the situation was reversed and the increase in aboveground litter was > 1.5 kg C m −2 in all grid cells ∼ 25 years after the mortality event, gradually decreasing afterwards. After about 75 years, the aboveground litter was lower in the outbreak simulations due to the reduced litterfall from the still recovering vegetation. The pre-outbreak equilibrium was reached about three centuries after the 10 mortality event. The monthly albedo (Fig. 3c) increased during the initial years as the needles fell from DSTs. The impact of snow cover was clearly apparent in the yearly cycle of albedo changes, with much higher albedo increases during winter months. The few points showing a decrease in albedo resulted from the earlier snowmelt in the outbreak simulations, a response that is illustrated for the central grid cell (Fig. 3d). While the snow amount 15 was slightly higher following the first snowfall events (barely visible in Fig. 3d), in the middle of winter the control grid cells generally had more snow. But above all the snowmelt started and finished much earlier in the outbreak simulations, by about three weeks in the case illustrated. 20 25 other disturbances (girdling or clearcutting) for the effect. We note, however, that the identification of appropriate control stands for field and satellite studies is not a straightforward task, which may partly explain why the qualitative impact (increase, no change, or decrease) of MPB mortality varied across studies for some variables. Furthermore, the level of stand mortality differed among studies or was not quantified and, except for a few modelling studies, was less than the 100 % mortality simulated in IBIS-MIM. These limitations prevented us from performing more quantitative comparisons. The comparisons covered 28 different variables related to carbon cycling and veg- 5 etation dynamics, energy exchanges, and the water cycle. These comparisons further spanned various timescales: annual (all variables related to carbon cycling and vegetation dynamics, albedo, evapotranspiration, runoff, and soil moisture), seasonal/monthly (all variables related to energy exchanges, evapotranspiration, transpiration, soil moisture, snow depth/amount, and snowmelt onset), and daily (peak flow, snow depth/amount, and 10 snowmelt onset). Among the 28 variables, IBIS-MIM prescribed only the snagfall dynamics of DSTs. IBIS-MIM results agreed with previous studies for most variables, thereby illustrating that the model constitutes an appropriate tool for studying the impacts of insect-induced plant damage on many inter-dependent variables spanning a large range of timescales. For most variables related to carbon cycling and vegetation dynamics, the qualitative 15 responses of IBIS-MIM changed over time for two reasons. First, as seen in Sect. 3.2, lower canopy biomass substantially increased following the canopy opening in the northern and central grid cells (but much less in the southern grid cell), eventually reversing the initial response for GPP, NPP, R a , R s , NEP, and total LAI (abbreviations are defined in Table 4) in these two locations ∼ 5-15 years after the MPB outbreak. Note that the simulated increase 20 in shrub biomass was marginal in the three grid cells, but that grass biomass increased substantially. Lower canopy fractional cover increased in the northern grid cell only, because this variable was already at its maximum value before the mortality event in the other two grid cells. Second, the prescribed snagfall dynamics of DSTs led to a carbon response over multiple timescales  that affected R h , NEP, and aboveground litter (see 25 also Fig. 3a and b). Among the variables related to energy exchanges, IBIS-MIM responses for temperature and albedo systematically agreed with previous studies. (The air temperature in field studies was measured close to breast height, a level at which IBIS-MIM does not estimate tempera-ture. As a proxy, we used the mean of the simulated temperature responses in the middle of the upper and lower canopies.) These responses became particularly strong and sustained after the complete fall of needles from DSTs. We note that the impacts on temperature variables in IBIS-MIM were generally opposite between winter and summer; unfortunately, none of the previous studies reported wintertime temperature changes. For latent and sensible 5 heat fluxes, however, IBIS-MIM differed noticeably from previous studies: after the year of mortality, summertime latent heat flux actually increased for three years in the southern grid cell and for much longer in the other grid cells. The pattern was the opposite for summer sensible heat, except in the southern grid cell where this variable did not show a systematic behaviour. We think that these responses for summer turbulent heat fluxes had two differ-10 ent causes. For 1-4 years after the mortality event, the higher summer latent heat flux in all grid cells came from a major increase in evaporation which, in turn, probably resulted from the combination of two pre-existing biases in the land surface module (LSX) that computes the exchanges of energy, water, and momentum within IBIS-MIM: (1) the overestimation of upper soil temperature in summer (El Maayar et al., 2001), which likely increased follow- 15 ing the mortality event; and (2) the overestimation of heat storage within stems (including DSTs in our simulations), leading to an overestimated nighttime evaporation flux when the heat is released (Pollard and Thompson, 1995;El Maayar et al., 2001). For ≥ 5 years after the mortality event, the increase in summer latent heat flux in the northern and central grid cells rather resulted from the strong growth of grasses mentioned previously. Indeed, NPP 20 and LAI of grasses in these grid cells were then large enough to overcompensate for the decreases due to tree mortality, resulting in higher total transpiration, evapotranspiration, and latent heat flux. In the southern grid cell, where the response from grasses was much smaller, summer latent heat flux decreased ≥ 5 years post-mortality.

Evaluation of performance
While acknowledging possible issues with these IBIS-MIM results, particularly 1-4 years 25 post-mortality, we want to underline the limitations from previous studies on turbulent heat fluxes and the closely related evapotranspiration. Three of the four modelling studies (Wiedinmyer et al., 2012;Mikkelson et al., 2013b;Chen et al., 2015) indirectly "forced" the responses they obtained by directly reducing the total canopy conductance without account-ing for the possible growth release of the surviving vegetation, while the fourth modelling study (included in the Mikkelson et al., 2013a, review) only computed the change in runoff and then assumed no change in soil moisture to estimate the change in evapotranspiration. The two satellite-based studies rest upon the highly-parameterized MODIS evapotranspiration dataset (Mu et al., 2011), which has not been developed and tested in the context of 5 MPB-killed forests. The only field-based study on evapotranspiration (included in the Mikkelson et al., 2013a, review) also neglected possible changes in soil moisture. Furthermore, other field-based studies -not included in our comparison due to their lack of control data -found very little change in evapotranspiration over various years following MPB mortality for sites located close to the northern grid cell Brown et al., 2014), 10 or found that evapotranspiration increased over three years despite an ongoing increase in MPB mortality at a site in Wyoming, US (Reed et al., 2014). For water cycle variables besides evapotranspiration, the agreement with previous studies was also fairly good. The soil water budget in IBIS-MIM is very sensitive to the distribution of precipitation events during each month, so the responses were highly variable 15 for runoff, peak flow, and soil moisture, particularly in the southern grid cell. Nevertheless, the responses provided in Table 4 were observed over the first ∼ 5 years. Afterwards, runoff remained higher in the outbreak simulation for the southern grid cell (resulting in part from the faster snowmelt), but became smaller in the other grid cells due to the increase in evapotranspiration as leaf area expanded. A field study on drought-induced tree mortality also 20 linked an unexpected decrease in annual runoff to a growth release of the lower canopy (Guardiola-Claramonte et al., 2011). Peak flow, on the other hand, remained overall higher in all grid cells for at least a decade. After an initial increase lasting ∼ 5 years, soil moisture showed a sustained decrease, likely caused by the snowmelt-related higher runoff in the southern grid cell and by the higher evapotranspiration in the other grid cells. Although of- 25 ten slightly higher at the beginning of the season, snow depth/amount overall decreased in IBIS-MIM (see Fig. 3d), contrary to most previous studies. This outcome likely resulted from the overestimated heat storage in DSTs and could lead to the simulated snow cover season ending too early. Yet areal snow coverage, which matters most for albedo, was equal for the control and outbreak simulations during most of the snow cover season and, most importantly, the earlier onset of snowmelt agreed with the majority of previous studies and was of reasonable magnitude.
We checked whether the outcomes presented in Table 4 were sensitive or not to the specific weather simulated by performing two additional replicates for each grid cell. We 5 found that the qualitative outcomes were the same for all variables, except for one minor difference: for one of the two replicates in the central grid cell, the post-outbreak fractional cover of the lower canopy increased slightly because it was not already at its maximum value, contrary to the case reported in Table 4. The quantitative results were also very similar across replicates, except for some water-related variables that are very sensitive to 10 the exact timing of precipitation events.
Finally, although assessing IBIS was not the point of this study and has already been done elsewhere (Foley et al., 1996;Delire and Foley, 1999;Kucharik et al., 2000Kucharik et al., , 2006Lenters et al., 2000;El Maayar et al., 2001, the results obtained for the three grid cells compared favourably to recent studies, with a small underestimation of biomass (Beau- 15 doin et al., 2014) and NPP (Gonsamo et al., 2013). Obtaining reliable data on soil carbon is notoriously difficult; when compared to the Harmonized World Soil Database (down to a depth of 1 m) as provided by Exbrayat et al. (2014), IBIS apparently overestimated soil carbon (down to a depth of 4 m), at least in the southern grid cell, even when accounting for the fact that a substantial fraction of soil carbon is found at a depth greater than 1 m 20 (Jobbágy and Jackson, 2000).

Discussion
Many previous studies have represented insect damage in DVLSMs or less comprehensive climate-driven terrestrial models (Randerson et al., 1996;Krinner et al., 2005;Seidl et al., 2008;Wolf et al., 2008;Albani et al., 2010;Schäfer et al., 2010;Edburg et al., 2011;Medvigy 25 et al., 2012;Mikkelson et al., 2013b;Chen et al., 2015). To our knowledge, however, our study is the first to assess, over daily to centennial timescales, the impacts from insect 20 damage on vegetation dynamics and the carbon, energy, and water cycles in an integrated way (see Sect. 3). We compared the qualitative impacts of a simulated MPB outbreak on 28 IBIS-MIM variables with many field-, satellite-, and modelling-based studies (see Table 4), finding an overall good level of agreement. Our results further suggest that the physical presence of DSTs can benefit vegetation regrowth due to their interactions with radiation. 5 A previous study also showed that falling DSTs can impact tree recovery through altered soil nitrogen dynamics . Since DSTs contribute substantially to the biogeophysical and biogeochemical legacies of insect outbreaks, they should be explicitly modelled when feasible.
We developed MIM to account for the major processes related to insect activity (Table 1), 10 including the dynamics of DSTs (Table 2) when applicable. The generic design of the module could serve as a template to represent other IFTs and/or DSTs, and should facilitate future developments such as replacing the prescribed intra-annual unfolding of insect activity with algorithms based on simulated insect phenology. Moreover, MIM could be modified by simulating the fall of DSTs probabilistically (e.g., as in FireBGCv2; Keane et al., 2011) or 15 enhanced by: simulating the fall of DSTs as a function of environmental conditions (Lewis and Hartley, 2005) or the size of DSTs (e.g., as in FVS; Rebain et al., 2010); reducing snow albedo when needles fall from DSTs (Pugh and Small, 2012); and accounting for changes in needle optical properties as they turn from green to red (Wulder et al., 2006). The simple structure of MIM should also facilitate the adaptation of the module to other 20 DVLSMs. Of course, MIM will then reflect many of the strengths and weaknesses of its host model. For example, the parameters of the boreal NE PFT in IBIS 2.6b4 were not based on lodgepole pine specifically. Furthermore, IBIS simulates a single boreal NE PFT, whereas different NE tree species can coexist in MPB-attacked stands (Hawkins et al., 2012). Since IBIS does not represent different age cohorts within the same PFT, the model cannot ac- 25 count for the fact that MPB generally targets the larger trees (Axelson et al., 2009;Hawkins et al., 2012;Hansen, 2014). For < 100 % mortality, the responses of surviving younger trees would likely differ from those of surviving mature trees and could enhance the recovery of the target species. Impacts on tree demographics might also lead to complex stand-level responses, for example increasing total biomass despite reduced productivity because of a strong decrease in competition mortality . Other shortcomings of IBIS that affected IBIS-MIM results came from the apparent overestimation of stem heat storage (Pollard and Thompson, 1995;El Maayar et al., 2001) and the absence of carbon-nutrient interactions (Edburg et al., , 2012Mikkelson et al., 2013a). On the other hand, IBIS 5 two-strata vertical vegetation structure and detailed biophysics computations, both inherited directly from the LSX land surface module (Pollard and Thompson, 1995), allowed the lower canopy growth release and the biogeophysical impacts of DSTs presence to be simulated more realistically than with many other DVLSMs. Finally, the strong link between climate and insect life cycles (Dukes et al., 2009;Bentz 10 et al., 2010) provides incentive for eventually enhancing MIM by including process-based representations of insect population dynamics in DVLSMs (Fosberg et al., 1999;Arneth and Niinemets, 2010), rather than prescribing insect damage through input data.

Conclusions
Insect damage to vegetation triggers major interacting effects on the cycles of carbon, 15 nutrients, energy, and water, and also affects atmospheric chemistry. Given that Dynamic Vegetation-Land Surface Models (DVLSMs) were designed to simulate coupled biogeophysical and biogeochemical fluxes within a consistent framework that accounts for changes in vegetation state, these models appear as good candidates to assess many of the consequences from insect-induced vegetation damage over a wide range of timescales. 20 Here, we presented version 1.0 of the Marauding Insect Module (MIM) developed to simulate, within the Integrated BIosphere Simulator (IBIS) DVLSM, the impacts of prescribed levels of annual insect damage. MIM currently includes three insect functional types (IFTs) broadly representing defoliators of broadleaf trees, defoliators of needleleaf trees, and bark beetles. The parameterization of IFTs was based on key outbreaking insects affecting North 25 American forests, but could be modified to represent other insect species, effects on other vegetation types (e.g., agricultural fields), and, with further adjustments, effects of some vegetation pathogens (e.g., Dietze and Matthes, 2014). Similarly, the fate of the insect-killed dead standing trees (DSTs) can easily be adjusted to go beyond the five cases currently implemented. Finally, MIM itself was designed in such a way that it should be transferable to other DVLSMs with limited adjustments.
We also illustrated the realism and usefulness of IBIS-MIM by simulating a 100 % mor-5 tality event caused by the mountain pine beetle at three locations within British Columbia, Canada. First, we looked at the impacts of the outbreak on a variety of processes spanning daily to centennial timescales. One interesting outcome from this assessment is that DSTs intercept radiation and therefore warm the surrounding air, which in a cold climate could be beneficial for tree recovery. Second, we found that IBIS-MIM agreed qualitatively 10 with the results from 38 field-, satellite-, and model-based studies for 28 different variables related to vegetation dynamics, and exchanges of carbon, energy, and water. These outcomes supported the idea that DVLSMs are valuable tools to study the consequences from insect-induced plant damage.
Insect outbreaks, but also less spectacular background-level vegetation damage caused 15 by insects, are part of the natural dynamics of terrestrial ecosystems worldwide. The use of IBIS-MIM and other similar process-based modelling tools suitable for climate-related studies should therefore help us better understand the wide range of possible impacts of insects on several processes in the Earth system, for past, current, and future conditions.

Appendix A: Soil spin-up procedure
The previous soil spin-up procedure lasted 150 years and was performed as follows: 40 iterations of the soil module were repeated each year during the first 75 years; then, during the following 25 years, the number of iterations per year decreased linearly from 40 to 1; and finally, during the last 50 years, soil carbon pools were brought to equilibrium under 5 a single iteration per year. The total number of soil module iterations under this procedure was around 3500.
The new soil-spin up procedure lasts 400 years and is performed as follows: 80 iterations of the soil module are repeated each year during the first 350 years; then, during the following 40 years, the number of iterations per year decreases linearly from 80 to 1; and finally, 10 during the last 10 years, soil carbon pools are brought to equilibrium under a single iteration per year. The total number of soil module iterations under this procedure is around 29 600.
Appendix B: Leaf-to-canopy scaling

B1 The "extpar" simplification
The net photosynthesis (A n (X), in mol CO 2 s −1 m −2 of leaf) for a leaf that is X units into the 15 upper or lower canopy (where X is the cumulative vegetation (leaf plus stem) area index, in m 2 of vegetation m −2 of ground, with X = 0 at the top of the canopy) is computed as: where A n (0) is the photosynthesis for a leaf at the top of the canopy and A, B, C, k, and h are coefficients computed in IBIS. Previously, this expression was simplified to: 20 A n (X) = A n (0) exp(−extparX) with: We therefore worked directly with Eq. (B1) and removed the "extpar" simplification from the code. Note that this simplification might have been required in version 1 of IBIS, which used a different leaf-to-canopy scaling approach than version 2 (Foley et al., 1996;Kucharik et al., 2000).

B2 Leaf-to-canopy scaling integral
The total canopy photosynthesis (A n,canopy , in mol CO 2 s −1 m −2 of ground) is given by the following scaling integral: where LAI is the total canopy leaf area index and XAI is the total canopy vegetation (leaf 10 plus stem) area index. Previously, the LAI/XAI factor was removed from the integral above and was included in the computation of the photosynthetically active radiation absorbed by leaves at the top of the canopy; the results for A n,canopy were then the same for light-limiting conditions, but not under Rubisco-limiting or CO 2 -limiting conditions. We therefore adjusted the code to work directly with Eq. (B4) under all conditions. Note that this adjustment and 15 the removal of the "extpar" simplification affected canopy transpiration, which is computed as a function of canopy photosynthesis.

Appendix C: Energy exchanges C1 Near-infrared optical properties of lower-canopy leaves
We modified the reflectance (unitless) from 0.60 to 0.40, and the transmittance (unitless) 20 from 0.25 to 0.30.

25
C2 Snow albedo vs. solar zenith angle IBIS increases snow albedo for solar zenith angles greater than 60 • , but these increases appeared too large for very high zenith angles. We therefore limited these increases to a maximum of 10 % above the value at 60 • for visible radiation and to a maximum of 15 % above the value at 60 • for near-infrared radiation.

C3 Visible and near-infrared snow albedo parameters
We decreased the following parameters related to snow albedo (unitless): low-temperature value in the visible (from 0.90 to 0.80), high-temperature value in the visible (from 0.70 to 0.60), low-temperature value in the near-infrared (from 0.60 to 0.50), and high-temperature value in the near-infrared (from 0.40 to 0.30).  n/a n/a total reflush 50 % of defoliation loss 15 n/a n/a duration reflush Typically ∼ 5 days a n/a n/a  a All transferred to litter on the year of mortality. b If some leaves/needles remain because mortality occurred with less than 100 % defoliation or reflush happened. c Rates are linear and start after the delay period; for stems, some cases have two consecutive linear periods showed on two lines: for each period, the duration (in years) and the total fraction transferred over the period (in %) are provided. 39 Table 4. Comparison of IBIS-MIM results for a simulated MPB outbreak with field-, satellite-, and model-based studies (increase: ↑; no change: -; decrease: ↓). Under the "Field", "Satellite", and "Model" columns, the numbers refer to the studies listed below. Under the "IBIS-MIM" column, the values in parentheses give the number of grid cells sharing the same qualitative results (only provided when the three grid cells differed).  Table 2). Mortality happened in year 0 and all values are for the end of the corresponding year.