More than the sum of its parts: how disturbance interactions shape forest dynamics under climate change

. Interactions among disturbances are seldom quanti ﬁ ed, and how they will be affected by climate change is even more uncertain. In this study, we sought to better understand how interactions among disturbances shift under climate change by applying a process-based landscape disturbance and succession model (LANDIS-II) to project disturbance regimes under climate change in north-central Minnesota, USA. Speci ﬁ cally, we (1) contrasted mortality rates and the extent of disturbance for four individual (single) disturbance regimes ( ﬁ re, insects, wind, or forest management) vs. all four disturbance regimes operating simultaneously (concurrent) under multiple climate change scenarios and (2) determined how climate change interacts with single and concurrent disturbance regimes to affect carbon stocks and forest composition. Under single disturbance regimes, we found that climate change ampli ﬁ es mortality, but did not substantially change the overall extent of disturbances. Tree mortality under the concurrent disturbance regime scenario was less than the sum of all single disturbance regimes, providing evidence of signi ﬁ cant negative feedbacks among disturbances, particularly under climate change. Finally, we found that climate change was the most critical driver of area harvested (via shifts in species composition), soil carbon, species composition, and diversity, while the disturbance regime (i.e., single or concurrent) had a larger in ﬂ uence on aboveground carbon and the relative dominance of conifers vs. hardwoods. In conclusion, our simulations suggest that disturbance interactions will be strongly mediated by climate change and will produce increasingly negative feedbacks, preventing worst-case disturbance outcomes. Our results underscore the importance of running simulations with multiple disturbances on the landscape concurrently rather than focusing on any one or two disturbances.


INTRODUCTION
Disturbances are an important driver of landscape change, and climate change is anticipated to cause more widespread and intense disturbances across a range of ecosystems (Westerling et al. 2011, Hicke et al. 2012. These changes in climate and disturbance regimes may cause widespread tree mortality that push ecosystems toward alternative stable states, potentially reducing forest health, diversity, or even forest land cover (Beisner et al. 2003).
Disturbances are often studied in isolation or as an interaction between two disturbances, focusing on either linked disturbances (i.e., a disturbance event that affects the likelihood or severity of subsequent events) or compound disturbances (i.e., overlapping disturbance events that affect system recovery rates; Buma 2015). The difficulty is that disturbance interactions can result in ecological surprises where system responses are unpredictable based on knowledge of each disturbance alone. For example, large bark beetle outbreaks in western North America prompted widespread concern that these infestations would lead to greater fire activity (Simard et al. 2011, Jolly et al. 2012, based on observed changes in fuel characteristics (Jolly et al. 2012, Klutsch et al. 2011, Simard et al. 2011). However, these concerns were later shown to be unsubstantiated. Douglas-fir, spruce, and western pine beetle outbreaks had no detectable effect on multiple measures of fire severity or extent in Wyoming, California, and Colorado (Kulakowski and Veblen 2007, Bond et al. 2009, Harvey et al. 2013, for exception, see Turner et al. 2000), with fire severity instead varying with weather, cover type, elevation, and/or slope position (Kulakowski andVeblen 2007, Harvey et al. 2013). In some cases, bark beetles may actually reduce overall wildfire potential if its host is also the primary source of ladder fuels and determinant of fire spread in the ecosystem (Scheller et al. 2017).
Nearly all forested ecosystems are affected by multiple disturbance types, but only a few studies have investigated three-way (or more) interactions among disturbance types (D'Amato et al. 2011, Bradford et al. 2012, and fewer still have examined how climate change will modify these interactions (Urban et al. 2000, Foster et al. 2002, Allen 2007, Gustafson et al. 2010. Climate change affects ecosystem dynamics at multiple levels, including vegetation growth, mortality, reproduction and competitive interactions, spatial patterns of resources such as water (Urban et al. 2000), and disturbance mechanisms such as ignition rates (Flannigan et al. 2009) and tree stress affecting insect damage (Raffa et al. 2008; Fig. 1). Disentangling the relative importance of interactions between disturbances is important, because both disturbance legacies and climate affect ecosystem recovery (e.g., recolonization, growth, succession).
In this study, we applied a process-based landscape disturbance and succession model (LANDIS-II) to investigate the independent and cumulative effects of disturbance on landscape dynamics across three climate change scenarios in north-central Minnesota, USA. Our objectives were twofold: (1) to contrast relative tree mortality rates and the extent of disturbance for four individual (single) disturbance regimes (fire, insects, wind, or forest management) vs. all four disturbance regimes operating concurrently (concurrent) under different climate change scenarios, and (2) to determine how climate change interacts with single and concurrent disturbance regimes to affect forest composition. We focused on a large landscape in north-central Minnesota (USA) because (1) all four disturbances are historically important, (2) climate change is expected to be substantial over the next century due to its relatively high latitude , and (3) there is a diverse suite of extant tree species that allows for emergent behaviors, as expressed through changes in succession via competition, facilitation, release, and compensatory growth (Xu et al. 2012).
Our simulation experiment evaluated three main predictions related to disturbance interactions. First, we hypothesized that climate change will amplify tree mortality, caused by increased insect mortality and larger fires (red arrows in Fig. 1). Second, tree mortality under the concurrent disturbance regime scenario will be less than the sum of all single disturbance regimes, because some mortality is compensatory, due to negative feedbacks (green arrows in Fig. 1). We expect compensatory effects to increase over time as disturbance likelihood generally increases with tree size/age for most of our disturbance types (i.e., harvesting, windthrow, and insects). For example, mortality caused by one disturbance type reduces the pool of survivors and changes forest composition, which can reduce the probability or severity of subsequent disturbances. Finally, at least over the next century, disturbances (green arrows) will modify forest composition more than climate change (red arrows); the relative importance of a disturbance type on forest composition will be in proportion to their extent and induced mortality.

Site description
Our study area consists of a 3.4 million-ha landscape in north-central Minnesota, which encompasses the Northern Minnesota Drift and Lake Plains Ecological Section (MDL or 212N), and the entirety of the Chippewa National Forest (CNF; Fig. 2). The climate is humid, continental, and cold temperate; mean temperatures are À15°C in January and 20°C in July; mean annual precipitation is 87 cm (PRISM Climate Group 2015). Soils are quite variable in texture, chemistry, stoniness, and drainage, resulting from heterogeneous glacial activity, bedrock, and climate across the landscape (Albert 1995). This variation in soils is reflected in the heterogeneous matrix of mesic forests, bogs, and hardwood forests that is also heavily influenced by management activities (e.g., red and jack pine (Pinus resinosa Ait., P. banksiana Lamb.) plantations). Mesic forests are widespread throughout the MDL with trembling and bigtooth aspen (Populus tremuloides and P. grandidentata Michx.), paper birch (Betula papyrifera Marshall), northern red oak (Quercus rubra L.), basswood (Tilia americana L.), and sugar maple (Acer saccharum L.). The eastern part of the MDL includes expansive bogs of black spruce (Picea mariana Mill.) and eastern larch (Larix laricina [Du Roi] K. Koch) and wetland forests of northern white cedar (Thuja occidentalis L.) and black ash (Fraxinus nigra Marshall). In the western MDL, there are Fig. 1. Conceptual framework for this study illustrating both direct (arrows) and indirect interactions among the forest structure and the four primary disturbance types under alternative climate change scenarios. Climate may affect forest composition and structure via its influence on the biophysical environment, affecting establishment and growth of tree species, and also effects on specific disturbances (i.e., fire and two-lined chestnut borer). These two disturbances types are fully dynamic, in that their regimes emerge entirely from underlying process. Fire disturbance is driven by the combination of climate and fuels (red + green arrows). The combination of forest tent caterpillar defoliation and drought (climate) triggers outbreaks of two-lined chestnut borer on its oak hosts (yellow + red + green arrows). The remaining disturbances are static-dynamic, in that the disturbance regimes (e.g., area, frequency) are predetermined, but their impacts (i.e., return green arrows) are dependent on the forest composition and structure. mixed forests of pine and boreal hardwood species, such as aspen and paper birch.
Fire was the dominant disturbance in this region during the presettlement period (1860-1890), affecting about 0.3% of the landscape annually (White and Host 2008). In contemporary times (i.e., by the 1970s), fire has been greatly reduced through fire suppression (Sturtevant et al. 2012), and instead, timber harvesting is the most prevalent disturbance in Minnesota with an annual harvest rate of 0.6%. Forest management is focused on pulp production; the most economically important species in the region are aspen and red pine (Minnesota Department of Natural Resources 2011). Insect defoliation is currently the most extensive natural disturbance-for example, the last outbreak of forest tent caterpillar (Malacosoma disstria Hubner) affected nearly 100% of the entire study area (Appendix S2). However, impacts by defoliation are host-specific and often affect growth more than mortality (Churchill et al. 1964

Simulation modeling
We used a simulation modeling approach because multiple disturbances are inherently difficult to disentangle using remote sensing data or to estimate from historic reference data (e.g., General Land Office). We chose a process-based framework that simulates climatic effects on vegetation and disturbance dynamics, while allowing for emergent behavior under multiple disturbances. We used LANDIS-II v6.1 (Scheller et al. 2007), a forest simulation model, to simulate the emergent effects of climate change, harvesting, wind, fire, and insects on forest succession and carbon and nitrogen cycling. In LANDIS-II, the landscape is comprised of interconnecting grid cells. Each cell is assigned to an ecoregion (wherein climate and soil properties are assumed to be homogenous), and within each cell, trees are represented as species-age cohorts, not individuals (Mladenoff 2004). Forest dynamics are projected by simulating the growth processes of cohort establishment, biomass growth, and competition and the degenerative processes of senescence, mortality, and disturbance. Species-age cohorts are therefore dynamic over time; there may be multiple species and age cohorts within each cell. Disturbances in LANDIS-II are stochastic with species' response to disturbance dictated by life history attributes and competition (Roberts 1996). To simulate succession and disturbance, we used eight extensions for this study: Century Succession, Dynamic Fire, Dynamic Fuels (collectively often referred to as Dynamic Fuels and Fire System), Biological Disturbance Agent, Insect Leaf Biomass, Biomass Harvest, Base Wind, and Linear Wind.
Our map of current tree species composition for the simulated landscape of 3.4 million hectares (Fig. 2) contained 34 tree species (Table 1) at a resolution of 4 ha (Lucash et al. 2017). The map was a compilation of plot data from the CNF and imputed U.S. Forest Service Forest Inventory and Analysis data (FIA; http:// fia.fs.fed.us/) using three maps: two published studies (Wolter et al. 1995, Wilson et al. 2012 and MN State DNR (unpublished data). Details of map assembly are outlined in the appendix of Lucash et al. (2017).

Forest succession and C dynamics
To simulate forest succession, we used the Century Succession extension (v4.0.2) of LANDIS-II , Lucash et al. 2017). This extension simulates aboveground and belowground growth of leaves, wood, fine roots, and coarse roots of each cohort on each site at a monthly timestep . It also simulates C and N cycling (described below) and contains a simple bucket model used to calculate water availability. To calculate growth, it uses algorithms that integrate species-specific life history attributes (e.g., longevity, shade tolerance), cohort age, competition, climate (e.g., air Red elm temperature, precipitation), soils (e.g., texture, drainage), and soil water and nitrogen availability. Species-specific, maximum growth rates are user-defined but can be reduced by soil moisture, available nitrogen, soil temperature, leaf area index, and maximum aboveground biomass. For example, growth limitation by soil moisture is a function of available water (determined through algorithms of precipitation of snow and rain, evapotranspiration, interception, runoff, and leaching; Lucash et al. 2017) and a user-defined curve of how production relates to soil water availability (specified by the extension's predecessor, the CENTURY soil model v 4.5; Parton et al. 1983). Temperature limitations to growth are determined by climate inputs (described in Climate data) and a user-defined soil temperature response curve (Poisson density function; Parton et al. 1983). The Century Extension also simulates tree mortality caused by senescence (continuous loss of trees and branches) and age (which accounts for the higher levels of mortality as a cohort approaches its life expectancy). In addition to growth and mortality, it also simulates the spatial pattern of regeneration via seeds or resprouting (vegetative reproduction) using life history attributes (e.g., age to sexual maturity and seed dispersal distances) and water and light availability (Scheller et al. 2007).
The Century Succession extension also simulates C and N cycling in pools of detritus (foliar, woody, fine roots, and coarse root detritus), soil (fast, slow, and passive soil pools), and vegetation biomass (leaf, wood, fine roots, and coarse roots by species and age; Scheller et al. 2011Scheller et al. , 2012. Decomposition rates are a function of litter characteristics (e.g., leaf C/N ratios and lignin content) and soil conditions (e.g., soil moisture, temperature, and soil texture; Parton et al. 1983). The nitrogen cycle in the Century Succession extension is dynamic and tightly coupled between the atmosphere (wet and dry N deposition), vegetation (N uptake), and soil (N mineralization and leaching; Lucash et al. 2014). Although the extension runs internally at a monthly timestep, model output was produced only every 10 yr. Details of Century Succession parameterization and the calibration and validation of Century are outlined in the appendices of Lucash et al. (2017).

Climate data
Climate data are required for spin-up to calculate the biomass associated with each species-age cohort (in effect growing the trees from their establishment date to their age at the model start time, 2010) and for simulating growth and establishment into the future. For spin-up and the historical climate scenario, we used the University of Idaho meteorological data at a 4 km resolution over the period 1979-2010 from the USGS data portal (http://cida.usgs.gov/gdp/) using areaweighted averages; we randomly sampled from this climate stream on an annual basis. To simulate climate change, we used 12-km projections from the Bias Corrected Constructed Analogs V2 Daily Climate Projections from CMIP5. Based on previous work where we evaluated 45 climate change models/RCPs (Lucash et al. 2017), we chose two of these future climate projections representing a warm (ACCESS RCP 4.5) and a hot (MIROC ESM RCP 8.5) climate change scenario to bracket the range of expected changes in climate.

Wildfire
Wildfire was simulated on an annual basis using the Dynamic Fuels and Fire System extension version 4.0 (hereafter referred to as DFFS) for LANDIS-II (Sturtevant et al. 2009). Dynamic Fuels and Fire System simulates fire in fire regions; fire ignition and spread are assumed to be relatively similar within individual fire regions, though fires can spread from one region to another. Both ignition and spread rates are influenced by fire weather and topography. Fuel types are defined on every site across the landscape and are derived from a combination of forest composition, age, and recent disturbance. Although fuel classifications do not directly influence the probability of ignition on a given site, fuel type parameters do help determine rate of spread (see Appendix S1) from successful ignitions. Crown fraction burned, which is determined by rate of spread, fine foliar moisture content, and fuel-type-specific parameters, is used as an indicator of fire intensity, the energy output of the fire. Fire severity, the effect on the forest (in this case, the number of cohorts killed) is a function of the fire tolerance and age of species cohorts that appear on a given site. Fire severity ranges from 1 (low-severity surface fire) to 5 (high-severity crown fire).
Dynamic Fuels and Fire System algorithms are based on the Canadian Forest Fire Weather Index System, and as such require the calculation of the Fire Weather Index (FWI) and fine fuels moisture code (FFMC). The DFFS extension was modified to internally calculate FWI by using daily temperature, precipitation, wind speed velocity, and wind direction derived from global circulation models (Lucash et al. 2017). We also modified the ignition probability equation in the DFFS extension as follows: where P is the probability of successful ignition, FWI is the Fire Weather Index value for each day, b 0 is an FWI shape parameter, and b 1 is an FWI scale parameter. This equation was derived from Beverly and Wotton (2007), who examined ignition and sustained flaming probabilities based on climate and different fuel types. This equation maintains the observed relationship between FWI and ignition probability and provides a continuous calculation of fire ignition probability, based on the relationship between fire weather and fuel ignition probability. The DFFS was applied using its durationbased option, where fire duration distributions are calibrated to match historical fire size distributions and fire rotation periods (Sturtevant et al. 2009). When applied in this way, the fire regime (frequency, extent, severity) within DFFS is an emergent property of weather conditions, fuels, and topography. Average fire rotation period was~2200 yr, which was about 10% lower than observed (2600 yr) for the study area. Details of our parameterization, calibration, and validation of DFFS are given in Appendix S1.

Insects
Insect outbreaks of jack pine budworm (JPB; Choristoneura pinus Freeman), two-lined chestnut borer (CB; Agrilus bilineatus Weber), and forest tent caterpillar (FTC; Malacosoma disstria Hubner) were modeled using two separate LANDIS-II model extensions. Biological Disturbance Agent (BDA; v 3.0; Sturtevant et al. 2004) was used to model JPB and CB outbreaks. Within BDA, temporal patterns of outbreaks are defined by frequency (i.e., mean time between disturbances) and the type of temporal pattern (i.e., periodic, random, chronic), and insect dispersal distance.
Host specificities within the BDA are species-specific and defined within a lookup table for each insect by tree species and age into three susceptibility classes (high, medium, and low). When outbreaks occur, the BDA extension uses host proportion at the cell and neighborhood levels to calculate the likelihood of a given cell being attacked by a given pest, with further constraints imposed by the user-defined dispersal distance. Cohort mortality is determined by an intensity class based on host resources available at the site. Details of our parameterization can be found in Appendix S2.
Jack pine budworm impacts were limited to jack pine, and its outbreaks were restricted to the western half of the study area where such outbreaks were historically prevalent (J. Almendinger, personal communication; Appendix S2). Its outbreaks were simulated as occurring periodically with mean frequency (i.e., length of time between outbreaks) of 8 yr. By contrast, CB outbreaks occurred when the host species (northern red oak) was stressed by the combination of drought and defoliation by FTC (J. Albers, personal communication). Drought which triggered CB dynamics occurred when the Palmer Drought Severity Index was less than À1.5. While patchy spot formation is typical of CB outbreaks (Michael and Jana Albers, personal communication), tree mortality by this agent is typically mediated by individual tree stress.
Insect Leaf Biomass v2.0 (Foster et al. 2011) is designed specifically to model the effects of defoliators informed by remote sensing analyses and was used to model how FTC affects both growth and mortality. Temporal dynamics of outbreaks are defined by the frequency and duration of periodic outbreaks. Spatial extent of individual outbreaks is specified by user-defined parameters, while the patterns of damage within these spatial extents are based on host density at a site and the previous year's defoliation, dictated by statistical relationships parameterized by multiyear remotely sensed outbreaks (Gustafson et al. 2017). Insect-host impacts are species-specific and, unlike BDA, are defined by cumulative defoliation functions which determine both growth reduction and mortality (Foster 2017). In the case of CB, the temporal pattern of outbreaks was defined by the spatiotemporal overlap of drought and host defoliation by FTC. We parameterized FTC outbreaks to occur with a frequency averaging every 16 yr, where the duration (yr) for a given outbreak was selected from a Weibull distribution (k = 1.2, k = 0.47; Appendix S2). Defoliation patch parameters were calibrated such that the area affected by a moderate-to-severe defoliation during a given event was 10-40% of the study area, based on analyses of aerial survey records from the Minnesota Department of Natural Resources. Moderate defoliation is defined as 34-66% leaf biomass (modeled) or area (mapped) consumed, while severe defoliation is defined as >66% leaves consumed. Host impacts were focused on deciduous tree species, where aspen, basswood, and oaks are the primary hosts.
Of the three insect pest species simulated, only CB was explicitly affected by climate (i.e., drought) and by another disturbance (i.e., FTC defoliation; Fig. 1). All insect pest disturbances were affected by the amount of available host. Both extensions were run at an annual timestep; details of our parameterization, calibration, and validation of all three insect disturbance regimes are given in Appendix S2.

Harvesting
The Biomass Harvest extension v. 2.2.2 (Scheller and Domingo 2015, 2016) was used to simulate forest management activities across the landscape. This extension simulates biomass removal (including partial biomass removal of individual cohorts) and the planting of tree species following harvesting. Harvesting is controlled by prescriptions targeted to specific forest types that specify how much biomass is removed from which species and ages within a forest stand (Gustafson et al. 2000). We calibrated acres harvested for each management area and prescription, based on current practices (USDA 2007, D'Amato et al. 2008; C. R. Blinn, personal communication), with details specified in Lucash et al. (2017). In brief, the harvest rate was 0.05% of the total landscape per decade distributed across three principle harvest prescription types (clear-cut logging, selective harvest, and thinning treatments); observed harvest rate was estimated to be 0.04% per decade. Our simulated interactions between harvest and other disturbance types were mediated by forest composition and age (Fig. 1) and did not include salvage logging.

Wind
The Base Wind Extension v2.1.2 (Scheller and Domingo 2003) was used to simulate small (≥4 ha) and moderate (up to 1000 ha) patches of microburst wind disturbance, while the Linear Wind Extension v 1.0 (Gustafson 2016) was used to simulate large, linear wind events such as derechos and tornados at a decadal timestep. Mortality by wind is age-dependent in both these extensions, and therefore, the oldest cohorts have the highest mortality when disturbed by wind. Both wind extensions were calibrated together under historic climate to match wind event sizes from Frelich (2002) and the wind return interval in White and Host (2008) using the procedures outlined in Lucash et al. (2017). In our simulations, the average size of a wind event was 54 ha, the maximum historical wind event was 2395 ha, and the mean wind return interval was 556 yr; the wind regime did not vary by climate scenario (Lucash et al. 2017). We assumed that climate change did not directly affect the wind disturbance regime, but rather altered forest composition and structure which influenced the spatial and temporal pattern of windthrow (Fig. 1).

Experimental design and treatment effects
Under each of the three climate scenarios, we ran each disturbance separately (i.e., fire, insects, harvest, and wind), referred to as single disturbance regimes, and then, all disturbances were simulated together (termed concurrent). Each climate-disturbance scenario was replicated three times to capture the spatial variation and stochasticity in fire, insects, windstorms, forest harvesting, and regeneration after disturbance.
We used the total net live biomass lost (Mg) across the landscape within a decade (timestep = 10) as an index of mortality. Mortality, as defined here, is likely underestimated, since it is the sum of biomass removed minus growth during the decade. We also estimated the extent of each disturbance, since defoliators like FTC affect a large area but their impacts are typically diffuse. We limited our assessment of disturbance extent by FTC to the area affected by high severity (i.e., >66% defoliation) at the conclusion of a given outbreak to be most consistent with area estimates by the other insect disturbance agents.
We compared mortality between single and concurrent disturbance regimes. Presumably, including more types of disturbances will result in more mortality and more area disturbed. However, we anticipated that mortality and extent disturbed would be less when all disturbances were simulated concurrently than when total mortality and disturbance extent were summed across all simulations when disturbances were simulated separately (i.e., concurrent < [harvest + fire + insects + wind only]), due to compensatory mortality from disturbance interactions. We further expected differences in total mortality between single and concurrent disturbance regimes to influence carbon stocks, quantified in terms of aboveground biomass and soil carbon.
We evaluated the relative importance of disturbances and climate change on forest condition via four metrics: community structure, successional stage, hardwood-to-conifer ratio, and tree species diversity. For community structure, we applied nonmetric multidimensional scaling using the vegan-community ecology package in R (Oksanen et al. 2013; R Development Core Team 2014). Using the approach outlined in Scheller and Mladenoff (2005), we created a species-by-ecoregion matrix of average aboveground live biomass for each ecoregion at the initial conditions (year 2010) and for each climate and disturbance scenario at year 2110. We only used one replicate, as variation in species biomass between replicates within scenarios was relatively small. We also computed Shannon's index to quantify landscape-scale tree species diversity. To approximate differences in successional stages across scenarios, we used the shade tolerance of each species (ranked in the model from 1 = shade intolerant to 5 = shade tolerant), weighted it by its relative species biomass, and then computed average shade tolerance at the landscape level for each timestep.

Interactions between disturbance and climate
The greatest cumulative mortality over the next century was human-caused under historical climate (i.e., no climate change); the tree mortality ranking was harvest > insects > wind > fire when disturbances were run separately (Fig. 3A). Under moderate climate change, the ranking shifted to harvest > wind > insects > fire, and under more extreme climate change, it was wind > harvest > insects > fire. The amount of area disturbed for each disturbance type was either unaffected (wind and fire) or declined (harvesting and insects) with increasing climate warming, with area disturbed by insects > harvest > wind > fire (Fig. 3B). Disturbance by insects-particularly FTC-was widespread, but had far less impact in terms of mortality per ha disturbed, reflecting the comparatively defuse and benign effects of defoliation relative to the other disturbance types.
Climate change increased tree mortality for all the disturbances, particularly fire, which increased by 360% between the historical and hot scenarios (Fig. 3A). However, fire remained a minor disturbance in this landscape (especially in terms of extent) even under climate change. Insectinduced mortality increased by 58%, despite the fact that the area affected either declined (À15% for FTC) or remained constant (i.e., CB and JPB) between the historical and hot scenarios (Fig. 3B). The enhanced mortality under climate change (Fig. 3A) was not caused by increases in extent (Fig. 3B), but rather by increasing biomass. Higher temperatures stimulated growth, regardless of the disturbances simulated, especially in the first 70 yr of the simulation (Fig. 4A), thereby increasing biomass available for mortality by subsequent disturbances (Fig. 4B). However, disturbances eventually caught up; mortality rates under climate change were over twice as high as rates under historical climate (Fig. 4B). Mortality rates were also a larger proportion of the biomass remaining on the landscape (Fig. 4C), confirming that mortality rates were exceeding growth rates. Ultimately, these high mortality rates reduced aboveground biomass under climate change to values equal to or lower than that of the historical climate scenario by 2100. Overall, cumulative mortality was positively correlated with remaining live biomass (P < 0.0006), explaining about 50% of the variation.
Looking more closely at the trends in disturbance extent through time, we found that the extent of harvesting was nearly constant through time under historical climate, but declined after 50 yr under climate change (Fig. 5). Since harvest prescriptions were predefined to target specific tree species and ages (and thus be hypothetically ❖ www.esajournals.org constant), the decline in extent under climate change indicates a change in forest condition that limited availability of forest types targeted for harvesting. These declines in harvesting after 2050 (Fig. 5) did not result in higher remaining biomass (Fig. 4), because growth was severely restricted under the high temperatures. The extent of harvesting was relatively unaffected by disturbances (Fig. 5), but its inclusion in the simulations made a huge impact on aboveground biomass (Fig. 4A). Natural disturbance regimes (wind, fire, and insects), when simulated independently, enabled the accumulation of aboveground biomass approaching or exceeding triple that of the initial conditions, whereas aboveground accumulation in the harvest only scenario was far more limited. We attribute this result to the prevalence of stand-replacing Fig. 3. Comparison of total landscape mortality (A) and extent (B) for four disturbance scenarios (harvesting only, fire only, wind only, and insects only) and three climate scenarios: historic, warm (ACCESS RCP 4.5), and hot (MIROC ESM2 RCP 8.5). Total landscape mortality (Mg) and extent (ha) were summed over the next century. Wind extent was separated into microbursts, derechos, and tornados; insect extent was separated into forest tent caterpillar (FTC), jack pine budworm (JPB), and chestnut borer (CB). Each disturbance was run separately. disturbances within the harvest regime that maintained a generally younger forest less susceptible to natural senescence. Without harvest, aging forests become increasingly susceptible to wind, hence the concentration of mortality within the wind disturbance scenario late in the simulations (Fig. 4B).
Disturbance by forest tent caterpillar, the most widespread insect pest, was far more variable than harvest (Fig. 6). It also interacted more with climate and other disturbances, as the area defoliated by FTC tended to increase when disturbance was generally more prevalent (i.e., concurrent vs. insects alone) and under the historic climate. By comparison, there were no temporal trends apparent in the extent of CB, JPB, wind, or fire disturbance regardless of the climate scenario (data not shown).

Compensatory mortality and carbon stocks
We compared mortality and area disturbed among the sum of the four disturbance regimes simulated singly vs. the total area disturbed when disturbances were simulated concurrently to look for evidence of compensatory mortality (i.e., concurrent < (harvest + fire + insects + wind only)). Under historical climate, we found that disturbances were additive in the first 60 yr of the simulation (Fig. 7A), as the scenarios with concurrent disturbances equaled the sum of all the single disturbances. After 2070, however, disturbances began to offset one another and mortality was lower when all disturbances were simulated concurrently, indicating compensatory mortality among interacting disturbances. This effect was magnified by climate change, with individual disturbances causing a doubling of mortality under the warm climate, and a quadrupling of mortality under the hot climate change scenario.
Although we found evidence that mortality was compensatory among disturbances when simulated concurrently, the extent of disturbance was not. Within individual climate scenarios, the sum of all disturbed area across individually simulated disturbances was similar to the area disturbed when disturbances were simulated concurrently, suggesting a lack of feedbacks between disturbance types with respect to disturbance extent (Fig. 7B). The trend in area disturbed under climate change scenarios was opposite that observed for mortality-the area disturbed began trending upward by the fifth decade under historical climate, but stayed essentially constant under each of the climate change scenarios. Variability in disturbance extent was also high relative to variability in mortality.
Aboveground biomass increased at a slower rate when there were concurrent disturbances operating on the landscape vs. the average across the single disturbance simulations (Fig. 8A). This result simply reflects the fact that greater standing biomass is possible when less area is disturbed (i.e., when disturbances were simulated separately). The effect of climate change on aboveground biomass changed through time, with initial increases in standing biomass over that of historical climate, but indicated declines by the end of the century. Soil C increased linearly across most scenarios, and the effects of disturbance amount and climate change on the rate of increase in soil C were more subtle than observed for aboveground biomass (Fig. 8B). Nonetheless, the higher biomass under single disturbances promoted the accumulation of soil C, while climate change reduced soil C more due to increased heterotrophic respiration.

Effects of disturbance vs climate change on tree species composition
The NMDS ordination of species biomass to discern the relative importance of climate change and disturbance at year 2100 resulted in two primary axes ( Fig. 9; R 2 = 0.95, iterations = 3). The first axis suggested alternative endpoints of succession, with shade-tolerant conifers (black spruce) at one end vs. oaks and northern hardwood species (sugar maple and basswood) on the other. The second axis was defined by earlysuccessional aspen, boxelder, and balsam poplar (Populus balsamifera L.) vs. mixed forests. Scenario groupings suggested NMDS1 was due primarily to within-landscape (i.e., ecoregional) differences, while NMDS2 represents primarily compositional shifts across scenarios, where the greatest separation was between the historical climate and the climate change scenarios (blue vs. yellow/red, Fig. 9A).
Compositional differences among different disturbance scenarios were more subtle than those observed for climate change. Nonetheless, when climate scenarios were analyzed separately, we found most of the separation was along the y-axis (NMDS2) with the wind, fire, and insects appearing more at the top and combined disturbances and harvesting scenario appearing lower (Appendix S3). As with the broader ordination, aspen played a large role in the separation among disturbance scenarios, with more aspen retained under harvesting and concurrent disturbance scenarios vs. the other three single disturbances (insects, fire, and wind).
Climate change was more influential in driving species diversity than the disturbance regime. Species diversity increased under historical climate, but plateaued and then decreased under warming climate (Fig. 10A). Across all climate scenarios, we found that diversity was initially higher when disturbances were run individually, but after~2070, diversity was higher when they were run together. Diversity had a negative relationship with mortality (P < 0.00001), with higher rates of mortality associated with lower diversity. Fig. 7. Comparison of total mortality (A) and total (B) disturbance extent (calculated in ha). Landscape mortality (Mg) was calculated was over the next century under three climate scenarios: historic, warm (ACCESS RCP 4.5), and hot (MIROC ESM2 rcp 8.5). In both graphs, the dashed line represents the sum of all individual disturbance runs (harvesting only + fire only + wind only + insects only), while the solid line represents concurrent disturbances with all the disturbances run together at the same time. We used the sum because we were looking for compensatory effects due to negative feedbacks.
Examination of other dimensions of forest composition did indicate an important role of disturbance. When disturbances were run individually, average shade tolerance was~32% higher than when disturbances where run concurrently, reflecting species responding to the greater light availability under multiple disturbances (Fig. 10B). Climate change also increased shade tolerance, but to a lesser degree (+13%), likely due to higher biomass and greater shading. Further, the disturbance regime had a large influence on whether conifers or hardwoods dominated across the landscape. Under multiple concurrent disturbances, the ratio of conifers to hardwoods decreased relative to the average of the single disturbance scenarios, with comparatively little influence of climate change (Fig. 10C).

DISCUSSION
Our results confirm our hypothesis that climate change exacerbates tree mortality, but not for the reasons we expected. We predicted climate change would increase the extent of disturbance (Dale et al. 2001, Ward andMasters 2007), but instead it had minimal effect on the total extent. Even disturbances that were explicitly linked to climate, such as fire and two-lined chestnut borer, did not noticeably increase in extent under climate change. Instead, we found that rising temperatures increased growth, and these early increases in biomass drove the subsequent changes in climate-induced tree mortality. Higher mortality over the past few decades of rising temperatures and drought has been reported at the global scale (Allen et al. 2010) and highlighted in the southwest (McDowell et al. 2015), but these occurred in concert with increases in the extent of soil drought. Our research indicates that climate change may exacerbate mortality in north-central MN as well, but that it reflects that mortality (when computed as biomass removed) is highly dependent on the temperature-induced increases in growth rates over the next century.
Our results highlight how the relative importance of different disturbances changes as the climate changes. Active forest management was the most important driver of mortality across this landscape, but harvesting became less important under climate change and was overtaken by wind. This is likely because the overall landscape is aging-particularly as a result of the great cutover >100 yr ago and lower harvest rates of secondary forests in recent decades (Niemi et al. 2016), making it more susceptible to wind mortality, as predicted by Rich et al. (2007). Changes in relative importance of disturbances were also driven by shifts in species composition. This was particularly true for aspen, which declined under climate change and caused a negative feedback to harvesting and FTC extent over time, since it is one of the primary hosts of FTC in this region (Parry and Goyer 2004).
As hypothesized, we found evidence for significant compensatory effects among disturbances, which were amplified under climate change. This result was in striking contrast to a global review, which found that disturbance interactions were primarily positive (i.e., disturbances facilitate the activity of other disturbance agents) across a range of disturbances (Seidl  In their study, only 16% of observations illustrated the dampening effect we observed. However, they found that insects and fire were more commonly associated with a dampening effect than other disturbances (closer to 30%), suggesting that their inclusion in a modeling study might be one of the main drivers of the compensatory effects we observed. In our simulations, running each disturbance separately overestimated the amount of mortality that would occur by~60%, compared to allowing all disturbances to co-occur on the landscape. This was exacerbated under climate change (~80%), indicating that it is even more important to account for all simulated disturbances to operate on the landscape concurrently when the climate is changing. In one study, disturbance interactions were ten times more sensitive to temperature than individual disturbances (Seidl and Rammer 2016), which also underscore the importance of considering disturbance interactions as a critical component of disturbance regimes under a changing climate. Further, we found the importance of climate on both mortality and carbon stocks was overestimated when disturbances were run separately, because running the disturbances concurrently dampens the climate change effect. This confirms work across a range of sites and models that illustrate that disturbances either dampen (Loehman et al. 2017) or cancel out productivity gains due to climate change (Reyer et al. 2017). Greater negative feedbacks under concurrent disturbances suggest that the typical approach to calibrating modeled disturbances separately may not be appropriate. Given the predominance of disturbance interactions we found, it is advisable to calibrate all disturbances simultaneously.
Finally, we found conflicting evidence for our hypothesis that the disturbance regime will be more important in shaping forest composition in our landscape than climate change. The disturbance regime had a stronger effect on successional stage (Fig. 10B) and the relative dominance of conifers vs. hardwoods (Fig. 10C), but climate change had a stronger influence on forest community structure (Fig. 9) and species diversity (Fig. 10A). As suggested by previous studies in the region , the importance of boreal tree species (i.e., aspen, balsam fir, spruce species) declined in favor of oaks and northern hardwoods, with the strongest community change related to aspen decline (Fig. 9). We suspect that species diversity was affected by the enhanced competition occurring under climate warming scenarios that may have depressed understory development over time. Indeed, high mortality and the start of biomass declines toward the end of the century (Fig. 4) suggest these processes may have consequences for the resilience of the forest to subsequent disturbance and senescence Scheller 2016, Lucash et al. 2017). Given the climate sensitivity of forest disturbance regimes, the ability to disentangle disturbance and climate may become even more difficult (Seidl et al. 2014, Westerling 2016) as the climate warms in the coming decades (IPCC 2013). Nevertheless, our results highlighting interrelated feedbacks illustrate the importance of capturing all principal drivers affecting forest landscapes, including multiple disturbance types, when projecting climate change effects on forested ecosystems.
Although we tried to capture as much uncertainty as possible by incorporating all the major disturbances (harvesting, wind, insects, and fire) and bracketing the range of climate projections, considerable uncertainty remains unexplored. New insect pests could be introduced (e.g., emerald ash borer) in the region and climate change may directly affect the wind regime over the coming century (Karnauskas et al. 2018), neither of which were explicitly included in our simulations. Also, drought is a significant driver of mortality in many ecosystems (Allen et al. 2010), but this was not explicitly simulated in our study. Our model simulated changes in soil moisture species biomass, C), and ratio of conifers:hardwoods over the next century between the sum of all individual disturbance runs (dashed line; harvesting only, fire only, wind only, and insects only) and all the disturbances run together (solid line) and under three climate scenarios: historic, warm (ACCESS RCP 4.5), and hot (MIROC ESM2 rcp 8.5).
( Fig. 10. Continued) ❖ www.esajournals.org which affected tree growth, but the Century Extension does not simulate tree mortality due to droughts. Simulating drought-induced mortality is also hampered by the fact that CMIP5 scenarios are ensemble means which dampen extreme precipitation events (IPCC 2013). In our study, forest management was simulated using a static set of rules defined by participating stakeholders ), but these guidelines might change as prominent timber species, like aspen, decline under rising temperatures. Also, our study did not consider salvage logging, housing development, or CO 2 fertilization, which could potentially lead to additional interactive effects that alter long-term trends in growth, mortality, and species composition.

CONCLUSIONS
Our study demonstrates that climate change effects on forest landscape dynamics will be mediated by multiple interacting disturbances within boreal and temperate forest systems. Projections of future disturbance rates and forest responses are complicated by the direct and indirect feedbacks among tree species (growth, competition, senescence, etc.), different disturbance types, and climate drivers (i.e., Fig. 1). While our results supported our prediction that climate change would amplify disturbance-related mortality, this result was the consequence of indirect feedbacks between climate and forest processes, rather than the direct effects of climate on disturbance regimes. As expected, projecting mortality and carbon storage implications of multiple disturbance types was not as simple as adding their impacts independently. Similarly, disturbances and climate change each had unique influence on forest landscape composition that, in turn, had important consequences for subsequent harvest and insect disturbance rates. Our results suggest that the effects of climate change on forest-disturbance interactions are as important as considering the effect of rising temperatures or a single disturbance like fire on tree growth and competitiveness.