1 Introduction

Spruce budworm (SBW; Choristoneura fumiferana (Clem.)) defoliation is one of the most influential disturbances affecting forest productivity in North America (Cooke et al. 2007; Morin et al. 2007). Periodic SBW outbreaks have caused widespread defoliation across greatly large areas, e.g., over 58 million ha during the 1970s and 1980s in the USA and Canada (Blais 1983; USDA Forest Service 2009). SBW activity has been rising in recent years in Maine, USA, and the neighboring regions, e.g., 9.6 million ha of forests have been defoliated in Quebec, Canada, by 2019 (Kanoti 2017; Ministère des Forêts, de la Faune et des Parcs 2019). SBW outbreaks in the past century usually lasted about 15 years (Blais 1983), during which defoliation occurred across vast landscapes and frequently shifted between severity conditions of harmlessly low to almost complete defoliation over rather large spatial and temporal extents (Peltonen et al. 2002; Chen et al. 2018a).

Given the importance and complexity of SBW defoliation, significant efforts have been placed on modeling the dynamics of SBW outbreaks, predominantly based on evaluations and syntheses of the many processes, e.g., dispersal, predation, and SBW-host relationship (Morris 1963; Royama 1984; Nealis and Régnière 2004), underlining SBW population dynamics (Sturtevant et al. 2015). However, the development and further refinement of these models have proven to be challenging due largely to the lack of comprehensive knowledge of the many processes operating at various scales, as well as difficulties in tracking and quantifying SBW population across the landscape for model development (Nenzén et al. 2017a). Consequently, there currently is no generally accepted theory or generalized modeling framework summarizing SBW population dynamics across the landscape (Pureswaran et al. 2016).

Instead of modeling SBW population dynamics, a few studies empirically evaluated large-scale defoliation caused by SBW (e.g., Alfaro et al. 2001; Hennigar et al. 2008; Nenzén et al. 2017b), which is relatively easier to measure and directly related to forest productivity. These evaluations related SBW defoliation to forest stand characteristics such as age, density, and species composition and formed the basis of defoliation risk mapping commonly used in forestry. As the periodic nature of SBW outbreaks indicates both rises and declines of defoliation over time, recent work has added a temporal component to improve these otherwise spatially and temporally static evaluations of defoliation dynamics. For example, Gray and MacKinnon (2006) and Zhao et al. (2014) summarized 27 and 28 temporal patterns of SBW defoliation, respectively, across eastern Canada. Furthermore, Chen et al. (2018a) developed a parametric model to explicitly predict individual tree annual defoliation through the 10+ years of the duration of an SBW outbreak. Regardless, these prior efforts have primarily focused on stand- rather than landscape-level trends in defoliation.

Despite the consideration of temporal dynamics, what is generally absent in the assessment of the dynamics of SBW outbreaks is the spatial component (Sturtevant et al. 2015). As growing evidence supports patterns of SBW outbreaks across the landscape such as the west to east development, synchrony, and hot spots of defoliation (Blais 1983; Irland et al. 1988; Pureswaran et al. 2016), it is reasonable to assume that defoliation at one location is affected not only by its biotic and abiotic characteristics but also defoliation at neighboring locations. This may be an important cause of the highly variable defoliation of similar stands observed in the same years (Chen et al. 2018a). Nevertheless, explorations of the spatial component in the dynamics of SBW outbreaks do exist. For example, Candau and Fleming (2005) related frequencies of defoliation, in terms of how many years a location was >25% defoliated between 1967 and 1998, to temperatures across Ontario, Canada; Bouchard and Auger (2014) used environmental factors over space and time to evaluate their significance on the initial expansion of an SBW outbreak in Quebec, Canada; and Nenzén et al. (2017a and b) developed a theoretical susceptible-infectious-recovered (SIR) model to showcase that SBW dispersal was necessary to generate cyclic outbreaks.

Despite advancing our understanding of the spatial dynamics of SBW outbreaks, a few limitations in these previous analyses may hinder their general applicability in forest modeling, management, and conservation. First, the variables of interest have primarily been the presence/absence of defoliation or its frequency. Not only are these variables not readily translatable to defoliation commonly measured in the field (e.g., % defoliation) and used in evaluations of forest productivity (Chen et al. 2017a), but how presence/absence is defined (i.e., how much defoliation is considered present) would also greatly affect the results of analyses. Second, there is currently no cross-scale approach as proposed by Senf et al. (2017) that simultaneously evaluates and validates the varied influences of both stand characteristics and environmental factors. This is particularly important considering that environmental gradients like elevation may have affected stand characteristics like species composition. Finally, explicit and data-based parameterization is generally not available, which may impede the validation of the findings, as well as utilizing them in future research and applications (Sturtevant et al. 2015; Nenzén et al. 2017b).

Given this current knowledge gap and availability of long-term data, the goal of this study was to develop a parametric statistical model flexible enough to explicitly evaluate the complex spatial and temporal dynamics of SBW outbreaks. Specifically, this model investigates how defoliation at every specific location and time is related to defoliation at surrounding locations across the landscape and during preceding years of an SBW outbreak. Without key assumptions about the many possible and diverse processes governing SBW population dynamics operating at various scales, our evaluation was based on extensive and detailed defoliation data derived from tree-level observations covering approximately 50,000 km2 and 10 years of the last SBW outbreak during the 1970s and 1980s in Maine, USA.

In addition to large spatial and temporal extents, these data also comprise a wide range of forest conditions and defoliation observations accompanied by detailed tree measurements, from which a variety of stand characteristics were derived and evaluated for their potential effects on defoliation dynamics. These stand characteristics were simultaneously evaluated in combination with environmental factors suspected of affecting SBW defoliation dynamics across the landscape such as wind (Anderson and Sturtevant 2011), temperature (Hardy et al. 1986), and landscape structure (Irland et al. 1988). Outcomes of our model were expected to be directly usable in existing forest growth and yield modeling frameworks and management decision support systems (MacLean et al. 2001; Weiskittel et al. 2011a) like the Forest Vegetation Simulator Acadian Variant (FVS-ACD) with SBW modifiers specifically developed for the region influenced by periodic SBW outbreaks (Chen et al. 2018b). In addition, the outcomes would also directly provide quantitative information supporting forest conservation and management in the face of a pending SBW outbreak on the landscape.

The specific research objectives of this study were to (1) identify influential factors affecting spatial and temporal dynamics of SBW defoliation across the landscape; (2) develop a flexible spatiotemporal parametric model explicitly predicting these dynamics; and (3) apply the developed model to simulate the development of an outbreak simultaneously over space and time under various outbreak scenarios. We hypothesized that SBW defoliation would rapidly and pervasively disperse across the landscape, which would result in stand characteristics being generally non-influential on the ubiquitous event of multiple-year defoliation of forests during an SBW outbreak. However, stand characteristics may potentially significantly affect the severity of defoliation and forests’ growth responses to defoliation as observed in prior analyses (e.g., Chen et al. 2017a, b).

2 Material and methods

2.1 Study area

The subject area of interest was the vast spruce-fir (Picea abies) forests in the state of Maine, USA, which are primarily distributed across the northern parts of the state. Specifically, the extent of our defoliation data is 44.94°–47.30° N and 67.30°–70.73° W, which covers an area of approximately 50,000 km2 and most of the spruce-fir forests in Maine (Fig. 1). Therefore, the developed model was believed to be representative and applicable to all of the spruce-fir forests in the state, where our simulations of the spatial and temporal dynamics of an SBW outbreak would be conducted.

Fig. 1
figure 1

The study area and its location in the USA, where black dots are sample plots

Full size image

The above study area is located in a transition zone between the temperate deciduous forest and the boreal coniferous forest and contains more deciduous trees than neighboring areas (e.g., New Brunswick and Quebec of Canada), which are further north and have also been affected by SBW outbreaks (Irland et al. 1988; McWilliams et al. 2005). Soils in this area are generally low in fertility, acidic, and low in permeability such that poorly drained, loamy soils underlain by loamy basal till are found on lower slopes and uplands, while well drained, loamy soils underlain by bedrock are found on upper slopes (Ferwerda et al. 1997).

Maine has a humid continental climate influenced by the humidifying and moderating effect of the adjacent Atlantic Ocean. Consequently, this area has warm and moist summers along with cold and snowy winters, which are milder than those in the middle of the continent. Periodic SBW outbreaks have occurred for at least hundreds of years in this region, while the last outbreak was in the 1970s and 1980s (Blais 1983; Fraver et al. 2007) and the area is currently on the verge of experiencing another outbreak (MacLean et al. 2019).

2.2 Data

Five independent sets of data were used in this study to derive various measures of defoliation as well as forest stand characteristics, weather conditions, and landscape structure during an SBW outbreak (Chen et al. 2021). While defoliation is the subject of this study, the others are all considered important factors directly related to the developments of defoliation across the landscape. All of these data were geographically referenced and projected to UTM zone 19N with the unit of meter using the datum of NAD83 in our analysis (Chen et al. 2021). Each of these five datasets is described in detail below.

2.2.1 Growth Impact Study

The Growth Impact Study between 1975 and 1985 collected data at 424 0.02-ha sample plots spread across the spruce-fir forests of Maine (Solomon and Brann 1992). Specifically, data from 376 of these plots with known coordinates, which occur throughout the study area and were generally representative of the forests (Fig. 1), were used in this study (accurate locations of the other plots could not be determined). These data contain annual individual tree measurements such as species, diameter at breast height (DBH), and height, from which metrics of stand characteristics such as dominant height (m; height of the tallest tree in a sample plot), relative stand density (additive stand density index/maximum stand density index; Woodall et al. 2005), and host tree percentage (balsam fir (Abies balsamea L.), red spruce (Picea rubens Sarg.), black spruce (Picea mariana (Mill.) B.S.P.), and white spruce (Picea glauca (Moench) Voss) by volume; Chen et al. 2017a) were computed. These metrics were used to represent the development stage, competition, and species composition of forest stands.

Each year between 1975 and 1985, current-year foliage on each host tree within each plot was visually examined for the degree of defoliation and categorized separately into one of five (before 1982; 0, 1–5, 6–20, 21–50, and 51–100%) or eleven classes (since 1982; 0–10, 11–20, 21–30, 31–40, 41–50, 51–60, 61–70, 71–80, 81–90, 91–99, and 100%) of defoliation. Meanwhile, stand-level defoliation (Def) was calculated as the weighted (by volume) mean of individual tree percentage defoliation and averaged between 7 and 18% annually through the study period.

The 11-year defoliation data from 1975 to 1985 covered most of the duration of the last SBW outbreak in Maine as an outbreak typically lasts about 15 years (Blais 1983). Since SBW is endemic to the region, the initiation of an outbreak is difficult to monitor and is considered a major data gap for SBW population ecology (Pureswaran et al. 2016). For example, there was an uptick of observed SBW population around 2016 in Maine but an outbreak has not ensued (Kanoti 2017). While above normal SBW activity was observed since 1972 (Irland et al. 1988), the outbreak was still in the expansion phase in 1975 as average defoliation was 7% across the 376 plots, and 128 of these plots had defoliation <5%. The relatively low defoliation is characteristic of the last SBW outbreak in Maine, which was generally lower than defoliation reported by some previous studies from neighboring regions (Chen et al. 2017a, 2018a).

2.2.2 Daily summaries of temperature

The temperature data are part of the Global Historical Climate Network data available from the US National Centers for Environmental Information, of which only those observed at the 62 stations located in the state of Maine were used in this study. Specifically, daily maximum and minimum temperatures (°C) in July (the timing of SBW moths’ eclosing and exodus flights; Irland et al. 1988) during the years of 1975 to 1984 were used to derive a metric of temperature: the number of days in July that had minimum temperature >14°C and maximum temperature <30°C (the suitable range for SBW moths’ exodus flights; Royama 1984). This measure of temperature was interpolated from the above weather stations to the locations of the Growth Impact Study sample plots using a k-nearest neighbors (k-NN) method, where k was selected to be six based on the lowest root mean squared difference and nearness was measured as Euclidean distances between weather stations and sample plots. This metric was tested as a potential predictor of defoliation in ensuing years (i.e., 1976–1985).

2.2.3 Monthly summaries of wind

An important feature of SBW, which probably has greatly contributed to the dynamic and widespread outbreaks of defoliation, is its strong ability of dispersal through flight, e.g., 200 km in 4 hours (Boulanger et al. 2017). This ability of dispersal is largely facilitated by the prevailing wind (Blais 1983; Sturtevant et al. 2013). The wind data are monthly aggregations of daily wind observations at 70.30°W and 43.64°N in Maine from the same source of the above temperature data. Because wind is known to change frequently at even very fine temporal scales, monthly aggregations were used to reduce noises in data and provide a more representative estimation of the prevailing wind during the period of SBW moth flights. Each year between 1975 and 1984, the fastest mile wind direction and velocity in July were used to develop a metric of wind (W) in the form of W = IW ∙  cos (Wa) ∙ Wv, where the indicator variable IW = 0 if Wa ≥ 90°, otherwise IW = 1, Wa is the angle between the direction of wind (Wd) and the direction of the line connecting the source and destination locations of SBW defoliation (panel a of Fig. 2; all pairwise measures of wind are included in this metric, e.g., between two locations of a and b, this metric is applied from a to b but also from b to a, i.e., all possible source and destination locations of SBW defoliation are considered in this metric and metrics below based on this idea), and Wv is the velocity of wind (m s−1). This metric was applied to all locations in the study area to predict defoliation in ensuing years (i.e., 1976–1985).

Fig. 2
figure 2

Diagrams of the metrics of wind, habitat suitability, and difference in elevation (Δelevation), where Wa is an angle in degree and ΔE is the difference in elevation (m)

Full size image

2.2.4 Land cover

The National Water-Quality Assessment Project refined the US Geological Survey historical land use and land cover data derived from aerial photographs from the 1970s. This refinement showed residential development between the 1970s and 1990s (approximately the same time period of the last SBW outbreak in Maine) and improved the accuracy of land cover classification (Hitt 1994). A 96-meter resolution map produced by that project classifies the landscape in our study area into various types including waterbodies, urban, agricultural, and forest lands.

In particular, forest lands were further divided into three categories, namely deciduous, evergreen, and mixed forest lands. SBW preferably feeds on the buds of a range of coniferous trees primarily of spruce-fir, which were the dominant species in our study area at the time of the last SBW outbreak (McWilliams et al. 2005), and deciduous trees generally were a minor component in the mixed forest lands (Irland et al. 1988; Table 1). We computed a metric of habitat suitability (L) based on the above classification of land cover as the ratio between numbers of host cells (covered by evergreen or mixed forests) and all cells along the line connecting any two locations in the study area (panel b of Fig. 2). In this metric, evergreen and mixed forests are considered corridors while deciduous ones are considered barriers of the dispersal of SBW. This metric is based on observations of near vertical descent of SBW moths from flights (Greenbank et al. 1980) and the proposal by Sturtevant et al. (2013) that SBW moths may be able to sense their hosts below and descend to suitable oviposition sites. Grant et al. (2007) also concluded that terpenes of host trees serve as general oviposition stimuli for SBW moths, and these stimuli are not species specific, i.e., SBW moths sense multiple evergreen host species rather than any one specific.

Table 1 A summary of the forest and environmental attributes in this study
Full size table

2.2.5 Digital elevation model

Variations in SBW defoliation across a range of elevations have been observed in previous studies, and elevation was thought to play a primary role in SBW moth dispersal (e.g., Osawa et al. 1986; Bouchard and Auger 2014). A 30-meter resolution digital elevation model (DEM) covering the state of Maine from the Maine Office of GIS was used to extract the elevation in the study area. This information was used to compute differences in elevation between any two locations (ΔE, positive if the destination location is higher in elevation than the source location of defoliation, and negative otherwise; panel c of Fig. 2) in order to show the vertical structure of the landscape and to supplement the above metric of habitat suitability. The range of elevation in this DEM is between −2 and 1603 m.

2.3 Analysis

Spatial and temporal dynamics of an SBW outbreak across the landscape were summarized based on predictions of percentage defoliation (an indirect measure of SBW populations) at every specific location and time. These predictions were attributed to local factors at the stand scale (i.e., intrinsic stand characteristics, e.g., species composition, to support local SBW population) and landscape factors (e.g., landscape structure and weather conditions affecting the dispersal of SBW). These two groups of factors were additively combined into a spatially and temporally weighted regression model, which influences of landscape factors from neighboring locations on a specific location were weighted by their distances in space and lags in time (weights were calculated using the g1 and g2 functions, respectively, defined below). Obviously, local factors at each specific location itself were not weighted.

This model shares similarities with the one proposed by Meyer et al. (2012) in the decomposition of a transmissible event into endemic (here called local) and epidemic (here called landscape) components, as well as predicting this event in continuous space and time. The latter feature provides the flexibility of utilizing data of various resolutions (e.g., the 30-meter DEM and 96-meter land cover data used in this study) and, more importantly, avoids possible significant effects of arbitrary cell sizes (of raster input and output) on model fitting and prediction.

However, our model is notably different in that percentage defoliation is explicitly used as the response variable. The interest of Meyer et al. (2012) is on the absence/presence of a disease so the response variable is conditional density (probability density) of a disease. SBW defoliation, however, regularly occurs in spruce-fir forests (usually at negligible levels) and affects the growth and survival of trees only when it is at relatively high levels (e.g., Chen et al. 2017a, b). Therefore, often what of interest is the severity (levels) rather than the presence of defoliation in forest ecology and management. Consequently, our model is in the following form:

$$ y\left({s}_i,{t}_j\right)=f\left(X\left({s}_i,{t}_j\right)\right)+{\sum}_{i,j}g\left({g}_1\left(\left\Vert s-{s}_i\right\Vert \right),{g}_2\left(\left|t-{t}_j\right|\right),Z\right) $$
(1)

where y(si, tj) is the percentage defoliation at location i (si) and time j (tj, year 1976–1985), f(X(si, tj)) is the function of the local component at the stand scale, and ∑i, jg(g1(‖s − si‖), g2(|t − tj|), Z) represents the landscape component of the dynamics of an SBW outbreak. Furthermore, the function of the local component is formulated as follows:

$$ f\left(X\left({s}_i,{t}_j\right)\right)=\alpha \bullet \left(D\left({s}_i,{t}_j\right)\bullet H\left({s}_i,{t}_j\right),H\left({s}_i,{t}_j\right),D\left({s}_i,{t}_j\right)\bullet C\left({s}_i,{t}_j\right)\right) $$
(2)

where X(si, tj) are stand characteristics at location i and year j, of which D is the relative stand density (0-1 ratio), H is the dominant height (m), C is the host tree percentage, and α is a vector of three parameters. The function of the landscape component is in the following form:

$$ {\sum}_{i,j}g\left({g}_1\left(\left\Vert s-{s}_i\right\Vert \right),{g}_2\left(\left|t-{t}_j\right|\right),Z\right)={\sum}_{i,j}\beta \bullet {e}^{\frac{-{\left\Vert s-{s}_i\right\Vert}^2}{2\bullet {d}_s^2}}\bullet \left({e}^{\frac{-\left|t-{t}_j\right|}{d_t}}\bullet def,L,W,W\bullet \Delta E\right) $$
(3)

where \( {g}_1\left(\left\Vert s-{s}_i\right\Vert \right)={e}^{\frac{-{\left\Vert s-{s}_i\right\Vert}^2}{2\bullet {d}_s^2}} \) and \( {g}_2\left(\left|t-{t}_j\right|\right)={e}^{\frac{-\left|t-{t}_j\right|}{d_t}} \) are the spatial and temporal kernel functions, respectively, in which ds and dt are range (lag) parameters estimated marginally on a null model using all available observations and having the values of 82 km and 1 year, respectively (Fig. 3), ‖s − si‖ are Euclidean distances between location i and all locations, and |t − tj| are lags between year j and all previous years. Z is the predictors of the landscape component of an outbreak, of which def is the annual stand percentage defoliation in each of the previous years, L is the habitat suitability (0-1 ratio), W is the metric of wind in July of the preceding year (m s-1), ΔE is the difference in elevation (m), and β is a vector of four parameters.

Fig. 3.
figure 3

Log likelihood of spatial and temporal kernel functions over values of range (lag) parameters, of which maximum values are indicated by black dots. While samples remain the same size in different years, their spatial distribution is presented along the x-axis

Full size image

In addition, Def, L, W, and ΔE are (376 ∙ 10) × 376 matrices, where 376 is the number of observations in each year and 10 is the number of years (1975–1984), i.e., these predictors for any single y(si, tj) are vectors of 376 measurements. def is weighted by both the spatial (376 × 376 matrix) and temporal (a vector of 10) kernel functions above, i.e., a matrix of 376 × (376 ∙ 10), while L, W, and ΔE are weighted by the spatial kernel function only. The reason here is that 1 year’s wind seems not directly affected by previous years’ wind (unlike defoliation), and landscape structure (represented by L and ΔE) does not change over a short period of time.

Simulations of spatial and temporal dynamics of an SBW outbreak were performed at a resolution of 10 km across an area of 160 km (longitudinal) × 240 km (latitudinal) shown in Fig. 4 based on the above model and scenarios defined in Table 2. The selection of a 10-km resolution is demonstrative of the methodology and because the interest of this study was in defoliation dynamics at a much larger scale across the complex forested landscape of tens of thousands of square kilometer. By no means was this simulation intended to ignore the complexities of SBW defoliation dynamics across a continuum of scales from branches, trees, and stands to landscapes (MacLean and Lidstone 1982; Gray et al. 2000; Hennigar et al. 2008; Chen et al. 2018a). In fact, this evaluation of defoliation dynamics across landscapes was made available at varying desired scales since the simulations can be performed in continuous space (i.e., at various resolutions) by our model. However, an increase in resolution from 10 to 1 km would result in approximately 104 times of computations with likely limited potential gains in general insight into the complex patterns that emerge at the landscape scale.

Fig. 4
figure 4

Simulated spatial (each grid) and temporal (each row) dynamics of an SBW outbreak in various scenarios (each column) introduced in Table 2

Full size image
Table 2 Scenarios in simulations of spatial and temporal dynamics of SBW outbreaks, where D is the relative stand density, C is the host tree percentage, Wv is the velocity of wind, Wd is the direction of wind, H is the dominant height, L is the metric of habitat suitability, Def is the defoliation, and ΔE is the difference in elevation
Full size table

In different scenarios of simulations, emphasis was rather placed on important factors of host tree percentage and wind, as well as the least influential factor of relative stand density in order to provide generally comprehensive comparisons of the dynamics of an SBW outbreak. Meanwhile, other factors like habitat suitability and elevation evolve at much slower paces. Therefore, any assumed values other than those observed would be less meaningful in reality and thus were not considered in the simulations. The foremost factor of defoliation itself was assumed to be initiated at locations of (68.17°W, 46.32°N), (69.56°W, 46.59°N), (69.57°W, 46.74°N), and (69.57°W, 46.95°N), where the highest-level defoliation of 75% (the same level used to initiate our simulations) was observed in 1975. Defoliation at these locations served as the only source of the simulated outbreaks, i.e., defoliation was set to be 0% initially at all of the other locations in our study area.

All of our analysis was conducted in R v3.2.2 (R Core Team 2015). Parameter estimations of our model were based on the maximum likelihood method. R packages of maps (Becker et al. 2016), maptools (Bivand and Lewin-Koh 2017), raster (Hijmans 2016), rgdal (Bivand et al. 2017), and yaImpute (Crookston and Finley 2007) were used to handle our data. The metric of temperature in July (see above) was not applied to our model because preliminary analysis indicated that it was generally not correlated with defoliation across the landscape and over time.

3 Results

3.1 Model fit and predictive performance

Parameter estimates for the predictors in our model are presented in Table 3 and were nearly all highly significant (p<0.01). The model has an R2 of 0.63 and an overall mean bias of +0.3%. When compared to observed values, 3611 (i.e., 96%) predictions fall within ±1 class of defoliation based on the classification before 1982, while 3099 (i.e., 82%) predictions fall within ±1 class of defoliation based on the stricter classification since 1982. The relationship between predicted and observed values was relatively consistent across the full range of defoliation present in this analysis (Fig. 6 in Appendix). Overall, the model fit statistics and general performance suggest robust behavior despite the high underlying variability of the available data.

Table 3 Parameter estimates, standard errors, and p values for the predictors in our model
Full size table

3.2 Influential factors of spatial and temporal dynamics of an SBW outbreak

Based on these estimates and observations in our data, the most influential factor in SBW defoliation dynamics across the landscape is defoliation itself (at a given forest location and its neighbors from a year before). It has a 31.7 times stronger influence compared to the least influential factor of relative stand density (as a reference, i.e., its influence was scaled to be one) when the mean distance between forest locations is 40 km (Fig. 5). Preceding-year defoliation has 13 and 61% more weight from a forest location itself than its neighbors 40 and 80 km away, respectively, in predicting current-year defoliation at this location, while defoliation from more than 1 year ago has minimal effects on current-year defoliation.

Fig. 5
figure 5

Relative importance of predictors in our model (at their mean values except ΔE is at preset values shown in the figure; with D as a reference of an absolute value of one; negative values indicate effects of reducing defoliation; Achen 1982) at two levels of mean distance between forest locations, where Def is the defoliation (%), L is the metric of habitat suitability (0-1 ratio), C is the host tree percentage (%), W is the metric of wind (m s-1), and H is the dominant height

Full size image

Habitat suitability is the second most influential factor (at distances shown in Fig. 5), while host tree percentage is the most influential local factor (intrinsic stand characteristic) of those affecting SBW outbreak dynamics evaluated in this study (Fig. 5). However, influences of landscape factors (e.g., wind and habitat suitability) decrease as distances between forest locations increase, while effects of local factors remain the same (Fig. 5). Increases in habitat suitability and relative stand density and decreases in elevation from source to destination locations of defoliation reduce the level of defoliation at the destination location (Table 3). However, the above effects only slightly delay, but not prevent defoliation, across the landscape through the course of an outbreak (Fig. 4). On the other hand, increases in defoliation in the landscape, host tree percentage, wind, stand dominant height, and differences in elevation are all positively related to subsequent intensification of an SBW outbreak (Table 3).

3.3 Simulations of spatial and temporal dynamics of an SBW outbreak

Host tree percentage and wind velocity have noticeably greater effects than the other factors on the spread and intensification of SBW defoliation across the landscape in simulated outbreaks presented in Fig. 4. In comparison, changes in relative stand density and wind direction only slightly modify how an outbreak develops compared to the reference scenario described in Table 2. SBW defoliation spreads over long distances (100+ km in cases) and results in relatively distinctive spatial patterns of defoliation in the first 2 years of simulated outbreaks in various scenarios.

However, defoliation generally becomes ubiquitous in 3 years despite it is initiated at only four locations in all scenarios (Fig. 4). The percentages of areas defoliated in each scenario (Table 2) in each of the 3 years after the initiation of an outbreak are shown in Table 4. The significant differences in various factors’ relative importance shown in Fig. 5 are mostly reflected in the severity (levels) rather than the presence of defoliation across the landscape as an SBW outbreak continuous to develop (Fig. 4). In addition, simulated defoliation levels are generally low, which is characteristic of the last SBW outbreak in Maine.

Table 4 Percentages of areas defoliated in each scenario (Table 2) in each of the 3 years after the initiation of an SBW outbreak
Full size table

4 Discussion

In this study, we developed a flexible spatiotemporal parametric model to explicitly account for various factors that can influence SBW defoliation in continuous space and time in order to evaluate the dynamics of SBW outbreaks across a complex landscape. The accuracy of this model is an improvement from the basic categorical predictions of low/medium/high or the presence/absence of defoliation (e.g., Hennigar et al. 2013; Rahimzadeh-Bajgiran et al. 2018). Simulations based on this model show that defoliation generally becomes ubiquitous in 3 years across an area of 160 × 240 km despite varying environmental and stand conditions over large ranges. Our defoliation observations fluctuate greatly over time (Chen et al. 2018a), while our model, specifically the estimate of the lag parameter dt (Fig. 3), indicates that defoliation has no real long-term memory such that current-year defoliation has a rather limited correlation with defoliation more than 1 year ago. This may be related to the role of SBW dispersal in sustaining defoliation across space and time and also suggests that defoliation risk mapping based on stand characteristics probably cannot sufficiently reflect the dynamic nature of defoliation.

After SBW defoliation is initiated, its development across the landscape seems to be primarily facilitated by wind, in terms of both the direction and speed of dispersal (Fig. 4). This is consistent with previous observations and analyses (e.g., Greenbank et al. 1980; Blais 1983; Anderson and Sturtevant 2011) and is also expected considering that long-distance dispersal is fundamental to population dynamics of many insect defoliators and SBW moths have limited fly ability without wind (Isard and Gage 2001; Sturtevant et al. 2013). However, this effect of wind is greatly influenced by landscape structure. For example, Bouchard and Auger (2014) found that the SBW outbreak during the last decade progressed more slowly in western Quebec compared to the rest of the province. Although this direction of dispersal was against the prevailing wind, they attributed their observation to the high abundance of hardwood trees. While it is obvious that forests composed of fewer host trees (e.g., those dominated by hardwood trees) would be less defoliated, it is relatively uncertain whether these less susceptible forests would also lower defoliation of other forests across the landscape. Our analysis suggests otherwise as less defoliation spreads from one location to another under the same wind condition when the landscape in between is more connected (i.e., composed of more host trees). A possible cause may be that SBW moth flights are more likely to land in suitable habitats midway and fewer will reach the other location. Obviously, this effect of more connected and abundant host trees only slightly delays long-distance dispersal of defoliation and may result in prolonged outbreaks.

In addition to habitat suitability, elevation has also been thought to affect the development of SBW outbreaks (Bouchard and Auger 2014), and this influence has been attributed to gradients in, e.g., temperature, soil, and species composition, which are all generally correlated with elevation (e.g., Blais 1958; Magnussen et al. 2004). Our analysis indicates that it is the difference in elevation rather than elevation itself that affects the development of defoliation. For example, Chen et al. (2018a) found that elevation had minimal correlation with defoliation after the effects of associated tree and stand characteristics have been accounted for. Specifically, compared to elevation at the source location of defoliation, defoliation appears higher at locations higher in elevation and lower at locations lower in elevation when all of the other conditions remain the same. This effect may be due to that SBW moth flights are concentrated between 300 and 800 m above ground level at the top of temperature inversion zones, where maximum horizontal wind speed is often recorded (Greenbank et al. 1980; Boulanger et al. 2017). Therefore, moth flights are more likely intercepted by high rises in elevation. This could also explain the much higher defoliation in Baxter State Park (Osawa et al. 1986), which is at the geographic center (i.e., crossroad of SBW dispersal pathways), and higher in elevation, compared to the rest of the spruce-fir forests during the last SBW outbreak in Maine (Chen et al. 2017a).

Forest stand characteristics (local factors in our analysis) like species composition, dominant height, and relative stand density also affect the dynamics of SBW defoliation, specifically in determining sustainable levels of defoliation and possible dispersal of defoliation by the emigration of SBW. Findings in this study are consistent with previous studies, e.g., Hennigar et al. (2008) and Chen et al. (2018a). What is unexpected is that the metric of temperature tested in this study is practically not related to defoliation dynamics, despite that SBW moths’ exodus flights only take place within a certain temperature range (Royama 1984). Since Boulanger et al. (2017) observed SBW moth flights of 200 km in 4 hours, it suggests that SBW moths may be capable of dispersing over rather long distances within a small time window when temperature is within the suitable range. Therefore, the number of days within this temperature range may be largely irrelevant to the dispersal of SBW and consequent defoliation across the landscape. However, it is also possible that SBW phenology differs enough across a relatively large region of this study to result in moths taking flight at varying times of the year (Anderson and Sturtevant 2011). Consequently, a single metric of temperature is likely insufficient to capture the varying SBW dispersal and defoliation dynamics.

This study is with some important limitations and uncertainties, which are mainly in the following three aspects. First, our defoliation data do not cover the full temporal extent, especially the early stage, of the last SBW outbreak in Maine. Specifically, above normal SBW activity was observed since 1972 (Irland et al. 1988), but the defoliation data only start in 1975. Therefore, the outbreak probably developed slower than our study suggests. In addition, the outbreak in Maine was not isolated from the outbreak across a much larger area in North America. However, there lack comparable defoliation data from neighboring areas to be evaluated for their effects on the development of the outbreak in Maine, e.g., defoliation measurements were primarily obtained through the rather coarse aerial survey in Quebec, Canada and were only available since 1984 in New Brunswick, Canada (MacLean and Erdle 1986; Gray et al. 2000). This limited spatiotemporal extent of our data may have added uncertainties to estimations of parameters in our spatial and temporal kernel functions. Second, aerial spraying of insecticide was applied to parts of the spruce-fir forests in Maine each year during the last SBW outbreak (Seegrist and Arner 1982). Although it was found that spraying did not provide significant foliage protection (Fleming et al. 1984; Lysyk 1990) such that there was no noticeable differences in defoliation compared to trees not sprayed (Chen et al. 2018a), the spatial and temporal patterns of spraying may have distorted those of an SBW outbreak to an uncertain degree. Finally, SBW, forests, climate, and their interactions continue to evolve. For example, SBW’s host trees have declined from 69% in 1975 to 44% in 2017 of all trees in the spruce-fir forests in Maine (both by stem counts; Chen et al. 2019). However, this does not necessarily indicate less severe defoliation in the future. The reason may be due partly to that short life cycles make SBW evolve much faster than its host trees, especially considering it has been present in this region for up to 8000 years, while the forests having undergone significant changes (Lorimer 1977; Morin et al. 2007). In addition, uncertain futures of climate change probably will also interact with the evolution of SBW, its host trees, and the forests in Maine, which will add even more uncertainties to the dynamics of future SBW outbreaks (De Grandpré et al. 2019).

Despite these key limitations, this study provides certain implications for forest conservation and management. In light of the rapid and pervasive development of SBW outbreaks, it may be more efficient to apply mitigation practices like highly targeted insecticide spraying early to initial spots (epicenters) of defoliation. Otherwise, the strong dispersal ability likely will enable SBW spreading to other suitable habitats capable of sustaining large populations and reinitiating the dispersal process again. This practice of early intervention is currently being tested across broad areas in New Brunswick, Canada, with positive results in the near term (MacLean et al. 2019). Furthermore, instead of hoping to escape defoliation that will quickly become ubiquitous across the landscape, management in the midterm probably should focus on improving forests’ resilience to withstand repeated defoliation by, e.g., altering species composition, as Fig. 4 shows that defoliation spreads and intensifies more slowly when host tree percentage is down to 50%.

5 Conclusion

Overall, our analysis provides an alternative yet robust approach to evaluate the spatial and temporal dynamics of SBW outbreaks by explicitly predicting defoliation in continuous space and time across the landscape, while avoiding the highly contested assumptions of the many processes governing SBW population dynamics (Pureswaran et al. 2016). The quantitative information generated by our model is both flexible in spatial and temporal scales and directly usable in existing forest growth and yield modeling frameworks and management decision support systems, which are useful tools for assessing health, productivity, and succession of forests influenced by SBW defoliation. Our analysis is readily extendable to evaluating spatial and temporal dynamics of other forms of defoliation across forest landscapes given the general robustness, flexibility, and strong performance of the approach.