The structure of boreal old-growth forests changes at multiple spatial scales over decades

Changes in the structure of boreal old-growth forests are typically studied at a specific spatial scale. Consequently, little is known about forest development across different spatial scales. We investigated how and at what spatial scales forest structure changed over several decades in three 4 km2 boreal old-growth forests landscapes in northeastern Finland and two in Quebec, Canada. We used canopy cover values visually interpreted to 0.1-ha grid cells from aerial photographs taken at three time points between the years 1959 and 2011, and error distributions quantified for the interpretation. We identified the spatial scales at which canopy cover changed between the time points, and examined the credibility of changes at these scales using the error distributions in Bayesian inference. Canopy cover changed at three to four spatial scales, the number of scales depending on the studied landscape and time interval. At large scales (15.4–321.7 ha), canopy cover increased in Finland during all time intervals. In Quebec, the direction of the large-scale change varied between the studied time intervals, owing to the occurrence of an insect outbreak and a consequent recovery. However, parts of these landscapes also showed canopy cover increase. Superimposed on the large-scale developments, canopy cover changed variably at smaller scales (1.3–2.8-ha and 0.1-ha). Our findings support the idea that the structure of boreal old-growth forests changes at discernible spatial scales. Instead of being driven by gap dynamics, the old-growth forests in the studied regions are currently reacting to large-scale drivers by an increase in canopy cover.


Introduction
Vast areas of boreal forests still remain outside of direct human influence (Gauthier et al. 2015). These forests play a crucial role in global biodiversity conservation (Thom and Seidl 2016), and terrestrial carbon cycling and storage (Bradshaw and Warkentin 2015). For efficient monitoring, conservation and management, it is important to understand the scaledependent dynamics of these forests with negligible human impact (Moussaoui et al. 2019). For example, the recent contradicting observations of tree biomass changes in the boreal zone (e.g., Girardin et al. 2016;Hellmann et al. 2016;Henttonen et al. 2017) could partly be related to varying scales of observation (Marchand et al. 2018), emphasizing the need of a better understanding of the scale-dependency of these changes (Estes et al. 2018).
Over a given time interval, the direction and magnitude of change in forest structure depend on the balance between three processes: the growth and mortality of existing trees, and recruitment of new ones. The drivers that influence these processes operate at various spatial and temporal scales. Consequently, changes in forest structure depend on the spatial and temporal scales of observation (Smith and Urban 1988). Traditionally, stand-replacing disturbances which caused tree mortality over large areas, particularly fire, were considered the most important processes that alter boreal forest structure (e.g., Zackrisson 1977;Bouchard et al. 2008;Wallenius 2011), while research during the past couple of decades has emphasized the role of small-scale disturbances (e.g., Kuuluvainen 1994;Pham et al. 2004;St-Denis et al. 2010). However, partial and patchy disturbances have been noted as important drivers of boreal forest dynamics (Bergeron and Fenton 2012;Girard et al. 2014;Kuuluvainen et al. 2014), and non-disturbance factors such as variation in topography may influence forest dynamics across spatial scales (Martin et al. 2018;Kulha et al. 2019). Consequently, understanding how boreal forests develop requires that forest dynamics is studied at multiple spatial and temporal scales, also at levels beyond the conventional gap-landscape -dichotomy.
Of the three processes that change boreal forest structure, scale-dependent drivers cause tree growth to increase tree biomass at varying rates within and between forest stands. Consequently, tree growth changes forest structure at different spatial scales within a forest landscape (Henry and Swan 1974;Martin et al. 2018;Moussaoui et al. 2019). For example, tree characteristics and the characteristics of the tree's neighborhood induce growth variation at within-stand scales , while changes in climate generally cause growth variation over large scales (Hellmann et al. 2016;Hofgaard et al. 2018). Certain factors, such as variation in soil water holding capacity or nutrient availability influence tree growth at a range of spatial scales (Hamel et al. 2004).
Similar to tree growth, but with an opposite effect on live tree biomass, tree mortality changes boreal forest structure from the scale of an individual tree up to the forest landscape (Kuuluvainen et al. 2014; Thom and Seidl 2016). The characteristics of tree mortality influence how forest structure changes. In time, mortality may occur as pulses (Bouchard and Pothier 2010) and/or continuously as background mortality (Kuuluvainen et al. 2014). In space, mortality can be similarly clustered, or random (Aakala et al. 2007). The prevalence of mortality events at different scales may vary between forest landscapes due to, e.g., differences in tree species composition, tree age and disturbance regime (Hennigar et al. 2008;Kuuluvainen and Aakala 2011;Girard et al. 2014).
Tree mortality creates openings in the canopy, enabling canopy recruitment in closed-canopy boreal forests where light availability is a limiting factor (Kuuluvainen 1994;Pham et al. 2004;Caron et al. 2009). Hence, the spatial scale at which recruitment occurs tends to depend on the sizes and spatial patterns of the mortality events. Similarly, the species and the number of individuals that ultimately reach the canopy largely depend on the characteristics of the opening (Gauthier et al. 2010). For example, small openings generally favor saplings of shade-tolerant species (Girard et al. 2014).
Compared to forests with closed canopy structure, the consequences of tree mortality for recruitment are more complex in open-canopy forests, such as the high-latitude boreal forests where light availability plays a smaller role on canopy recruitment (Hofgaard 1993). Instead, tree mortality influences recruitment by providing suitable regeneration microsites such as nurse logs or mineral soil exposed due to tree fall (Grenfell et al. 2011) and by decreasing belowground competition for water and nutrients (Kuuluvainen 1994). However, non-mortality related processes such as variation in ground microtopography (Hörnberg et al. 1997), the properties of ground vegetation, and the production of viable seeds (Zackrisson et al. 1995) also influence recruitment in these open-canopy forests.
Besides defining the direction and magnitude of how forest structure changes, the balance between tree growth, mortality and recruitment also dictates whether a forest landscape is in an equilibrium or non-equilibrium state (e.g., Smith and Urban 1988;Turner et al. 1993). Contrasting changes that are small relative to landscape area may create forest landscapes which are in a dynamic equilibrium state (i.e. shiftingmosaic steady state), where the overall state of the landscape remains constant over time, while a nonequilibrium landscape develops stochastically, for example due to disturbances that are large compared to the landscape area and the following succession (Baker 1989). However, forest structure changes at multiple spatial scales (Smith and Urban 1988). Compared to the equilibrium and non-equilibrium paradigms, the changes that occur at multiple spatial scales are better explained by the hierarchy theory and the hierarchical patch dynamics paradigm, which suggest that for this specific reason, forest landscapes consist of nested but discernible hierarchical patches at various successional stages (O'Neill et al. 1986;Wu and Loucks 1995;Wu and David 2002).
The equilibrium and non-equilibrium paradigms, and the hierarchical patch dynamic paradigm have been widely applied to characterize boreal forest dynamics. For example, the fire-created forest mosaics in Fennoscandia (Zackrisson 1977) and the old-growth forests of northeastern North America show characteristics of equilibrium landscapes and hierarchical patch dynamics (Gauthier et al. 2010), while forests where crown fires are common could be considered as non-equilibrium systems (e.g., Kafka et al. 2001). However, how forest dynamics appear depends on the spatial scale of observation and on the timing of observation period (Estes et al. 2018;Marchand et al. 2018). Therefore, such conceptual conclusions of forest dynamics necessitate the consideration of the spatial scale of observation and observation period (Smith and Urban 1988;Turner et al. 1993;Estes et al. 2018).
A major issue in addressing scale-dependent changes in forest structure has been the lack of suitable material that would enable the multiscale change analysis over large areas and long periods of time (Ohmann et al. 2014). While remote sensing provides records that cover large spatial extent and currently span decades, the lack of spatially explicit ground-truth values limits their usability in long-term change analysis (Lechner et al. 2012). Here, we overcame this limitation by using data from Kulha et al. (2018), where tree-ring based reconstructions of forest stand development were used to calibrate remotely sensed data from the past several decades and, importantly, to quantify the measurement error in these data. The calibrated time series and the scenespecific error estimates allowed us to examine changes in forest structure at multiple spatial scales, and to assess the credibility of changes at these scales over different periods of time, thus negating this major shortcoming in the analysis of archived remote sensing data. Incorporating the scale-dependency of forest development into change analyses allows for novel insights into the scales at which forest ecosystem changes should be examined, and into the processes that drive forest dynamics. Further, considering the scale dependency of forest development is required to assess whether forest landscapes are in equilibrium or non-equilibrium state.
In the current study, we examined the changes in boreal forest structure at different spatial scales over the past several decades. In accordance with the hierarchy theory and the hierarchical patch dynamics paradigm, we tested the idea that the structure of boreal old-growth forests changes at hierarchical spatial scales that are discernible. Specifically, we asked (1) whether we can identify specific scales at which forest structure changes, and if these scales are identifiable, (2) how the direction and magnitude of change varies among different scales, and (3) over different observation periods. Ultimately, this analysis intends to test the validity of small-scale gap dynamics as the main driver of forest dynamics, as is currently considered for boreal old-growth forests.

Study area
To test the premises derived from the hierarchy theory and hierarchical patch dynamics paradigm, we studied five boreal forest landscapes (2 km 9 2 km each) on two continents (Fig. S1). These landscapes differ, for example, in their tree species composition and disturbance history and thus enable testing the premises in various old-growth boreal forests. In northeastern Finland we examined three landscapes, hereafter denoted as Hirvaskangas, Pommituskukkulat and Hongikkovaara (67°44 0 N, 29°33 0 E). In the North Shore region of eastern Quebec, Canada, we examined two landscapes, denoted Lac Dionne and Pistuacanis (49°38 0 N, 67°55 0 W).
Of the study areas, northeastern Finland has a subcontinental climate. Mean annual temperature is -1°C, and the lowest monthly mean temperature is -13°C (highest ?13°C). The average annual precipitation sum is 570 mm. The North Shore region has a humid climate, with an annual mean temperature of 0°C. The mean temperatures for the coldest and warmest months are -18°C and ?14°C, respectively. The average annual precipitation sum is 1100 mm (climate data are averages from years 1970-2000Fick and Hijmans 2017).
The studied landscapes are characterized by mosaics of mineral soil forests, waterbodies, and forested and open peatlands. Unstratified glacial tills and gently rolling low mountains with treeless upper slopes are typical in northeastern Finland. Here, the elevation of the studied region ranges between 200 and 500 m above sea level (asl). In the North Shore region, soils in gentle slopes and depressions are mostly of unstratified glacial tills. Glaciofluvial sand deposits are common in valleys, as are rocky outcrops on moderate summits and slopes. Elevation in the Quebec study region ranges from 300 to 500 m asl.
Both studied regions have low tree species diversity. Pinus sylvestris (L.) dominates in Hirvaskangas, and Picea abies (L.) Karst, and Betula pubescens (Ehrh.) in Pommituskukkulat. No one particular tree species dominates in Hongikkovaara. Instead, all three abovementioned species occur variably within the landscape. The main tree species in Quebec are Picea mariana (Mill.) in Lac Dionne and Abies balsamea (L.) in Pistuacanis, A. balsamea also dominating in parts of the Lac Dionne landscape.

Visual interpretation, bias correction and interpretation error
We used the corrected canopy cover values from Kulha et al. (2018) to analyze changes in canopy cover over three time intervals in each studied landscape. We define canopy cover as the proportion of forest floor covered by the vertical projections of tree crowns. In short, the canopy cover values were visually interpreted for the forested parts of the landscapes (excluding, e.g., waterbodies), using stereopairs of aerial photographs taken at three time points and a grid of 0.1-ha cells (4096 0.1-ha cells per landscape; see Kulha et al. 2018 for full details on the used photographs, and their processing and interpretation). Henceforth, we call the photographs from the three time points the oldest, middlemost, and newest photographs (Table 1), and the time interval between the oldest and middlemost photographs the first time interval, the time interval between the middlemost and newest photographs the second time interval, and the time interval between the oldest and newest photographs the whole time interval.
Briefly, in Kulha et al. (2018), tree-ring based reconstructions of forest stand development in randomly selected interpretation grid cells (n = 66) were used to correct visual interpretation bias and to quantify random interpretation error. For this, the growth histories of all living and dead trees with a diameter of C 10 cm at a height of 1.3 within the selected cells were determined from tree-ring samples using standard dendrochronological methods . The growth histories were used to backcalculate the field-mapped tree sizes (Aakala et al. 2016) to correspond to the years the aerial photographs were taken. This was achieved by first back-calculating tree diameter change and then converting the change in diameter to change in crown size using the relationship between tree diameter and crown area. Dead trees were resurrected at their cross-dated year of death during the back-calculation. After resurrection their crown size changes were back-calculated similar to live trees, but assuming circular crowns. Last, the proportional canopy covers used to correct visual interpretation bias and to produce interpretation error distributions in regression modelling were calculated as the sum of non-overlapping crown areas within a cell divided by cell area.
We compiled the corrected canopy cover values from Kulha et al. (2018) into raster maps that show canopy cover in 0.1-ha cells at three time points in the studied landscapes. We subtracted the sequential maps to produce maps depicting canopy cover change in the studied landscapes over the first, second and whole time intervals. We denote these differences maps of canopy cover change.
Scales and spatial patterns of canopy cover change We used the maps of canopy cover change to examine how and at what spatial scales boreal old-growth forest structure has changed over the past several decades, and the distributions of visual interpretation error from Kulha et al. (2018) to assess the credibility of these changes. According to the hierarchy theory, the spatial scales at which the most salient changes occur are discernible and identifiable (O'Neill et al. 1986). Hence, we first aimed to identify such scales of canopy cover change using a method referred to as the scalederivative analysis (Pasanen et al. 2013). Then, for examining the changes and their credibilities within the identified scales, we used Bayesian scale space analysis for images (iBSiZer; Holmström and Pasanen 2012;Pasanen and Holmström 2015). In iBSiZer, we smoothed the maps of canopy cover change based on the identified scales, and used Bayesian inference for detecting credible changes in canopy cover at each identified scale. Smoothing helps to reveal changes in canopy cover at multiple spatial scales, as smoothing with a high smoothing level evens out the small-scale details, revealing locally average behavior (i.e. changes at large spatial scale), and smoothing with a low smoothing level maintains all but the smallestscale variation (i.e. changes at small spatial scale).
The scale-derivative analysis, utilized here for identifying the spatial scales at which the most salient changes in canopy cover occurred, is an objective approach that bases on the concept of 'scale-derivative'. Briefly (see Pasanen et al. 2013 for full details), scale-derivative is the derivative of the smooth with respect to the logarithm of the smoothing level, and in the scale-derivative analysis, the scales at which the most salient changes in canopy cover occurred are detected based on the smoothing levels that minimize the scale-derivative vector norm. For example, in a signal that is a sum of two components with different scales, the location of a local minimum represents a level at which the smaller scale is smoothed out, and the large-scale component, not yet affected by smoothing, is revealed. Here, we defined sequences of smoothing levels using such local minima of the scale-derivative norm for each map of canopy cover change. Henceforth, we call the identified local minima as scale breaks (sensu Wu 1999). We obtained patterns of canopy cover change within the identified scales by applying a Nadaraya-Watson smoother with a Gaussian kernel (e.g., Wand and Jones 1994) to each map of canopy cover change based on the scale break sequence identified for the particular map.
In scale-derivative analysis, we placed certain scale breaks manually by considering for weaker signs such as saddle points or changes in slope of the scalederivative norms, instead of automated identification. We verified that canopy cover changes occurred at multiple spatial scales by comparing the original scale-derivative norms with those of the permuted data (Fig. S2). Next, we used iBSiZer to evaluate the credibility of canopy cover change in each smoothed map of canopy cover change. For this, we first developed posterior predictive distributions for each map of canopy cover change by using the distributions of interpretation error by Kulha et al. (2018). We then approximated the posterior distributions of the smoothed maps of canopy cover change by producing a sample for the distributions. We drew a large sample from the posterior predictive distributions for each map of canopy cover change, subtracted the sampled maps of subsequent time points, and smoothed these subtractions according to the scale break sequence identified for the particular landscape and time interval. Then, we identified cells with credible canopy cover increase or decrease in each landscape, each time interval and at each identified scale using joint inference over all cells of a particular map (highest point-wise probabilities; Erästö and Holmström 2005;Holmström and Pasanen 2012). We flagged cells where the joint posterior probability exceeded a threshold of 95% as cells with credible canopy cover increase or decrease. The posterior means of canopy cover change for cells with the lowest credible canopy cover increase and decrease were reported as the credibility threshold values to illustrate the change in magnitude needed for credible change at different spatial scales. If no credible change occurred, credibility threshold was not reported.
Last, to produce comparable numerical information on the scales of canopy cover change, we estimated the characteristic spatial size of the canopy cover change at the identified scales using a combination of scalederivative analysis (Pasanen et al. 2013) and the diameter of the representative circle approach (cf. Pasanen et al. 2018). First, we decomposed the maps of canopy cover change into scale-dependent components using the identified scale break sequences (Holmström et al. 2011;Kulha et al. 2019). Second, we used the maximum of the scale-derivative vector norm of the scale-dependent components and the concept of full width at half maximum to estimate the diameter of a circular feature representative for the particular scale, and converted the diameter into area (cf. Pasanen et al. 2018). Within a scale, the first scale break represents the grain and the second the extent of the particular scale. Thus, features with a range of sizes exist within an identified scale, and the size estimation is of a feature representative for a particular scale.

Spatial scales at which canopy cover changes occurred
During the whole study interval (38 to 52 years, depending on the landscape), we identified canopy cover changes at three spatial scales in each landscape, based on the scale-derivative analysis. During the first time interval (16 to 32 years), we identified canopy cover changes at four spatial scales in Hirvaskangas and Pommituskukkulat, and at three spatial scales in the other three landscapes. During the second time interval (20 to 24 years), four spatial scales of canopy cover change were identified in all landscapes except Pommituskukkulat, where three scales were identified. Henceforth, we denote the identified spatial scales small, mid, large, and landscape scale. Small-, mid-, and large-scale changes were identified during all time intervals. When identified, landscape-scale changes only emphasized the large-scale patterns. Hence, they are only included in the supplementary material (Figs. S3-S4).
Representative feature size at the smallest identified scale was 0.1 ha (grain of the data) for each landscape and time interval. Mid-scale feature sizes were 1.3-2.8 ha (mean mid-scale feature size 2.0 ha; see Table S1 for full details). Feature sizes ranged from 15.4 to 125.6 ha (mean 41.6 ha), and from 180.9 to [ 321.7 ha (mean 272.6 ha) at the large, and landscape scales, quantified for all landscapes and time intervals, respectively.
For the whole time interval, the scale breaks between mid and large scales were automatically identified in Hirvaskangas, Hongikkovaara, and Lac Dionne. For the first time interval, the scale breaks between the mid and large scales for Hirvaskangas and Pistuacanis, and all the scale breaks between the large and landscape scales were detected automatically. For the second time interval, the analysis automatically identified the scale breaks between the mid and large scales in Hirvaskangas, Pommituskukkulat, and Pistuacanis, and between the large and the landscape scales in Hirvaskangas. Other scale breaks were placed manually at locations where saddle points or slope changes appeared in the scale-derivative norms. We verified the existence of the identified scales by comparing the scale-derivative norm of a canopy cover change map to the scale-derivative norms of corresponding permuted maps (Fig. S5). The only scale of canopy cover change that could be identified from the permuted data corresponded to small-scale changes in the original data (i.e. the 0.1-ha scale). The identification of only a single scale of change in the permuted data confirmed that canopy cover changes occurred at multiple spatial scales in the original data.
Canopy cover changes and their credibility at the identified spatial scales The rate at which annual canopy cover change occurred varied between the identified spatial scales (Fig. 1). Over the whole time interval (38 to 52 years in Finland), the Finnish landscapes showed mostly credible canopy cover increase at all scales, the increase being strongest at large scale (Figs. 1C, 2A1-A9, Table 2). Few individual small scale cells showed canopy cover decrease, whereas only positive changes were detected at the mid and large scales ( Fig. 2A1-A9). Consequently, 0% of the forested area showed a credible canopy cover decrease in the Finnish landscapes, independent of spatial scale, while the proportions of landscapes with canopy cover increase varied between 12 and 100% ( Fig. 2D1-D9, Table 2).
Over the whole time interval (46 years in Quebec), the Quebecois landscapes showed both canopy cover increase and decrease variably at the identified scales (Figs. 1, 3A1-A6, Table 3). Similar to Finland, canopy cover increased most at large scale, while both increase and decrease occurred more evenly at small and mid scales (Fig. 1C, Table 3). More credible canopy cover decrease occurred in Pistuacanis landscape than in Lac Dionne landscape, independent of spatial scale (Fig. 3D1-D6). The proportion of landscapes with canopy cover increase exceeded the proportion with canopy cover decrease in both landscapes (Table 3).
During the first (16 to 32 years in Finland) and second (21 to 23 years in Finland) time intervals, similar to whole time interval, the Finnish landscapes showed the most canopy cover increase at large scale and the least at the smallest identified scale (Figs. 1, 2B1-C9, Table 2). Of the individual landscapes, most credible canopy cover increase occurred in Pommituskukkulat during the first (Fig. 2E6) and in Hongikkovaara during the second time interval (Fig. 2F9). Contrary to the whole time interval, also credible canopy cover decrease occurred at small and mid scales in each Finnish landscape (Fig. 2E1-F8), Hirvaskangas landscape showing canopy cover decrease also at large scale (Fig. 2E3, F3).
During the first time interval (22 years in Quebec), annual canopy cover changes were mostly negative in Quebecois landscapes (Fig. 1, Table 3), and more negative in Pistuacanis than in Lac Dionne, independent of spatial scale (Figs. 1, 3B1-B6). Consequently, both landscapes show mostly credible canopy cover decrease at all identified scales during the first time interval. However, individual patches also showed canopy cover increase (Fig. 3E1-E6). Of the identified scales, larger proportion of the landscapes showed canopy cover decrease at large than at mid scale, and at mid scale than at small scale (Table 3). During the second time interval (24 years in Quebec), the annual canopy cover change was mostly positive in the Quebecois landscapes, independent of spatial scale (Figs. 1, 3C1-C6). Consequently, canopy cover credibly increased in these landscapes, while especially Pistuacanis landscape also showed individual patches where canopy cover credibly decreased (Fig. 3F1-F6). Similar to Finnish landscapes, larger proportion of the landscapes showed canopy cover increase at large than at mid, and at mid than at small scale (Table 3).

Discussion
We detected canopy cover changes at multiple spatial scales in boreal old-growth forests, using the scalederivative analysis (Pasanen et al. 2013). We identified changes at three to four hierarchical scales, the number of scales depending on the studied landscape and time interval. The nested hierarchical scales at which changes occurred were identifiable and decomposable. This is consistent with the hierarchy theory (O'Neill et al. 1986) and the hierarchical patch dynamics paradigm (Wu and Loucks 1995;Wu and David 2002), where landscape dynamics are seen as a sum of changes that occur hierarchically and in a nested manner at various spatial scales that are decomposable.
Considering the whole time interval (38 to 52 years, depending on the landscape), canopy cover changes occurred at three spatial scales. The credible large-scale changes indicated that canopy cover increased in the Finnish landscapes and in Lac Dionne, Quebec, while Pistuacanis landscape showed both credible increase and decrease in canopy cover. These large-scale patterns reflect the increase in mean canopy cover previously detected in the same four landscapes, and the zero mean change detected in Pistuacanis landscape (Kulha et al. 2018).   (Aakala et al. 2016), or by seedling establishment (Hofgaard 1993). In open-canopy boreal forests, factors such as the availability of seedbeds or the production of viable seeds (Zackrisson et al. 1995) and the availability of suitable regeneration microsites (e.g., fallen dead wood, exposed mineral soil; Grenfell et al. 2011) limit seedling establishment. As we studied changes in the canopy cover of the overstorey trees, the observed canopy cover increase being due to recruitment, in addition to growth of existing trees, requires that the established seedlings also successfully recruit to the canopy. The field data used in this study indicated continuous canopy recruitment in the Finnish landscapes during the whole time interval (Aakala 2018), and the current forest structure in the Finnish landscapes showed fewer signs of recent disturbances compared to the Quebecois landscapes (Kulha et al. 2019). For these reasons, we consider continuous tree growth and canopy recruitment together with limited tree mortality a plausible explanation for the observed large-scale canopy cover increase in the opencanopied Finnish landscapes.
At large scale, parts of the Hirvaskangas landscape in Finland showed credible canopy cover decrease during the first time interval . This indicates that due to a disturbance event, tree mortality locally surpassed the influence of tree growth and recruitment in this landscape. A number of similarly aligned trees were visible in the aerial photographs at this location, suggesting that this decrease in canopy cover was due to a storm. Hence, while disturbances directly acting as drivers of forest dynamics at these larger scales were rare in the Finnish landscapes, disturbances continue to influence these forests despite the near-complete elimination of wildfires from the region approximately a 100 years ago (Wallenius 2011;Aakala 2018). However, as the area influenced by the particular disturbance was small in relation to the landscape area, the disturbance negated the trendlike canopy cover increase only locally (cf. Baker 1989).
In Quebec, the role of disturbances clearly differed from the Finnish landscapes. Here, canopy cover at large scale predominantly decreased during the first  and increased during the second (1987-2011) time interval. This temporal pattern was consistent with the well-documented spruce budworm Choristoneura fumiferana (Clem.) outbreak that occurred in the region from the 1970s to the mid-1980s (Bouchard and Pothier 2010), and the consequent recovery of the forest after the end of the outbreak. In line with this inference, canopy cover decrease was more prominent in A. balsamea-dominated Pistuacanis than in P. mariana-dominated Lac Dionne landscape, as A. balsamea is more susceptible to the insect (Hennigar et al. 2008). In Lac Dionne, canopy cover decreased predominantly in A. balsamea-dominated parts of the landscape, while the open-canopied areas dominated by P. mariana showed canopy cover increase similar to Finnish landscapes over the whole time interval. In Quebec, tree species composition is reflective of disturbance history and soil characteristics, the overstorey proportion of A. balsamea being higher in mesic sites and typically increasing with stand age (Hamel et al. 2004;De Grandpré et al. 2000;Gauthier et al. 2010). Hence, the factors that influence tree species composition together with species-specific disturbances drive large scale changes in these forests. In addition to the changes at large spatial scales, we identified changes in forest structure at two smaller b Fig. 2 Annual canopy cover change in Finland during the whole time interval (column A), the first time interval (column B), and the second time interval (column C), and their credibilities in the respective order (columns D-F). Nonforested cells (e.g., waterbodies) appear as dark gray in all maps  The spatial scale at which the mid-scale changes occurred corresponds to patch-scale dynamics that have been identified from boreal forests (Kuuluvainen et al. 2014) and attributed to disturbances Bergeron and Fenton 2012). However, our results depict changes at these scales even in landscapes where patch-scale disturbances were mostly absent (Kulha et al. 2019). Further, changes at mid-scale characterized forest landscapes on both continents, despite differences in tree species composition and disturbance regimes (cf. Pham et al. 2004;Kuuluvainen and Aakala 2011). These findings supported the notion that boreal old-growth forest structure changes also at spatial scales other than the commonly proposed gap-and landscape-scales (Zackrisson 1977;Pham et al. 2004;Caron et al. 2009). At these scales, factors such as variation in soil properties and topography that cause differences in tree demographic rates, as well as recovery from past disturbances could partly explain the features detected here (Gauthier et al. 2010;Martin et al. 2018;Kulha et al. 2019).
The smallest scale at which we identified changes in forest structure were equal to the size of the grid cells (0.1 ha). Changes at this scale were evident in all landscapes, and during all time intervals. This scale corresponds to gap dynamics that are driven by mortality of individual or small groups of trees and the recovery from such events (e.g., Henry and Swan 1974;Caron et al. 2009;St-Denis et al. 2010). Our choice for the grain size of the data excluded the possibility to analyze within-stand dynamics. However, the significance of changes that occur at\ 0.1-ha scales is well-known in boreal forests (Kuuluvainen 1994;McCarthy 2001), and has been documented in both studied regions (Pham et al. 2004;Aakala et al. 2016).
The area with credible change at the 0.1-ha scale was low compared to changes at larger spatial scales. This could partly be related to the size of the grid cells. The magnitude of canopy cover change due to the death of an individual tree, or the infilling of a gap by an individual tree may not have been sufficiently large to surpass the visual interpretation error (Kulha et al. 2018). Furthermore, as smaller absolute change is needed to exceed credibility thresholds at larger than at smaller scales (Holmström and Pasanen 2012), the importance of changes at a particular scale cannot be directly deduced by comparing the number of credible cells at different scales (Kulha et al. 2018). However, changes in the smoothed maps of canopy cover change were stronger than changes in the permuted and smoothed maps of canopy cover change (Fig. S5), indicating that true large-scale changes were stronger than those randomly generated. While it is possible that the analyses underestimate the importance of gap dynamics, our results support the notion that largescale processes are currently driving the development of the old-growth forest structure in the studied regions.
Canopy cover increased and decreased variably at small scales, but either predominantly increased or decreased at the large scale, depending on the landscape. This suggests a scale dependency in the direction of how forest structure changes. Due to this scale dependency, differences in observation scales may cause differences in how forest ecosystem change appears (Estes et al. 2018;Marchand et al. 2018) and could underlie the contrasting changes in forest structure reported even from the same localities of the boreal region (Girardin et al. 2016;Hellmann et al. 2016;Henttonen et al. 2017). Our results emphasize that comprehensive understanding of forest dynamics necessitates that the scale of observation is explicitly considered when analyzing forest ecosystem change and the underlying drivers.
Because the direction and magnitude of change depend on both the observation period and the spatial scale of observation, whether a forest landscape can be considered an equilibrium or non-equilibrium system depends on scale and timing of observation (Smith and Urban 1988;Turner et al. 1993;Estes et al. 2018).
Here, at small scale the studied landscapes appeared either to be in equilibrium (e.g., canopy cover increase and decrease in Hongikkovaara during the first time interval) or in non-equilibrium (e.g., canopy cover decrease in Pistuacanis during the first time interval). However, at the large scale each landscape except for Pistuacanis appeared to be a non-equilibrium system where disturbances influenced a relatively small proportion of the landscape, consequently leading to directional landscape development (Baker 1989).
Certain scales of canopy cover change were not automatically identified in the scale-derivative analysis. Potential reasons for the difficulties in scale identification include notable patch size variability within a scale level and small patch size difference between the successive levels (Pasanen et al. 2013). This indicates that the scales at which ecological phenomena occur are difficult to detect (Scholes 2017). Still, to understand the behavior of forest ecosystems at the salient scales, identifying these scales should be targeted in forest change analysis. In particular, the identification of salient scales at which forest structure changes allows for speculation of the factors driving these changes (Wu and David 2002).
The Finnish landscapes and the P. mariana-dominated parts of Lac Dionne landscape showed synchronous large scale canopy cover increase over the whole time interval, potentially resulting from multiple driving processes. Recent research conducted in the same regions suggests that changes in forest structure after the late 19 th century could have resulted from changes in temperature and precipitation (Hofgaard et al. 2018;Sulla-Menashe et al. 2018). However, other factors such as prolonged fire return interval (Aakala 2018) and atmospheric nitrogen deposition (Henttonen et al. 2017) have also influenced forest growth and dynamics during the same time period. Consequently, attributing canopy cover increase to a particular driver is more complex than attributing canopy cover decrease to disturbances (Emmett et al. 2019).
Understanding the scale dependencies in forest development requires observations with broad spatial and temporal coverage typical for remotely sensed records. However, the usage of these records is limited by the lack of spatially explicit ground-truth values (Lechner et al. 2012). Here, the retrospectively produced ground-truth values enabled the analysis of these scale dependencies at spatial scales beyond those at which forest development is typically studied (i.e. plot, stand or landscape scale; Marchand et al. 2018). Another possibility for ground-truthing would be to integrate remotely sensed records with forest inventory data, given that issues that are related to scaling the different datasets permit such integration (Ohmann et al. 2014). We identified the characteristic scales at which forest structure changed using the scalederivative analysis that identifies these scales uniformly over the entire landscape (Pasanen et al. 2013). Obtaining calibration data limited our study extent to 4 km 2 , thus excluding potentially relevant changes in forest structure over larger areas due to, e.g., standreplacing fires that may occur in Quebec, and the following succession (De Grandpré et al. 2000;Bouchard et al. 2008). However, the approach can be widely applied to explore multiscale changes in any raster-form data and potentially at much larger scales, if suitable calibration data is available.

Conclusions
The exceptionally long time series of calibrated aerial photographs and the quantified visual interpretation error enabled the use of Bayesian inference in multiscale change analysis of boreal old-growth forests. The decadal-scale analysis revealed that the structure of boreal old-growth forests changed at hierarchical spatial scales that were discernible, and indicated scale dependencies in how forest structure changed. As expected from the current theory of boreal old-growth forest dynamics, we detected changes at scales corresponding to gap and patch dynamics, and discovered large scale changes plausibly caused by episodic disturbances in regions where such disturbances are prevalent. Instead, areas with a minor influence from larger-scale disturbances showed a trend-like increase in canopy cover that contrasts with expectations from this disturbancedriven system, and demonstrates the potential of topdown drivers (such as climate warming) in currently driving the development of boreal old-growth forest ecosystems.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Author contributions TA developed the original idea and study design, and NK, LP, TK, LD and TA refined the approach. NK interpreted the aerial photographs, and TA, TK, and LD collected the field data. LP and LH developed the analysis methods, and NK and LP conducted the analyses. NK wrote the first draft of the manuscript. All authors contributed to writing the final version.