Understanding and Modeling Forest Disturbance Interactions at the Landscape Level

Disturbances, both natural and anthropogenic, affect the configuration, composition, and function of forested ecosystems. Complex system behaviors emerge from the interactions between disturbance regimes, the vegetation response to those disturbances, and their interplay with multiple drivers (climate, topography, land use, etc.) across spatial and temporal scales. Here, we summarize conceptual advances and empirical approaches to disturbance interaction investigation, and used those insights to evaluate and categorize 146 landscape modeling studies emerging from a systematic review of the literature published since 2010. Recent conceptual advances include formal disaggregation of disturbances into their constituent components, embedding disturbance processes into system dynamics, and clarifying terminology for interaction factors, types, and ecosystem responses. Empirical studies investigating disturbance interactions now span a wide range of approaches, including (most recently) advanced statistical methods applied to an expanding set of spatial and temporal datasets. Concurrent development in spatially-explicit landscape models, informed by these empirical insights, integrate the interactions among natural and anthropogenic disturbances by coupling these processes to account for disturbance stochasticity, disturbance within and across scales, and non-linear landscape responses to climate change. Still, trade-offs between model elegance and complexity remain. We developed an index for the degree of process integration (i.e., balance of static vs. dynamic components) within a given disturbance agent and applied it to the studies from our systematic review. Contemporary model applications in this line of research have applied a wide range process integration, depending on the specific question, but also limited in part by data and knowledge. Non-linear “threshold” behavior and cross-scaled interactions remain a frontier in temperate, boreal, and alpine regions of North America and Europe, while even simplistic studies are lacking from other regions of the globe (e.g., subtropical and tropical biomes). Understanding and planning for uncertainty in system behavior—including disturbance interactions—is paramount at a time of accelerated anthropogenic change. While progress in landscape modeling studies in this area is evident, work remains to increase model transparency and confidence, especially for understudied regions and processes. Moving forward, a multi-dimensional approach is recommended to address the uncertainties of complex human-ecological dynamics.


INTRODUCTION
Reciprocal interactions between disturbances and forested landscapes have been a cornerstone of landscape ecological research and modeling for decades (Baker, 1989;Mladenoff and Baker, 1999;Seidl et al., 2011;Perera et al., 2015). A current frontier in this line of research is the spatially explicit investigation of disturbance interactions across spatial scales (Buma, 2015;Burton et al., 2020;Canelles et al., 2021). Within forested ecosystems, natural disturbances, anthropogenic disturbances, and climatic and anthropogenic drivers interact across a range of spatial scales to shape forested landscapes in term of patterns, processes, and functions ( Figure 1) (Turner, 2010). Modeling disturbance interactions can be particularly challenging when and where the cumulative effects of such interactions, including non-linear and threshold behavior, result in catastrophic mega disturbance (Millar and Stephenson, 2015). Hence, a combination of empirical and modeling studies is needed to understand forest ecosystem dynamics that emerge from the interactions of multiple disturbances as well as biophysical and demographic drivers within forested landscapes (Fraterrigo and Rusak, 2008;Johnstone et al., 2016;Davis et al., 2018).
Natural disturbances (e.g., wildfire and insect outbreaks) have traditionally been investigated separately within different subdisciplines of ecology (fire ecology and entomology, respectively) and further separated from the effects of human disturbances, focusing therefore on the properties of "disturbance regimes" (patch size distributions, severity, frequency or rotation length; Figure 2) and with an emphasis on stochasticity. Despite the apparent stochastic nature of natural disturbances, disturbance regimes generally emerge from feedbacks between internal system processes and external drivers across scales in time and space (Peters et al., 2011). Comparatively, anthropogenic disturbances (e.g., harvesting; Figure 2) have been traditionally viewed as deterministic (i.e., under human control). Yet, human systems are subject to analogous uncertainty and surprise caused by economics, social pressures, and political change that directly impact our ability to implement harvest and other land use plans at the temporal scale of forest rotations (Messier et al., 2019). Interactions among natural disturbances, human disturbances, and vegetation responses to those disturbances further influence system predictability. Proactive and adaptive management practices that embrace system uncertainty are therefore needed to respond to emerging "surprises" in disturbance regimes (Foley, 2005;Peters et al., 2011;Parrott and Meyer, 2012;Allen et al., 2014).
Several conceptual advances have been proposed to help disentangle the emergent properties of disturbance interactions (Foley, 2005;Fraterrigo and Rusak, 2008;Turner, 2010;Peters et al., 2011;Buma, 2015;Messier et al., 2016;Kane et al., 2017;Davis et al., 2018;Ratajczak et al., 2018). These conceptual advances underscore the interplay between pattern and process in natural disturbance dynamics, deterministic and stochastic elements of anthropogenic disturbance, and uncertainty due to climate change that need to be accounted for in modeling frameworks of processes and disturbances across scales (Keane et al., 2015;Urban et al., 2016). Here, we review the conceptual advances and empirical approaches that help disentangle the apparent complexity of disturbance interactions. We further conducted a systematic review of forest landscape simulation modeling studies including more than one disturbance type published since 2010. A class of models known as forest landscape models (FLMs) dominated this field of study. We therefore overview the general design of FLMs, showing how recent developments have shifted from statistically-based disturbance regimes (i.e., static) to process-based methods where disturbance regimes and ecosystem responses are emergent behaviors (i.e., dynamic), and further expanded the ability to choose the degree of system feedbacks depending on the question at hand (Seidl et al., 2011;Keane et al., 2015;Perera et al., 2015). Studies were then categorized according to the specific disturbance interaction questions investigated (as clarified by recent conceptual advances) and the relative balance of static to dynamic model components across disturbance types. We argue for minimum standards in documentation -particularly as model complexity increases -for increased transparency and confidence in model results. We conclude highlighting the current limitations, frontiers, and directions in the understanding and modeling of disturbance interactions at landscape levels.

CONCEPTUAL ADVANCES
Disturbances act upon the components of an ecosystem in a way that changes the structure defining the system (Pickett and White, 1985;Rykiel, 1985;Lindenmayer et al., 2017;Newman, 2019). Within forested ecosystems, a disturbance typically disrupts the functioning of its dominant life form (trees) via physical or chemical damage impacting growth and survival. The effects of disturbance can range from "pulse" disturbances that are concentrated in space and time and lead to abrupt change (Jentsch and White, 2019) to "press" disturbances that are diffuse in space and time and lead to cumulative system stress (Bender et al., 1984). The intensity of a disturbance is measured in terms of force, energy, or analogous quantity (e.g., density of insects), while the severity of the disturbance is the consequence FIGURE 1 | Disturbances at the mesoscale (landscape-level) are affected processes and drivers (bottom-up) at the microscale (stand-level) and higher-order drivers (top-down) at the macroscale (biome-level). Landscape disturbances interact in space and time with reciprocal feedback evident at the landscape level (double-arrow dotted lines). Cross-scale interactions occur where either both bottom-up processes and/or top-down drivers amplify or attenuate disturbance processes within landscapes via threshold behavior in time and space.
of that disturbance intensity on the state of the system (Keeley, 2009). Diversity in life history traits among tree species and other life forms mediate the relationship between disturbance intensity and severity via adaptation (Noble and Slatyer, 1980). Defining characteristics of a disturbance also depend on the scale of observation (Allen and Hoekstra, 2015). For example, a disturbance may be locally abrupt and severe (e.g., killing individual trees) but broadly diffuse and mild if killed trees are thinly dispersed across a large area. Once disturbed, forest system traits such as crown closure, height, and composition can take a long time to recover via the processes of recolonization, growth, and succession, respectively. Disturbances therefore affect landscape spatial heterogeneity including tree-species composition, age structure, and configuration (e.g., James et al., 2011b;Sturtevant et al., 2014).

Disturbance Interaction Process and Terminology
The opportunity for disturbance interaction occurs when one event follows another. Kane et al. (2017) observed that the likelihood of such occurrence depends on both the frequency and size of each disturbance, where opportunities for disturbance interactions increase by chance alone as their respective frequency and size increase. The nature of the specific interaction may be determined by a range of factors such as the disturbance mechanisms at play, the precise sequence, extent, and arrangement of the overlap, and the time since the previous disturbance (Kane et al., 2017; Figure 3). Interactions may include one or more types defined by the constituent components of the disturbance eventits incidence/extent, its intensity, and (or) its effects (Kane et al., 2017). Any of these components may be positively or negatively affected, or otherwise unaffected by the previous event. The nature of the interaction can also take multiple forms-for example the effect of one disturbance on another, or the combined effect of two disturbances on an ecosystem property or function.
From a disturbance-interaction standpoint, the legacy from one disturbance becomes the new system state for the subsequent disturbance (Peterson, 2002). A "linked" disturbance interaction occurs when the legacy of one disturbance is a critical driver for another, affecting the likelihood, extent, or severity of the second disturbance (Buma, 2015). As elaborated by Kane et al. (2017), the specific components underlying the interaction are important to the nature of that interaction. For example, logging disturbance does not generally affect either the incidence (frequency and extent) or intensity of windstorms. In this regard, a wind event following logging disturbance might be considered a "random co-occurrence" (sensu Burton et al., 2020). However, the severity (i.e., effects) of a given wind event can be strongly linked to logging practices -for example, taller trees (generally older) tend to be more susceptible to wind damage than shorter (generally younger) trees, hard edges created by specific logging practices can increase the susceptibility of trees along that edge, and tree species with different rooting habits or wood densities can have differing responses to winds of a given intensity (Quine et al., 1995). Hence, patterns of wind effects may well be linked to the legacy of logging practices. Conversely, if salvage logging is a commonly applied practice, then the incidence and extent of salvage logging will be linked to the disturbance patterns created by wind events.
The occurrence, relative strength, direction, and duration of linked disturbance interactions are mediated by both internal ecosystem processes and external drivers. For example, insect disturbance leading to tree mortality can create pulses of heavy fuels that can set the stage for higher severity fire once the new fuels become combustible, but the period over which that legacy persists as a factor affecting fire depends on fuel decomposition rates (Fleming et al., 2002; but see Section "Empirical Approaches, Challenging Long-standing Dogma"). In some boreal systems, harvesting can transition highly flammable conifer forests to fire-resistant deciduous forest (Carleton and Maclellan, 1994;Cumming, 2001)-a legacy mediated by internal ecosystem recovery processes. Opportunities for disturbances to interact at a given location are therefore sensitive to the timing and sequence of when and where a given disturbance event overlaps a subsequent one, contingent on any relevant timelag effects (Burton et al., 2020). For example, while tree-killing disturbances can influence soil stability underlying the likelihood of landslide disturbances, this linked disturbance interaction is constrained by topography (i.e., slope angle and position) (Carabella et al., 2019).
Certain disturbances are more sensitive to forest composition or age structure than others, and this sensitivity will determine the relative strength of the feedback between vegetation, the specific disturbance, and (by extension) linked interactions with other disturbances. For instance, most outbreaking insects are limited to a few host tree species, leading to relatively strong feedbacks between vegetation and insect disturbance impacts across all three interaction types (incidence, intensity, and effects; Figure 3) (Ohmart, 1989;Hennigar et al., 2008;Bentz, 2019). When landscape spatial pattern is shaped by the cumulative effects of one disturbance (e.g., logging; Sturtevant et al., 2014), the consequent spatial pattern can either enhance or limit subsequent disturbances in terms of their intensity, duration, and size (e.g., insect outbreak; Robert et al., 2018). Such disturbance "legacies" often accumulate across space and time, producing "landscape memory" (Peterson, 2002) with important consequences for future disturbance regimes and consequent landscape dynamics (Foster et al., 1998;Buma and Wessman, 2011).
"Looped disturbances" form a subset of linked disturbance interactions in the form of recurrent disturbances with feedbacks that create a metastability of the system (Burton et al., 2020). Classic looped disturbances are often associated as a component of disturbance regimes supporting a systems' historic range of variability (HRV; Landres et al., 1999). For example, boreal conifers with fire-dependent life-history traits (e.g., serotiny) are both adapted and conducive to infrequent but extensive crown-fire events. Likewise, many temperate pine systems once dominant across much of the southern United States (Hanberry et al., 2020) are adapted (e.g., thick bark, high canopy) and conducive (e.g., grassy flammable understory) to frequent surface fire events (Figure 4). Looped disturbances of the same FIGURE 3 | The factors and processes specific to disturbance interactions. A series of interaction factors define the degree and strength of the disturbance interaction, which may manifest through one or more interaction types any of which may be positive, negative, or neutral. Feedback between interaction types (and also drivers; Figures 1, 5) may manifest as cross-scale interactions. Ultimately, the effects of the disturbance interaction determine the ecosystem response ( Figure 5). Reprinted from Publication: Forest Ecology and Management Vol 405, Characterizing interactions between fire and other disturbances and their impacts on tree mortality in western U.S., p188-199, 2017, with permission from Elsevier. type can occur wherever life history traits of the primary vegetation are adapted to given type of disturbance so that it replicates itself -such as the build-up of fir regeneration in the understory that is released by the periodic defoliation by eastern spruce budworm (Choristoneura fumiferana Clemens) in the North American boreal forest (Baskerville, 1975). More complex looped disturbances across disturbance types are also possible when those types help reinforce each other in sequence (as illustrated by Burton et al., 2020).
"Compound" disturbance interactions occur when the combination of two or more disturbances have a multiplicative effect over the state of the system (Buma and Wessman, 2011), often via the recovery processes affecting ecosystem resilience (Buma, 2015;Coop et al., 2020). In the above example for boreal conifers, two fires in rapid succession can burn the tree regeneration of a serotinous species before they have sufficient time to produce new seeds-dramatically delaying subsequent system recovery (Figure 4). Indeed, starkly different vegetation communities resulted from different combinations and sequences of fire and wind in Minnesota as the life history strategies of different tree species were affected quite differently by each disturbance type and their compound disturbance effects (Johnstone et al., 2016). By contrast, regeneration failure due to more frequent fire is a feature of pine savannah (Figure 4) (Hanberry et al., 2020). Decreasing fire frequency enhances regeneration success, leading to increased stem density and canopy connectivity that make the systems more susceptible to catastrophic disturbances such as bark beetles, crown fire, or both in combination (Negrón and Fettig, 2014;Clarke et al., 2016). The consequence of the restriction of a fundamental structuring agent in this example (frequent surface fires) is therefore a bifurcation of the system to the extremes: either closed canopy forest or grassland (Figure 4). Analogous patterns have been observed in many regions (e.g., western North America; Coop et al., 2020). The potential for compound disturbance is likely increasing as climate change can directly and indirectly affect disturbance processes (Soranno et al., 2014;Becknell et al., 2015), human land use intensification is disrupting natural disturbance regimes on a global scale (Foley, 2005), and increasing mobility and transportation increases novel introductions of pests and disease (Liebhold et al., 2017). Such broad-scaled changes can have non-linear "cascading effects" resulting in disturbance regime shifts and transforming the system to alternative states (Buma and Wessman, 2011;Buma, 2015;Johnstone et al., 2016;Ramsfield et al., 2016;Burton et al., 2020).
The effects of the disturbances and their interactions depends on a larger context of system dynamics. This larger context includes (1) external drivers operating at regional scales [e.g., biome constraints (Burton et al., 2020), anthropogenic pressure (Turner, 2010)] (Figure 1), (2) states of the system (e.g., vegetation composition, age structure, soil conditions, etc.) that determine the behavior and effects of (3) disturbance mechanisms (e.g., defoliation, stem breakage, and combustion) that in turn determine the new state of the system (Peters et al., 2011;Figure 5). The new state is considered the legacy of the disturbance, upon which internal ecosystem processes (e.g., recolonization, growth, and succession) enable subsequent system recovery (Peterson, 2002, Johnstone et al., 2016. The type of ecosystem response, as defined by Kane et al. (2017; Figure 3), depends upon how the disturbance legacy relates to the internal processes of the system. If the essential structure (i.e., tree density, size, and composition) remains effectively intact, the ecosystem response is one of resistance. If the system requires a longer recovery period but remains within the natural variability of the broader ecosystem, the ecosystem response to the disturbance interaction is one of resilience. If the recovering system falls outside that natural variability the disturbances compounded to result in an alternate system state ( Figure 5). External drivers such as climate and climate change interact with disturbances via three pathways  -directly via its effect on a FIGURE 4 | Compound disturbances are those affecting system recovery (i.e., resilience), where the Historic Range of Variability (HRV) serves as the context for ecosystem resilience. By example, serotinous boreal conifers are adapted to high intensity crown fire but require sufficient time to develop a canopy seed bank. Non-serotinous pines often require frequent fire to dampen regeneration to maintain open savannah structure. Fire frequencies intermediate to these extremes can consume the seed producing potential in either system, leading to a conversion to new systems (see text for references).
disturbance property such as incidence or intensity, indirectly via effects on internal processes that modify system properties that in turn affect the disturbance process, and through a secondary interaction with another disturbance process that in turn affects the focal disturbance (i.e., the focus of this paper).

Cross-Scale Interactions
To understand processes affecting multiple disturbances and their interactions, integration of processes and disturbances operating at more than one scale is needed (O'Neill et al., 1996;Wu and David, 2002;Peters et al., 2007;Soranno et al., 2014) (Figure 1). Living and non-living elements of a system tend to interact most strongly with other elements that have characteristic process rates and spatial scales similar to their own (Allen and Hoekstra, 2015), leading to more closely-linked feedbacks within a single domain of scale (Messier et al., 2016). For example, at tree and stand-levels, life history traits of tree species affect demographic (mortality and recruitment), ecosystem (nutrient cycling), and disturbance (e.g., flammability, plant defense, etc.) processes that can include neighborhood effects that may also be self-reinforcing (Frelich and Reich, 1999). Nonetheless, the statistical properties of disturbance regimes (i.e., Figure 2) are not defined by processes at any one scale, but rather they emerge from the interactions of processes across scales (Figures 1, 5). So while compartmentalizing system dynamics into scale-specific components is a useful way to disentangle system complexity (Allen and Hoekstra, 2015), cumulative changes at one scale can either amplify or attenuate system behaviors across scales-a phenomenon known as a "cross-scale interaction" (Peters et al., 2004(Peters et al., , 2007 Figure 3). For example, when local fuels become interconnected across landscapes to the point that firestorms (fire events that create their own weather) are possible, then discontinuities that previously served as fuel breaks may become irrelevant (Peters et al., 2004).
FIGURE 5 | Disaggregating disturbance mechanisms, drivers, and internal processes affecting ecosystem response to disturbance interactions. Drivers (yellow) refer to the microscale and macroscale factors affecting mesoscale dynamics in Figure 1. Hourglass symbols represent system changes that can either be abrupt (redrapid disturbance effects), gradual in time (blue -succession), or gradual in space and time (black -natural variability in local system states). Depending on the degree of change in system properties, the legacies that remain, and the recovery processes, the mesoscale ecosystem response can be one of resistance, resilience, or (if overlapping disturbances compound to impede recovery processes) transformation to an alternative system. (Adapted from Peters et al., 2011). Cross-scale interactions often underlie dramatic change as a system shifts from a strong negative effect to a strong positive one. More precisely, different processes operating at different scales can work synergistically to amplify disturbances across scales. For example, tree-killing bark beetles subject to strong, non-linear and density-dependent feedbacks (e.g., Raffa et al., 2008) are often linked to other disturbance types. Under lowdensity endemic situations, individual healthy trees easily defend against attacking beetles. Outbreaks are generally triggered by drivers that enable beetle populations to grow sufficiently to overwhelm the resources of healthy trees. In Europe, windstorm events provide a pulse of food resources in the form of weakened or fresh-killed trees (i.e., a disturbance legacy) that enables European spruce bark beetle (Ips typographus L.) populations to increase to the point they can overcome adjacent healthy tree defenses, such that a strong negative feedback (strong defense) changes to a strong positive feedback (ample food) (Kausrud et al., 2012). Such an effect can be further amplified if wind events overlap periods of drought, where trees are already stressed and defenses weakened (Kausrud et al., 2012). Climate warming in central Europe serves as a common driver increasing the frequency and intensity of windstorms, the frequency and intensity of drought events, and decreasing generation time of I. typographus to dramatically increase population growth potential (Senf and Seidl, 2018, Figure 6). Climate-driven synergies underlying explosive outbreaks of this insect have been further exacerbated by a long history of spruce plantation establishment in the region (Jansen et al., 2017), leading to relatively continuous food resources (Seidl et al., 2016b, Figure 6). While this example appears like a "perfect storm" event, it is a pattern replicated across many forest insect systems across the globe, particularly in bark beetle systems subject to similarly strong non-linear feedback processes (Burton et al., 2020;Kneeshaw et al., 2021).

Insights From Conceptual Advances
Current conceptual advances in disturbance interaction investigation include formal disaggregation of disturbances into their constituent components, embedding the disturbance process into system dynamics affected by both external drivers and internal processes (Figure 5), and explicitly recognizing the interaction factors, types, and ecosystem responses (Figure 3). While early modeling investigations of disturbance interactions focused only simple binary interactions (Figure 7; Mladenoff et al., 1996), the desire to include more processes structuring landscapes has spurred development and application of more complicated and complex models including many interacting disturbances having direct and indirect effects (Figure 7) -a theme we will return to in Section "Modeling Approaches." Furthermore, under the Anthropocene, socio-economic drivers are crucial to include as well in forest disturbance models (e.g., Figure 6). Importantly, the same drivers affecting disturbance mechanisms may also influence the internal processes affecting ecosystem response ( Figure 5). The sum total of the drivers' effects on disturbances, internal processes, and their interactions can lead to changes in system properties and system behaviors that can be either gradual or abrupt, depending on the direction and form of the effect (Ratajczak et al., 2018). Empirical disturbance interaction studies (see Section "Empirical Approaches") should therefore take care to define the nature and form of the interaction explicitly, while models (see Section "Modeling Approaches") should effectively capture the form, FIGURE 6 | Direct and indirect interactions of disturbances and drivers contributing to the amplification of a disturbance agent, illustrated using the European spruce beetle example from central Europe (see Section "Cross-Scale Interactions").
direction, and magnitude of interaction processes, including accounting for processes that transcend spatiotemporal scales.

EMPIRICAL APPROACHES
In the last decade, there have been several excellent papers reviewing interactions among disturbances. For example, Kleinman et al. (2019) reviewed compound disturbance interactions reported in nearly 300 journal articles, finding that such studies tended to be concentrated in North America, and focused primarily on fire, wind, and salvage logging. Salvage logging interactions have been more extensively reviewed by Leverkus et al. (2018) and Thorn et al. (2018), demonstrating a concentration of such studies in North America and Europe. Kane et al. (2017) reviewed fire interactions with other disturbance types concentrating on studies from the western United States (Figure 3). Kolb et al. (2016) reviewed the responses of biotic disturbances (insects, disease, and parasitic plants) to drought in the United States. Emerging from their review was a series of observed trendlines between different classes of insects and disease and their impact response to drought severity. Canelles et al.'s (2021) systematic review focused on disturbance interactions between insect disturbance and other forest disturbance types, and included a synthesis revealing different processes underlying those interactions.  performed a meta-analysis of the direct, indirect, and interaction effects of climate change on disturbance processes. Among the commonalities of these different synthesis articles is the geographic concentration of studies primarily in North America and Europe, an emphasis on certain disturbance types, with important knowledge gaps in mass-movement disturbances such as avalanches and landslides, meteorological events beyond wind and drought (e.g., ice storms), and disease. The reviews further noted an emphasis on the amplifying effects of disturbance interactions, with fewer empirical examples of the buffering effects of such interactions. While most of these reviews contained landscape modeling studies, they were not the primary focus of the research. Here, we present different empirical approaches to investigating disturbance interactions as a fundamental source of knowledge underlying the conceptualization and implementation of disturbance interactions within landscape models ( Table 1).

Serendipitous Studies
The study of disturbance interactions in forested ecosystems has unique challenges because it is nearly impossible to manipulate for the overlap of two (or more) disturbances in the field. Consequently, there are many serendipitous correlational studies in the disturbance interaction literature. For example, the extensive bark beetle epidemics affecting western North America in the last two decades, coincident with both widespread drought , and highly active fire seasons  enabled a broad range of opportunistic FIGURE 7 | Challenges and approaches to simulating multiple disturbance interactions. When disturbance interactions are simulated as direct (A), the number of interactions can multiply quickly to become increasing complicated, especially if the interactions are bidirectional (i.e., double-arrows). Indirect interactions are those mediated through ecosystem state variables, such as composition or fuel (B). Depending on model design, it may be more elegant (i.e., fewer interactions) to explicitly model the indirect processes affecting disturbance processes, but such models may also be more interdependent and therefore more complex.
disturbance-interaction studies (Kulakowski and Veblen, 2007;Kulakowski et al., 2012;Stevens-Rumann et al., 2015;Agne et al., 2016;Millar and Delany, 2019), including the widespread study of salvage logging and its effect on ecosystem response and recovery (Leverkus et al., 2018). Such studies are often descriptive in nature, documenting patterns rather than understanding the processes producing them. Nonetheless, as the Anthropocene increasingly enhances disturbances under climate change, such disturbance overlap may become increasingly frequent (Seidl et al., 2016b;Burton et al., 2020;Kim et al., 2020).

Retrogressive and Prospective Studies
One common approach for investigating disturbance interactions is "retrogressive" (sensu Simard et al., 2011), where investigators reconstruct disturbances over a longer period of time, using a variety of data sources such as archived satellite imagery, air photos, dendrochronology, and disturbance records (e.g., Bebi et al., 2003;Howe and Baker, 2003;DeRose and Long, 2012;Sturtevant et al., 2014;Hart et al., 2015;Meigs et al., 2015;Robert et al., 2018Robert et al., , 2020. A challenge facing retrogressive approaches is that the quality of data defining underlying drivers, such as forest composition or age structure, generally degrades the further back in time one investigates (McKenzie and Robert et al., 2020) (Table 1). Conversely, empirical data defining a given disturbance legacy (e.g., fuel) are often used to drive models to infer future disturbance behavior (e.g., flame length) (DeRose and Long, 2009;Simard et al., 2011;Meigs et al., 2015;Seidl et al., 2016a). Such "prospective" studies, including those referenced above, typically apply models at the plot or stand scale as their method to infer consequent system behavior.

Experimental
A number of recent manipulative experiments have focused on disturbance mechanisms and their interactions on a variety of ecosystem responses. For example, Royo et al. (2010) examined the interactive effects of canopy gaps, ungulate browsing pressure, and fire via girdling, fencing, and prescribed fire, respectively. The authors found substantial synergistic effects leading to greater herbaceous diversity under the combined disturbances relative to any one disturbance alone. Quentin et al. (2012) applied a combination of irrigation and mechanical defoliation to investigate the interactive effects of water supply and leaf consumption on eucalyptus (Eucalyptus globulus Labill) plantations in Tasmania. Sparks et al. (2018) investigated the interactive effects of drought stress and fire on western larch (Larix occidentalis Nutt.) seedlings and found that severely drought-stressed seedlings had higher survival following fire relative to moderately drought-stressed seedlings, in part due to the physiological response to drought. Cannon et al. (2014) simulated windthrow disturbance by the static winching of Retrospective Disturbances reconstructed over longer time frames from either archived spatial data (remote sensing, air photos) or disturbance surrogates (tree-ring data).
Can more directly address temporal factors, may or may not contain a wide breadth of severities and ecosystem states.
Data quality and quantity degrade as one goes back in time; Inference of causal relationships not possible. Bebi et al., 2003;Howe and Baker, 2003;DeRose and Long, 2012;Sturtevant et al., 2014;Hart et al., 2015;Meigs et al., 2015;Robert et al., 2018Robert et al., , 2020 Prospective Empirical data used to parameterize models of disturbance mechanisms to make DI inferences.
Data quality and quantity can be high (present day); enables insight into underlying process.
Availability of supporting process-based models can be uneven among disturbance types; aspatial.
DeRose and Long, 2009;Simard et al., 2011;Meigs et al., 2015;Seidl et al., 2016a Experimental Controlled manipulative studies that focus on mechanisms underlying DI and ecosystem response.
Control for external factors to get at causation; replication; repeatability.
Limited scope, scale, and tree size/life stage. Royo et al., 2010;Quentin et al., 2012;Sparks et al., 2018;Cannon et al., 2014;Westlake et al., 2020 Holistic syntheses Combine lines of evidence from multiple data sources, scales, and models to holistically infer underlying process and consequence of DI.
Comprehensive; big picture enables detection of threshold behavior and cross-scaled interactions.
Case studies may be limited in scope to a system and/or location; generalizations are qualitative. Allen, 2007;Raffa et al., 2008;Anderegg et al., 2015;Soranno et al., 2014;Cobb and Metz, 2017;Ratajczak et al., 2018 Advanced statistics Advanced statistics such as Bayesian hierarchical modeling as well as data-mining techniques (e.g., machine learning) applied to large and complex data sets (remote sensing and time series).
Scope, scale, and range of severities and system states can be very broad; detection of patterns (including threshold behavior) not possible through traditional methods.
Risk of overfitting, patterns may or may not correspond with underlying process; "black-box" data mining may limit interpretation. Fleming et al., 2002;Sturtevant et al., 2014;Seidl et al., 2016b;James et al., 2017;Mezei et al., 2017;Candau et al., 2018;Itter et al., 2019;Robert et al., 2018;Brice et al., 2020;Robert et al., 2020 trees to investigate consequent fire behavior, burn intensities, and consequent effects of those intensities related to elevated fuels on vegetation composition. Such studies can support the development of more process-based models enabling the simulation of disturbance interactions based on first principles (Figures 3, 5). The challenges in their application relate to their inherent limitations, such as the extrapolation of disturbance or stressor responses across tree life stages (Máliš et al., 2016;Dayrell et al., 2018), or the scale or intensity of experimental vs. natural disturbances (Cannon et al., 2014) (Table 1). Indeed, it is often the spatial dynamics of disturbances and their interactions that simulation models are designed to investigate, precisely because empirical measurement of such processes are difficult to quantify explicitly in space. Nonetheless, experimental studies can inform even spatial processes, such as the spatial concentration of ungulate herbivory associated with burned locations, and the cascading effects on plant species diversity across space (Westlake et al., 2020).

Holistic Syntheses
Despite the value of process identification and quantification, a whole systems approach is often necessary to fully understand the nature and behavior of disturbance interactions. For example, Anderegg et al. (2015) offered a data-driven conceptual model to capture the processes underlying interactions between temperature, precipitation, tree stress, and the different physiological, chemical, and ecological constraints influencing herbivory by bark beetles and defoliators that can lead to widespread mortality under climate change. Likewise, Cobb and Metz (2017) adapted the heuristic "disease triangle", traditionally applied to forest and plant pathology, to understand the landscape ecology of tree diseases. The disease triangle consists of the pathogen, its host, and the environment conducive to the infection of the host by the pathogen. From the perspective of disturbance interactions, past disturbance events from various agents might contribute to the "environment" point of the triangle that can either facilitate or impede infection rates. Such syntheses can serve as the foundation for process-based landscape simulation models. Case study syntheses further enable the investigation of disturbance interactions from a more holistic systems perspective. Evidence for the relative strength of specific interactions are built from diverse data sources-often at widely different scales-and gaps in knowledge may be addressed through a combination of prospective modeling and speculation based on systems and scaling theory (Peters et al., 2007). For example, Allen (2007) investigated cross-scaled interactions among drought, bark beetles, fire, and erosion in northern New Mexico based on a conceptual model that included non-linear dynamics that could either amplify or ameliorate the propagation of the different disturbances at play. Synthesized data sources spanned seven distinct studies ranging in spatiotemporal scale from weekly dendrometer (i.e., physiology) measurements to regional maps of inter-annual forest dieback developed via remote sensing. Such syntheses may be uniquely poised to address the sharp non-linearities associated with cross-scaled disturbance interactions that can lead to fundamental disturbance regime shifts. However, broader inference from such case studies will remain limited until the greater body of underlying theory can be rigorously tested (Allen et al., 2014) (Table 1).

Advanced Statistics
Analysis of ecological data in recent decades has shifted from the use of analytically tractable statistical methods constrained by parametric assumptions, such as ordinary least-squared regression and analysis of variance, to more advanced but computationally intensive methods such as mixed effects models, Bayesian inference, and machine learning techniques (Touchon and McCoy, 2016;Holmes and Huber, 2018). Such recent methods are better able to address the complexities of ecological datasets that are rapidly accumulating, including data collected via remote sensing and other automated sensors, as well as time-series data from long-term monitoring and experiments.
The combination of advanced statistical techniques and rich data resources enabled the modeling of underlying processes, legacy interactions, and driver roles related to disturbance interactions in time and space. For example, For example, Mezei et al. (2017) examined temporal patterns in disturbance rates of the European spruce beetle as a function of wind disturbance, prior beetle disturbance, and climatic variables affecting brood frequency. The authors used generalized additive models to allow detection of non-linear relationships, and applied multimodel inference (Burnham et al., 2011) to select from a set of candidate models reflecting plausible alternative hypotheses of underlying processes. Hierarchical and hidden process modeling better accommodate the variability and uncertainties inherent in ecological research (Gimenez et al., 2014), as well as nonlinear responses due to the interactions among disturbances and climatic drivers (Brice et al., 2020), leading to more complete understanding of complex phenomena like disturbance interactions. As such, Itter et al. (2019) applied Bayesian hierarchical modeling to tree ring data in eastern and western regions of boreal Canada to evaluate the separate and combined impacts of tree defoliation and drought on tree growth patterns. They found a synergistic effect of defoliation plus drought on the "ecological memory" of the disturbances in terms of their temporal recovery period. Thus, advanced statistical modeling techniques can be applied to complex datasets to determine the relative importance of disturbance types and factors that affect the state of forest ecosystems.
Similarly, machine learning techniques are often more appropriate to modeling complex and non-linear relationships (Olden et al., 2008). For example, Candau et al. (2018) applied random forest methods to define spatial domains in central Canada where fire and eastern spruce budworm defoliation were likely to interact, based on compositional (i.e., host) and climatic gradients (i.e., moisture). Deep learning methodologies such as convolutional neural networks have been applied to extract useful patterns from complex data, such as imagery and LiDAR, although their application has focused on extracting patterns rather than understanding the processes underlying those patterns (Brodrick et al., 2019). Machine learning techniques are both diverse and evolving, and do have limitations such as propensity to over fit models and (particularly in the case of neural nets) lower transparency in the model-building process (Olden et al., 2008;Liu et al., 2018). Nonetheless, such methods may detect relationships or domains of scale that are poorly understood and worth more explicit attention.

Challenging Long-Standing Dogma
Regardless of the empirical approach, a critical role of empirical studies is to challenge long-standing assumptions of disturbance interactions that often find their way into landscape models. For example, the persistent dogma of "insect disturbance increases fuel load and therefore fire risk"-has been confronted with conflicting empirical evidence. Analysis of spruce eastern budworm outbreaks suggested that they increased forest fire probability in the first 20 years post-outbreak (Fleming et al., 2002), while Meigs et al. (2015) found that fire probability was lower after a western spruce budworm (Choristoneura freemani Razowski) event. These two studies were similar in that they were retrospective studies that evaluated whether the area burned increased during periods of budworm outbreaks, using a range of lag periods (i.e., time since outbreak). By contrast, James et al. (2017) investigated a similar question for spruce budworm in a similar geographic area to that of Fleming et al. (2002) but included weather drivers as covariates in the analyses. They found that the magnitude of defoliation effects was lower than those associated with weather. Likewise, a literature review of bark beetle-fire interactions in western North America concluded the effects of mountain pine beetle (Dendroctonus ponderosae Hopkins) disturbance on fire occurrence and severity were ambivalent (Simard et al., 2008). A subsequent study applied prospective modeling to focus on a particular fire process (i.e., probability of active crown fire; Simard et al., 2011). The authors considered the changes in fuel characteristics with time since disturbance (i.e., dead needles in canopy, deposited needles, fallen trees, and understory development), weather characteristics affecting fire behavior (wind speed and moisture) and interactions between the disturbance legacy and weather drivers (e.g., micrometeorology below canopy). As in Simard et al. (2011), James et al. (2017 found that the potential amplifying effects of the fuel disturbance legacy from the outbreak were constrained by fire weather conditions. Hicke et al. (2012) found that categorizing studies according to a conceptual model of expected lags in different fuel components as they related to different stages of bark beetle outbreaks led to more consistency in results and consequent fire implications across studies.

Insights From Empirical Studies
Empirical study of processes and drivers governing disturbance interactions are prerequisites of plausible models that simulate their dynamics. The different empirical approaches to the investigation of such interactions each have their role to play in the development and application of such models. Empirical study should justify when and where the added complexity of disturbance interactions is warranted, and confront persistent dogma based on perception rather than data. While interest in this line of research has increased in recent decades, disturbance interactions remain complex and challenging to quantify. Recent conceptual advances (see Section "Conceptual Advances") are helping to refine the questions and terminology to channel empirical research, while advances in statistical methods, coupled with increasingly detailed data sources, provide additional tools to probe the characteristically messy data for useful insights. As evidenced by a series of recent reviews, research in this area is active but also unevenly distributed across geography, biomes, and disturbance types. By extension, such deficiencies limit the scope of landscape modeling in this line of research.

MODELING APPROACHES
The application of landscape models with multiple disturbances has exploded with the availability of advanced landscape models and increasing accessibility to computing power. We systematically reviewed the literature using the online search engine ISI Web of Science. Our initial review (access date August 15, 2020) applied the keywords "forest * ", "landscape * , " "simulation * , " and "disturbance * "; where we searched abstracts for landscape simulation studies including least two simultaneous disturbance agents. Landscape simulation studies were identified as those falling within the range of spatial scales falling within the landscape to regional extent (i.e., 5.0 × 10 1 − 5.0 × 10 7 ha) defined by , with model designs that allowed for spatial interactions among finer-scaled entities (i.e., grid cells or polygons). We further limited our review to papers from 2010 forward to focus on the most recent studies with the latest technology. Our initial search yielded 81 papers, however, we found it missed a number of relevant studies. We therefore repeated the literature review of the same search engine, using the names of the modeling software identified in the first review as keywords. This second review, with a final access date of June 06, 2021, identified an additional 66 papers meeting the above criteria (total = 146). The full details of the systematic review and the resulting citation list are provided in Supplementary Materials.
We note that our review methodology had consequences for the final set of papers included in the review. Landscape simulation studies that included more than one disturbance agent were included regardless of whether the investigation of disturbance interactions was stated as an explicit objective. In reality, investigators include multiple disturbance types for a wide diversity of questions that fall along a spectrum of explicit to implicit attention to the interactions of those disturbances within the simulations. Conversely, we did not include studies focused on a single disturbance agent, such as the fire frequency interactions illustrated in Figure 4. Doing so would retain all landscape simulation studies that simulated single disturbance regimes. However, we did retain studies investigating interactions between different agents of the same type, for example, between multiple insect species, or between wildfire and prescribed burning practices. Further, certain classes of models were excluded that have relevance to the broader question of disturbance interactions, but were not landscape models per se. Prospective modeling studies (Table 1) of disturbance interaction were commonly investigated using plot or stand level models (e.g., Simard et al., 2011). Many such studies covered a large geographic extent (e.g., Jain et al., 2020), but could not account for contiguous spatial processes fundamental to landscape ecology. Similarly, aspatial forest optimization models aggregate local dynamics across broad spatial units without spatial interaction within or between units. This is also the case for dynamic global vegetation models (DGVMs) that aggregate disturbance processes within course-resolution cells without any spatial interactions within or between cells. We also identified a couple fine-resolution land surface model applications that lacked spatial processes. All such aspatial model applications were excluded from our review. It is also possible that the initial use of the term "landscape" may have neglected models from some disciplines (e.g., hydrology) that use alternative terms (i.e., watershed, basin, etc.) to define their geographic extent.
Human land management was included in most of the disturbance interaction studies surveyed (79%), though only 7 studies included land use/land cover (LULC) change as a disturbance agent. By comparison, the most frequently investigated natural disturbance type from our review was fire (65%), followed by wind (38%) and insects (37%). Less frequently investigated disturbance types included water balance (drought or flooding; 15%), mammalian browsing/grazing (10%), and mass movement (i.e., erosion, debris flows, landslides, avalanche; 7%). Only two interaction studies included a simulated disease, perhaps due to a lack of empirical studies (Cobb and Metz, 2017).
Consistent with past reviews of disturbance interactions more generally (see Section "Empirical Approaches"), studies identified by our systematic review were unevenly distributed across continents and biomes ( Table 2). Nearly two thirds of the studies were from North America, primarily from the United States and Canada, including nearly all of the boreal examples. Close to a quarter of the studies were from Europe, focused on temperate and alpine systems, but also including some Mediterranean examples. Studies from Asia accounted for about 10 percent of the total, distributed a little more evenly across biomes, including several studies from China. We found only 3 studies from Australia, one from South America, and none from Africa. Clearly there are underrepresented regions of the world within this field of research, most notably within tropical and subtropical systems. However, some of this uneven distribution may be attributed to our search methodology. We suspect the inclusion of the keyword "savanna * " would have expanded the geographic scope a bit further. Searching on formal model names, while clearly essential in terms of the numbers of studies identified, may have biased the results to regions of the world where those models are applied. Custom models implemented in more generic software such as R were identified within the first search, but not the second name-specific search; this method may have also missed studies from under-represented regions. Nonetheless, our systematic review represents a strong cross-section of the landscape modeling investigations of disturbance interactions.

Forest Landscape Models
A specific class of models known as forest landscape models (FLMs) dominated the literature of landscape  Numbers under each disturbance type reflects the different levels of process integration applied for each type ( Table 4). Note that the absence of a value under a given disturbance type does not mean the software is not capable of simulating that type, it simply was not present in the pool of studies from our review. Bold underscore indicates dominant mode in applications. Numbers in parentheses indicate an external coupled model.
disturbance interactions, representing 92% of studies from our systematic review (Table 3). FLMs were first developed in the late 20th century (Baker, 1989;Mladenoff and Baker, 1999;Gustafson, 2013) and have greatly diversified since then (He and Mladenoff, 1999;Schumacher et al., 2004;Lischke et al., 2006;Scheller and Mladenoff, 2007;Keane et al., 2011;Seidl et al., 2012;Wang et al., 2014, among others). Landscape modeling software (Fall and Fall, 2001; e.g., Bolte et al., 2007) further enabled the rapid development of customized FLMs. Despite this diversification, a common FLM structure has emerged integrating vegetation dynamics (succession, growth, mortality, senescence, etc.) with spatiallyexplicit disturbances affecting those dynamics. In some fashion, both vegetation and disturbance dynamics may be influenced from the "bottom up" by biophysical drivers such as soil characteristics and terrain, and from the "top down" by climatic drivers such as precipitation and temperature patterns, and in some cases atmospheric conditions (e.g., CO 2 concentrations, flux of photosynthetically active radiation). Model development during the last two decades improved the architecture of the models, in terms of modular design and scaling methodologies, the relative balance between phenomenological patterns and simulated processes, and the diversity of disturbance agents that can affect system dynamics. Active development and widespread application of FLMs precipitated several reviews in the last decade (Seidl et al., 2011;Keane et al., 2013Keane et al., , 2015Shifley et al., 2017).
Disturbance processes simulated within FLMs and other landscape models share common elements of initiation, spread or contagion, intensity, and effects, where opportunities for interaction can occur within each of these processes (Figure 3 and Table 3). Direct interactions between disturbances are often approximated through some surrogate for underlying process. For example, time since fire may serve as a proxy the process of fuel build-up. Indeed, Mladenoff et al. (1996) extended this approximation to enable one of the first formal disturbance interactions within an FLM, where a recent wind event increases the intensity of a subsequent fire event. Such logic becomes increasingly complicated, however, as one simulates multiple overlapping disturbances or other disturbance types, such as harvesting or insects (Figure 7). To what degree does the surrogate (in this case, either time since fire or time since disturbance) approximate the process (fuel accumulation) when multiple disturbances are at play? There is the potential for compounding error as additional direct disturbance interactions are added (Figure 7). This issue relates to any of the interaction types that can occur within a system (Figure 3). The complexity of state and transition models, for example, multiply with the number of disturbance processes that are included, because transition probabilities must be specified for every combination of possible disturbance sequences (James et al., 2011a;Daniel et al., 2017).
The solution to the increasing intractability of "direct" interactions via surrogates is to model the underlying process, such as fuel dynamics, as a state variable. Here, the interaction becomes indirectly mediated by the internal processes underlying the dynamics of the state variable, that in some cases, may be more elegantly modeled than attempting to keep track of overlapping surrogate variables (Figure 7). Such methods have other benefits, such as allowing for external drivers such as climate to further mediate the interaction (Figure 5), and to model the lagged effects appropriate for a given process and location (Burton et al., 2020). Simulation of multiple interacting disturbances is likely related to the decline of the state and transition model approach within FLMs in favor of tracking individual tree species that can respond individualistically to a wider range of disturbance types and their various combinations. Only 10 percent of the studies in our review applied the state and transition model approach (Table 3). A more recent innovation in FLMs is the simulation of tree stress via carbon balance (e.g., Keane et al., 2011) or even more explicitly in the form of nonstructural carbohydrate reserves (Seidl et al., 2012;DeBruijn et al., 2014). Such innovations enable trees to respond to the accumulation of multiple stressors that include both disturbances and environmental drivers.
A related model design choice relates to trade-offs in "static" versus "dynamic" disturbance regimes. A static disturbance regime is one parameterized according to its statistical properties (e.g., rotation, event size distribution, intensity distribution, etc.; Figure 2). The focus of such applications is typically the response of the system to a given disturbance regime, without feedback to the regime itself. A fully dynamic disturbance regime is one where the statistical properties of the regime are an emergent property of the system dynamics . In practice, any of the individual components of a given disturbance type (i.e., Figure 3) can be simulated as either static or dynamic. Most contemporary FLMs simulate disturbance effects dynamically, where mortality is contingent on the relative susceptibility of forest types, tree species, age classes, or some combination of these variables. Many introduce stochastic elements that may be related to process (e.g., likelihood of insect damage or fire ignitions related to vegetation characteristics). Some models offer different options where the user can specify the relative balance of static vs. dynamic components (e.g., Sturtevant et al., 2009;Keane et al., 2011;DeJager et al., 2017).
We developed a classification scheme to address the degree of process integration (i.e., balance of static vs. dynamic components) within a specific disturbance agent within a given landscape disturbance application (Table 4), and then applied that scheme to the different disturbance agents simulated across the studies identified by our systematic review ( Table 3). While the classification scheme was both logical and readily applied given proper documentation -we found this supporting documentation was sometimes lacking within individual studies. Many referenced earlier model applications, or original model publications, requiring more extensive review to fully understand the nature of the disturbance interactions simulated. Other applications required assumptions on our part where the documentation was insufficient, including vague descriptions of the specific model options applied, or references to model code requiring literacy of the underlying programming language. Nonetheless, we were able to document ranges of process integration across disturbance types and model software platforms (Table 3), including dominant modes of process integration, and instances where more detailed process integration was possible via the coupling of models (see Section "Coupled Models"). We found that while the trend in disturbance modeling is toward more explicit modeling of disturbance processes and their emergent effects on system dynamics and landscape structure (Gustafson, 2013;Keane et al., 2015), static disturbance regimes are still commonly applied. The relative balance of static vs. dynamic components in application was also unevenly distributed across different disturbance types and model approaches -including non-FLM applications ( Table 3).
Some contemporary FLMs are better described as model frameworks, where different trade-offs in process vs. surrogates, static vs. dynamic disturbance regimes, and the number and types of disturbance mechanisms simulated may be customized according to both the question addressed and the available supporting data. The most frequently applied model software from our review was LANDIS-II, an FLM model framework accounting for 59% of all studies ( Table 3). A case study pair of applications of this model in north central Minnesota (United States) illustrates the framework approach (Lucash et al., 2017(Lucash et al., , 2018 (Figure 8). The authors applied a highly mechanistic succession option (Scheller et al., 2011a) where tree-species cohorts respond individualistically to different disturbance types and compete for above and belowground resources (light, water, nitrogen) via a set of life history traits and parameters. The first study focused on how the combination of climate change and forest management affected forest resilience following wind events (i.e., forest recovery rates), applying a factorial experiment that crossed climate scenarios with forest management scenarios (Lucash et al., 2017). The second study focused on the degree to which multiple disturbance regimes affected tree mortality rates (i.e., effects). Results of the latter supported antagonistic interactions, where total tree mortality caused by multiple disturbance regimes applied simultaneously was less than the sum of the mortality from each of the respective disturbance types applied independently. Closer examination of the simulation methods reveals a combination of process integration methods were applied to different disturbance types (Lucash et al., 2018). Wind disturbance was simulated as "mostly static" regimes with characteristic size, shape, and rotations of three classes of wind patterns (microburst, tornado, and derecho), each with species and age-specific impact patterns. Such static regimes are consistent with the lack of feedback between vegetation and wind event frequency or intensity, but cannot address spatial factors such as the effect of forest edges on tree susceptibility to windthrow (see Section "Disturbance Interaction Process and Terminology"). Fire regimes were simulated as mostly dynamic functions of the vegetation, fuel dynamics, and climaterelated fire weather. Limitations of this approach in terms of emergent fire behavior are shared across virtually all FLMs: rules governing fire spread do not account for critical feedbacks between fire behavior and weather that fundamentally change those rules (such as long-distance dispersion of fire embers) (see Section "Cross-Scale Interactions"). Recent wind or insect disturbance could modify fuel attributes as direct interactions (Figure 7). The three different insect disturbance agents included intermediate blends of dynamic and static regime elements, where the extent of each disturbance was directly related to availability of vulnerable host, but the frequency and intensity of outbreaks was predefined. The exception was the patterns of mortality in oak (Quercus spp.) species by the two-lined chestnut borer (Agrilus bilineatus Weber) that occurred in instances where forest tent caterpillar (Malacosoma disstria Hübner) defoliation overlapped in time with drought events, reflecting a surrogate for tree stress defined by direct interactions with defoliation and climate, respectively (Figure 8).

Questions and Approaches Across Disturbance Types
We summarized the class of explicit disturbance interaction questions addressed across the different disturbance types (Table 5), using the conceptual advances summarized in Figures 3, 5. Frequencies of studies across questions and disturbance types were limited to those who either quantified the effect via their experimental design, or explicitly defined the nature of the simulated interaction within their methods. "Ecosystem response" was the one exception (see below). Linked disturbance interactions were specified according to which disturbance types were affected by another type. While looped disturbance types of the same disturbance (e.g., the effect of fire frequency on future fire likelihood) were not included in Table 5, instances where disturbance of a given type included multiple agents of that type (e.g., more than one insect agent) were included in the tabulation. Linked disturbances generally related to questions of disturbance resistance (Buma, 2015), but some studies aggregated responses across multiple disturbance types. In these cases, we tallied the study under "resistance (general)" ( Table 5). Articles addressing resilience interaction questions were defined as studies that quantified the rates of return to a reference condition, shifts to alternative vegetation states, direct comparisons with HRV, or some combination of these responses. Climate-mediated disturbance questions refer to those studies where a given disturbance type was either modified to be consistent with future climate, or dynamically linked to climate. In the case of land management, climate-mediated disturbance involved some form of intentional climate adaption, such as assisted migration via planting. Cross-scaled disturbance frequencies refer to those studies that quantified cross-scaled interactions via their experimental design. Ecosystem response frequencies refer to those studies quantifying state variables other than those expressly linked to the disturbance interaction questions above (e.g., compositional trends). We note that most studies, whether self-identified as disturbance interaction studies or not, have an interest in capturing system dynamics under multiple interacting disturbance regimes. Frequencies of ecosystem responses in Table 5 therefore reflect the number of studies that included ecosystem responses regardless of whether the specific effect of a given disturbance type could be quantified in those response variables via the experimental design.
The most commonly investigated linked disturbance was fire (Table 5), with most investigations focused on how land management actions such as fuel treatments or logging affected fire severity or extent. Far fewer studies investigated fire as linked to insect disturbance. Insect disturbance was investigated as linked to a broad range of disturbance types in half the studies for which it was included. Wind was less often investigated as linked to other disturbances, with the largest exceptions from studies in central Europe. Studies including the major disturbance types (management, fire, wind, and insects) investigated questions related to system resilience in roughly a quarter of their respective 4 | An index reflecting the relative degree of process integration applied for a given disturbance type and more complete description, with index values corresponding to the numbers applied within disturbance type columns in Table 3.

Degree of Integration
Description (see text for details) (1) Susceptibility index Surrogate variable such as host abundance or fuel distributions (2) Fully static regime All components predetermined, may be applied stochastically (3) Mostly static regime Most components predetermined/stochastic, but with dynamic effects (4) Intermediate static-dynamic Intermediate mix across incidence, extent (spread), and intensity components. Can include simple climate constraints.
(5) Mostly dynamic regime Dynamic across all disturbance components, but without explicit mechanisms for cross-scaled interactions. Climate or weather often a driver of behavior.
(6) Fully dynamic regime Reserved for spatialized physics models and/or individual-based models that enable emergent behavior based on simulated state variables and human and environmental drivers. Explicit cross-scaled interactions may be simulated.
FIGURE 8 | Example modeling study investigating multiple disturbance interactions, where disturbance interactions are primarily indirect via site-level state variables such as composition, age structure, and fuel loading. Climate drivers directly affect just small subsets of disturbances directly, but indirectly affect state variables that can lead to emergent system behavior as a function of two scenario factors (i.e., climate and management). (Adapted from Lucash et al., 2018).
cases. For less-studied disturbance types, questions addressed varied widely, though the sample size was low. Drought was primarily investigated as a climate-mediated disturbance (16 of 22 studies) and related questions of forest resilience (6 of 22 studies). Six of 15 studies including mammalian browse investigated questions of ecosystem resilience, often related to the loss of forest cover. Almost all of the studies including mass movement disturbances focused on linkages to other disturbance types that either removed forest cover or otherwise destabilized soil in steep terrain.
Half of the studies from our review included some form of climate mediation of disturbance regimes ( Table 5). Of those, 61 studies evaluated the degree to which forest management scenarios could mitigate the less desirable effects of climate change, such as elevated mortality (i.e., resistance) or transition to alternative system states (i.e., resilience). For example, Cadieux et al. (2020) investigated how projected climate change interacted with logging and fire disturbance to affect bird species habitat in the boreal forests of Alberta, Canada. In this case, the focus was on the vegetation response rather than on disturbance processes. They therefore applied mostly static fire and drought regimes parameterized differently for each climate scenario, based on empirical studies (Boulanger et al., 2014;Chen et al., 2016, respectively). The combined effects of climate change in combination with the associated elevated fire and drought mortality resulted in transitions from conifer-dominated Linked disturbance interactions specify which disturbance types (columns) were affected by another type (rows). Disturbance questions relate to the framework defined by Figures 3, 5. See text (Section "Questions and Approaches Across Disturbance Types"). ** Frequency of studies in bold italics text refer to the subtotal for the question type (rows) or disturbance type (columns). Subtotals are generally less than the sum of the respective column or row due to inclusion of multiple disturbance types and questions within individual studies.
forests to either deciduous forest or open habitats, with important implications for bird habitat. How this change to alternative forest states might interact with the fire regime itself could not be quantified using this approach. Keane et al. (2019) explicitly focused on emergent fire dynamics related to climate change, and how such behavior might be mitigated via either fire suppression or fuel treatment rates. In this study, they used simulations under historic climate and without human intervention to define the multivariate HRV for three contrasting landscapes in the northern Rocky Mountains (United States). This approach, combined with their process-based modeling framework (Keane et al., 2011), enabled them to draw inferences with regards to tipping points in system behavior that indicated loss of system resilience relative to the benchmark HRV, and the degree to which human action might affect system resilience in the face of climate change.
The relative balance of static to dynamic process in disturbance agents define the degree to which reciprocal interactions among those agents may be investigated, and therefore must be specified carefully with regard to the questions asked and addressed. For example, Chen et al. (2011) investigated the effects of larch caterpillar (Dendrolimus superans Butler) defoliation of different intensities on fire regime. However, key assumptions governing the fire regime, such as whether fire size and rotation was an input or an output of the model, were not provided, making it impossible to evaluate their results regarding insect effects on response variables such as fire frequency. Sturtevant et al. (2012) investigated reciprocal interactions between spruce budworm and fire disturbance and demonstrated that budworm defoliation disturbance reduced live ladder fuels, and therefore decreased fire extent over long time scales (centuries), despite increasing fire spread and intensity at shorter (decadal) scales. Specifically, the authors examined how an intermediate static-dynamic insect disturbance regime affected a mostly dynamic fire regime. While they could directly address the question of budworm effects on fire regimes, they could not fully address the reciprocal question (i.e., how does fire influence the insect disturbance regime?). Temperli et al. (2015) coupled a spatially-explicit bark beetle model with the LandClim FLM (Schumacher et al., 2004) to investigate interactions among spruce bark beetles (Dendroctonus rufipennis Kirby), fire, wind, and climate. In this case, wind disturbance was simulated as a mostly static regime, fire was intermediate (i.e., ignitions and size distributions were static inputs, but spread and effects were dynamic), drought stress was linked to climate inputs, and beetles responded dynamically to vegetation, drought, host abundance, and windthrown host resources. Analogous to Sturtevant et al. (2012), the authors could address the effects of wind and fire on bark beetle disturbance behavior, but not the reverse. They found that climate change-driven drought events would promote beetle outbreaks, but the combination of the two mortality agents progressively reduced long-term susceptibility to bark beetle outbreaks. As a final example, Seidl and Rammer (2017) simulated emergent dynamics of bark beetles responding to wind disturbance and climate change in the Austrian Alps, concluding that climate change amplified the outbreaks initiated by wind events, in large part by shifting the optimal development rates of the beetles to higher elevations where its spruce host was more concentrated. Since the authors applied disturbance agents that each responded dynamically to both climate drivers and ecosystem state variables, they could make inferences about reciprocal interactions among disturbance agents, as well as the mediating role of climate.
Some of the constraints on the types of disturbance interaction questions that may be addressed via landscape modeling are imposed by the design constraints of the selected model software (e.g., Table 3). Nonetheless, model designs are ultimately constrained by supporting knowledge and data. In the insect examples above, there is greater consensus on the landscape factors controlling bark beetle dynamics than there are for defoliator dynamics (Kneeshaw et al., 2021). Indeed, bark beetles were included in all four of the disturbance interaction studies that explicitly addressed cross-scaled interactions ( Table 5). Of those, one approximated non-linear responses of bark beetles to host abundance via percolation theory (Seidl et al., 2016a). The remaining three explicitly modeled emergence and spread of bark beetle agents as they interacted with environmental drivers such as temperature, food resources including both pulses provided by windthrow and host concentrations impacted in part by drought stress (Temperli et al., 2015;Dobor et al., 2020;Honkaniemi et al., 2020) (see Figure 6). While a few model designs may have allowed for cross-scale interactions in other studies (Table 3), they were not explicitly quantified so their importance cannot be evaluated.

Modeling Anthropogenic Forest Disturbance
Forest management was frequently investigated as a linked factor affecting other disturbance types, however, the reverse (effects of natural disturbance on land management) was far less common ( Table 5). Eleven of the studies linking management to other disturbances included explicit salvage logging and/or sanitary logging, representing post-disturbance harvest of felled trees and preventative treatments to reduce insect impacts, respectively. Given the quantity of empirical studies investigating salvage logging, it is somewhat surprising this process is less investigated via landscape modeling. FLMs most often model forest management practices as partially dynamic functions of landscape conditions (Table 3). Harvest locations are often selected via ranking algorithms such as oldest harvested first or biased toward species or forest types of interest (Gustafson et al., 2000). Harvest intensity and effects can likewise be defined to reflect a wide range of management prescriptions. However, the rates (i.e., extent) of any given treatment is generally predefined within scenarios. By contrast, forest planning models represent a different class of models that approach the question of how much to harvest dynamically based on the temporal changes and sometimes spatial distribution in the timber resource (Pukkala, 2013). Disturbance interaction questions from a forest optimization point of view often revolve around accounting for natural disturbances that can add uncertainty to the optimization process (Savage et al., 2011). There is a well-developed literature on spatial forest planning models that focus on the drivers of timber management such as growth and yield, access, road-building, and especially costbenefit trade-offs that FLMs generally lack (Baskent and Keles, 2005;De Pellegrin Llorente et al., 2017). Nonetheless, such approaches have struggled with the stochastic and contiguous nature of natural disturbances, and therefore either simplistically account for such processes as aspatial rates of disturbance or ignore them altogether (De Pellegrin Llorente et al., 2017). Occasionally these divergent approaches have been blended to account for both human and ecological drivers of disturbance interactions. For example, forest optimization models have been used to define harvest rates for scenarios implemented in an FLM (e.g., Daniel et al., 2017). Raulier et al. (2014) used the SELES modeling platform to apply an iterative harvest optimization process -one that included spatially contagious processes such as road building and salvage logging where existing roads enabled rapid access. They used the model to simulate reciprocal interactions between fire and timber harvest to better quantify risk to timber supply under climate change in boreal Canada.
More recent methods for integration of socioeconomic drivers that contribute to both reciprocal interactions with forest landscape systems include demand and allocation models and also agent-based models (An, 2012;Egli et al., 2019). Developers recently incorporated agent-based logging and salvage logging within the iLAND model (Rammer and Seidl, 2015). Such development enabled simulation experiments examining reciprocal interactions among wind, European bark beetle, and forestry practices in the face of climate warming (e.g., Dobor et al., 2020) as illustrated in Figure 6, and account for the studies linking land management to wind and insect disturbance in Table 5. Agent-based forest management approaches are now being implemented in other FLM platforms (Spies et al., 2017;Sotnik et al., 2021). Land-use intensification and other anthropogenic disturbance activities have profound effects on landscape legacies and disturbance regime shifts (Watson et al., 2014). An emerging hotspot of research activity therefore lies in the investigation of reciprocal interactions between human activities and ecosystem processes-falling broadly under the scope of coupled human and natural systems (Wimberly et al., 2015). Anthropogenic and natural disturbance interactions might be modeled using integrated modeling platforms (Verburg and Overmars, 2009), where the combined effects of disturbance events and climate change can be explored via experimental designs that include either contrasting or a range of climatic and land-use intensification scenarios (Becknell et al., 2015). While six studies from our review used outputs from human land use change models to drive land cover change within FLMs (e.g., Duveneck and Thompson, 2019), none coupled land used change dynamically with forest harvest or natural disturbances to investigate their reciprocal interactions. One study applied an FLM to estimate forest productivity under alternative climate change scenarios and coupled those outputs with that of yield models from other economic sectors (agriculture and livestock) to drive an economic land allocation model to simulate land conversions under the twin drivers of climate change and human incentives (Briner et al., 2012). This example represents the closest approximation to the complete system portrayed in Figure 1 of any included in our review.

Coupled Models
As illustrated by the coupling of FLMs with land cover change models in the previous section, some authors of articles in our review extended the capabilities of a given model by coupling them with one or more additional models. Different models may be linked via a "loose coupling, " where the output of one model is used as input for another (i.e., the most common form, see Section "Modeling Anthropogenic Forest Disturbance"), or "tight coupling", where the dynamics of one model interact with the dynamics of another (Sturtevant et al., 2007). Loose couplings of models are widely used as a scaling technique, for example by deriving the input parameters of a coarse-scaled model from a finer-scaled process model (Urban et al., 1999;Xu et al., 2012;Boulanger et al., 2018;Huang et al., 2018) or down-scaling vegetation or disturbance behavior from coarsescaled models such as DGVMs (Halofsky et al., 2014;Turner et al., 2015) (Table 3). FLM outputs are also commonly used as inputs for wildlife habitat or viability models that may be considered a dimension of ecosystem response to disturbance interactions (e.g., Longru et al., 2010;Scheller et al., 2011b;Nitschke et al., 2020).
Apart from human factors such land cover change or forest optimization, use of coupled models to address disturbance interactions centered on better representation of hydrologic and mass movement processes currently lacking within FLMs (Table 3). Research by Saksa et al. (2020a,b) coupled an ecohydrology model with a stand growth and yield model (for vegetation dynamics) and a fire behavior model to evaluate the integrated effects of fuel treatment patterns on fire risk and hydrologic outputs from small watersheds of the Sierra Nevada Mountains in the western United States. A planar surface model was used to define the flood disturbance regime for an FLM application to a Mississippi floodplain system in the central United States (DeJager et al., 2019). The combination of forest cover, disturbances, and impervious surface outputs from a coupled land cover change and FLM model system were used to project flood risk via a simple hydrologic model (Thompson et al., 2016). Scheidl et al. (2020) used the projections of an FLM to drive a physically-based spatial hydrologic model coupled with a slope-stability model to evaluate the integrated effects of climate change, forest management, and natural disturbance processes on shallow landslide frequency within the eastern Alps of Austria.

Hybrid Empirical-Simulation Modeling
Hybrid empirical-modeling studies leverage the synthetic knowledge within process models to relate data representing measurable system components to processes underlying disturbance interactions that are either difficult or impossible to measure directly. For example, McGuire and Youberg (2019) empirically investigated the effects of repeat burns on both soil water infiltration rates and the occurrence of debris flow disturbances in dry montane forests of southeastern Arizona, United States. They then coupled remote sensing and terrain inputs with a high-resolution hydrologic model to demonstrate that differences in soil infiltration between once and twice burned soils were sufficient enough to affect threshold rain intensities necessary to generate debris flows. Similarly, Hancock et al. (2012) measured windthrown trees and consequent erosion patterns following a cyclone impacting northern Australia. They then used this data within the geomorphic landform model SIBERIA (Willgoose et al., 1991) to investigate ecohydrological controls (including windthrow patterns) on long-term erosion and catchment evolution. Økland et al. (2016) investigated a critical transition in bark beetle outbreak dynamics where beetle populations move from trees injured or killed by wind events to attack healthy trees (i.e., windfall-driven to patch-driven dynamics). They compared alternative models representing different hypothesized drivers within an individual-based model framework and found that the model driven by beetle aggregation was the most consistent with their empirical data on the spatial progression of an outbreak in Tatra National Park in Slovakia. Analogous hybrid modeling by Seidl and Rammer (2017) now serves as the foundation for the individual-based model of beetle disturbance contained within the iLAND FLM.

Model Transparency, Uncertainty, and Confidence
A primary challenge to conducting a systematic and consistent review of landscape models featuring disturbance interactions was the lack of standards by which to classify the implementation of a given disturbance agent in terms of its relative degree of process integration -or for that matter, documentation standards for FLMs more generally. There is an increasing call for standardization of methods in analogous classes of models, including agent-based model designs (Grimm et al., 2010) and species distribution models (Zurell et al., 2020). Table 6 represents a partial example of how documentation standards relevant to disturbance interactions can convey the nature of the interactions being simulated and investigated. Minimum standards should include the disturbance interaction type simulated (Figure 3), the response variable quantifying its effect, and the dynamic components of a given disturbance agent that dictate which interaction types are justifiably investigated. This documentation would facilitate a rapid classification -like that in Table 4 -so that the nature of the simulated interaction is transparent to the reader. Other essential documentation should include a clearly articulated study design and basic information such as the spatial and temporal grain and extent of simulated processes ( Table 6). The latter is coming increasingly important to clarify as the temporal and spatial resolution (grain) has become more specific to individual processes as model designs become more refined ( Table 6). Additional relevant documentation includes the level at which disturbance effects are applied (patch/cell, age class, species, cohort, tree, etc.), the spatial implementation of the disturbance agents (spread, dispersion, etc.), and how the various disturbance agents are interspersed in time. All such documentation should be specific to the model application at hand, rather than the potential capabilities of a given model software, as this was sometimes conflated in the articles we reviewed. Outputs from spatial models are complex enough without confusing artifacts of modeling assumptions with emergent behavior based on underlying process.
We do not assert that complex models are better than simple models. Rather, that one should account for the processes underlying the question being investigated to the extent that data and knowledge can support them. The contemporary argument for design of models based on first principles is that such models should make more robust predictions under the increasingly novel conditions of changing climate, atmospheric composition, and disturbance agents (Cuddington et al., 2013;Gustafson, 2013). We note also that many landscape modeling studies did not explicitly investigate disturbance interactions per se, but rather included multiple disturbance agents to simply account for all the relevant structuring processes within their forest system. Regardless, it is important to acknowledge and account for the cumulative effects of model error propagation as the number of interacting components increase and from one scale to another (Rastetter et al., 1992). To be more explicit about the relative importance of multiple processes and errors, Cressie et al. (2009) promoted the use of hierarchical statistical modeling that allows the quantification of the relative importance of uncertainty stemming from the data model (i.e., data error), the process model (i.e., model specification error), and the parameter model (i.e., parameter estimation error). There are, however, multiple additional sources of uncertainties-ranging from the data, models, processes, knowledge of the processes, and climate scenarios-that will affect the precision and accuracy of model projections.
To quantify the relative importance of these various types of uncertainty, one could compare the outcomes of several models (Cheaib et al., 2012). A recent comparison of four "state of the art" FLMs demonstrated wide divergence in model outcomes, in large part due to design choices affecting vegetation response to disturbance (Petter et al., 2020). This result suggests a certain degree of humility when reporting the results of landscape models -particularly where interactions by multiple disturbance agents may compound error associated with cumulative vegetative responses. Scheller (2018) suggested that parameterization and "model confidence" are twin challenges facing contemporary landscape modeling studies. The increased burden of input parameter specification that generally comes with replacement of simple statistical relationships with more mechanistic detail may potentially undermine the perceived robustness of process-based models if those parameters are either poorly understood or quantified (McKenzie et al., 2019). Relevant knowledge acquisition necessary for parameter estimation will require multidisciplinary research crossing different domains in scale, integrating experimental research, monitoring of earth systems, advanced statistical analyses, and next generation models assembled from these insights to determine process-based responses to future climatic conditions. Knowledge gaps may further require more novel hybrid modeling frameworks (e.g., Talluto et al., 2016) that integrate process-based and statistical models to include both bottom-up processes and drivers and top-down drivers and account explicitly for feedbacks across scales (Figure 1). Such integrative models require knowledge about ecological responses based on targeted experiments in the field and carefully designed data analysis of existent data (Poisot et al., 2016).
Model confidence refers to the degree to which independent readers trust the utility of the model results -in this regard independent validation data remains a fundamental issue underlying confidence in FLMs (Shifley et al., 2017;Scheller, 2018). Validation issues may be particularly acute in the area of disturbance interactions given system sensitivity to the timing and order of disturbance mechanisms and the individualistic responses of the vegetation. While there were examples within our review that evaluated their results using empirical observations (Henne et al., 2013;Karam et al., 2013;Schwörer et al., 2014;Thom et al., 2018), these examples were rare. Validation of model components in combination with targeted uncertainty and sensitivity analyses may also improve model confidence (Cressie et al., 2009;Xu et al., 2012;Scheller, 2018;McKenzie et al., 2019). More specifically, model uncertainty refers to the degree to which model inputs are known as well as their associated measurement error, whereas sensitivity analyses focuses on the degree to which model inputs have influence on the model outputs (McKenzie et al., 2019). Guidance on effective uncertainty and sensitivity analyses can be found across disciplines for simulation modeling more generally (e.g., Saltelli et al., 2006;Wagener and Pianosi, 2019).

Synthesis: New Insights
We found a broad spectrum of questions related to disturbance interactions in forested landscapes have been addressed by an equally diverse range of modeling approaches (Tables 3, 6). Notably, it is the more intractable of these questions that are also the most challenging to quantify empirically, and these require the support of simulation models to investigate. Integration of recent conceptual advances (i.e., linked disturbances, compound disturbances, cross-scaled interactions, and cascading effects) often require more detailed processes to simulate. Among the most important are non-deterministic tree-species competitive interactions and disturbance responses via life history traits, explicit modeling of stand structural characteristics, physiological traits of tree species that enable dynamic responses to climate and biophysical variables, and various integrations of these factors (Keane et al., 2011;Scheller et al., 2011a;Seidl et al., 2012;DeBruijn et al., 2014;Wang et al., 2014). Depending on the questions asked, modeling disturbances require ecological detail necessary to effectively capture their incidence, extent, intensity, and/or effects (see examples within Section "Questions and Approaches Across Disturbance Types"). We are approaching a point where some of the most process-based disturbance models, such as physical fire models, individual-based insect models, and agent-based human disturbances are possible either within FLMs or coupled model architectures (Rammer and Seidl, 2015;Bentley and Penman, 2017;Seidl and Rammer, 2017). Nonetheless, tradeoffs between model elegance and complexity remain (McKenzie and Perera, 2015). For example, the parameter demands of highly complex models may exceed the data available to parameterize them (Scheller, 2018) (see also Section "Model Transparency, Uncertainty, and Confidence"). We also observe that the most mechanistic examples reviewed here were comparatively limited in the extent at which they were applied (e.g., thousands versus millions of ha, years versus centuries), illustrating practical tradeoffs between ecological detail and scope remain despite advances in computing power ( Table 6).
The study of disturbance interactions through landscape modeling reflects the asymmetry in systems and disturbance agents of the empirical studies (see Section "Empirical Approaches"). Studies were concentrated primarily and temperate, boreal, and alpine regions of North America and Europe ( Table 2) and were almost absent from the southern hemisphere. Perhaps it is no coincidence that the prevalent disturbance regimes of subtropical and tropical regions (e.g., mass movement, hurricane/cyclones, mammalian grazers; Burton et al., 2020) were less commonly investigated in the context of disturbance interactions ( Table 5). We found that human disturbance -by far the most commonly simulated disturbance type -was primarily simulated as a driver rather than an interactive component of the simulations. Cross-scale interactions remain an important concern that are particularly challenging to model and quantify. Each of these deficiencies suggests fruitful areas of future research. Finally, we echo the plea made by previous authors to improve model transparencyperhaps via documentation standards such as that developed for other classes of models.

FUTURE DIRECTIONS
Our review demonstrates substantial progress in the study of disturbance interactions necessary to model their effects on forested landscapes. We now have far greater clarity in the conceptual underpinnings of disturbance interactions, a wealth of empirical studies to draw from, and more advanced model frameworks capable of simulating multiple disturbance types across a range of complexities, relative to what was available at the turn of the millennium when such questions were gaining traction. Questions asked by ecologists are becoming increasingly more ambitious, in part due to the expanding data resources and tools available, and in part due to urgently expanding needs for reliable forecasts of natural and managed ecosystems under increasing anthropogenic stressors. Multiple disturbances and drivers as well as their interactions increasingly insert uncertainty and surprise as forested systems adjust to novel conditions (Urban et al., 2016), complicating our ability to project their behavior. It is therefore paramount to unravel the cumulative effects of multiple disturbances and factors (e.g., edaphic, environmental, and climatic) on vegetation change, with the understanding that as the number of interactions increases, so does the uncertainty associated with ecosystem responses. Such complexity underscores the need to include both conceptual advances and empirical approaches in a review that is ultimately about effective modeling of forest landscape disturbance interactions. As so eloquently stated by Seidl (2017), "to model or not to model, that is no longer the question. . .". Conceptual understanding, empirical study, and simulation modeling should continually reinforce one another if we are to unravel the complexities of disturbance interactions in time and space.
The "art" of modeling landscape disturbance interactions calls for balanced approaches and scales: while macroscale, empirical models are used to approximate processes, the integration of micro-and mesoscales mechanisms into process-based models is the way to adapt and plan for the novel conditions of the "Anthropocene" (Bodner et al., 2020). The simulation studies cited in this review demonstrate a broad range of modeling approaches that reflect the diversity of ways investigators balance trade-offs in scale, scope, and ecological complexity when modeling disturbance interactions (e.g., Tables 3, 6). Given the importance of transparency to model confidence (Scheller, 2018), it is critical that modelers be explicit about basic assumptions and methods. We note that many landscape simulation studies in the literature cite the potential capabilities of a given modeling platform while glossing over the specifics of the options they applied. We further caution against loading simulation studies with multiple disturbance types simply because a given model platform allows it. Careful consideration of how the disturbances interact within the context of the model design and the questions being addressed are recommended to avoid compounding model artifacts that can produce misleading results. Furthermore, computing technology may enable modeling of processes that lack the empirical studies necessary to parameterize them-hence the need for a multi-dimensional research agenda. Turner (2010) emphasized the importance of explicitly accounting for human-induced changes and societal landuse intensification demands, as well as societal responses to disturbance regime shifts and rapid land cover changes. The effects of landscape change resulting from disturbance regime shifts can be investigated by comparing alternative scenarios of plausible system states (Titeux et al., 2016;Urban et al., 2016). By contrast, the more human-focused modeling studies tend to focus on a single ecological response variable (Egli et al., 2019). Blending the state-of-the-art in coupled human natural systems modeling and FLMs remains a current frontier, leveraging the advances in each respective discipline (e.g., Briner et al., 2012).
To sum, modeling disturbance interactions across scales epitomizes the broader challenges plaguing ecological research more generally. The dynamic nature of earth systems makes a strong argument for process-based modeling grounded in a mechanistic understanding of the dynamics of individual system components (Cuddington et al., 2013;Gustafson, 2013). Yet the compounding uncertainty of the interactions between system components-including disturbance processes-that can be both non-linear and scale-dependent lie at the heart of the intractability of complex human-ecological systems (Turner, 2010). Clarity in concepts (knowledge), an empirical foundation (data and analysis), and model designs (synthesis and software) will ultimately enhance understanding of complex systems, but the choice of which processes to model explicitly and which processes to aggregate remains a fundamental challenge of our time.