Topographic Variation in Forest Expansion Processes across a Mosaic Landscape in Western Canada

Changes to historic fire and grazing regimes have been associated with the expansion of tree cover at forest–grassland boundaries. We evaluated forest expansion across a mosaic landscape in western Canada using aerial photos, airborne laser scanning, and field transects. The annual rate of forest expansion (0.12%) was on the low end of rates documented across North America and was greater from the 1970s to the 1990s than from the 1990s to 2018. Most forest expansion occurred within 50 m of established forests, and 68% of all tree regeneration in grasslands was within 15 m of the forest edge. The intensity of cattle grazing did not affect the tree regeneration density. Despite the slow pace of land cover change, grassland areas near the forest edge had an average of 20% canopy cover and 9 m canopy height, indicating the presence of tall but sporadic trees. The rate of forest expansion, density of tree regeneration, and tree cover within grasslands were all greater at lower elevations where trembling aspen (Populus tremuloides) and white spruce (Picea glauca) were the dominant tree species. We conclude that proportions of forest–grassland cover on this landscape are not expected to change dramatically in the absence of major fire over the next several decades.


Introduction
Forest expansion into grasslands has been observed around the globe at varying rates [1][2][3][4]. Forest expansion is a subcategory of woody plant encroachment and can be defined as an increase in density, cover, or biomass of trees in vegetation types where they were not previously found [4][5][6]. Rates of woody plant encroachment into grasslands in North America range from <0.1 to 2.3% per year [5] and there is concern for the long-term implications of such land cover change. Grasslands are one of the most at-risk biomes for habitat loss with over half of the world's temperate grasslands already developed or converted to other land uses [7]. There is also uncertainty in how aggregate shifts from herbaceous to woody plant cover affect terrestrial carbon stocks [8]. A better understanding of the factors that influence forest expansion processes would help in projecting future land cover change as well as its effects on carbon stocks [5] and habitat availability.
Rates and patterns of forest expansion vary across ecoregions, and even within landscapes [6,9]. Mean annual precipitation broadly determines the proportional coverage of different vegetation types (e.g., woody and herbaceous species). Sankaran et al. [10] identified a precipitation threshold of 650 (+/− 134) mm/yr, where mosaic landscapes with annual precipitation levels greater than this value were climatically determined. Areas with lower precipitation were considered to be disturbance-driven, meaning fire, grazing, and browsing were required to maintain the presence of both trees and grasses. Moisture availability also influences tree establishment and development in semi-arid climates [11][12][13], and it has been hypothesized that moisture and temperature limitations help to maintain forest-grassland boundaries [14][15][16][17].

Study Area
The western portion of Cypress Hills Interprovincial Park encompasses a 348 km 2 area in southwestern Saskatchewan and southeastern Alberta, Canada (49 • 37 N, 110 • 0 W) ( Figure 1). This area ranges from 1100 to 1450 m in elevation, rising an average of 250 m above surrounding agricultural and range lands. The average annual temperature and precipitation are 3.3 • C and 600 mm, respectively [48]. The warmest month is July with an average daily temperature of 15.9 • C. June and July are the months with the highest rainfall (103 mm and 60 mm, respectively), and 42% of all precipitation comes from snowfall.    The Cypress Hills landscape is a mosaic of forests and grasslands. Three tree species dominate the park's forest stands: trembling aspen (Populus tremuloides Michx.) and white spruce (Picea glauca (Moench) Voss) stands occupy the valleys and lower slopes, and lodgepole pine (Pinus contorta Dougl.) stands are found at higher elevations. The two main grassland types are rough fescue-dominated grasslands on upper plateaus, and mixedgrass prairie (Stipa spp. L., Bouteloua gracillis (Kunth), Koeleria cristata Pers., Agropyron dasystachyum (Hook.)) on drier slopes and at lower elevations [49,50]. The oldest forest stands originated from fires in 1885 and 1889, which were the last major fires in the park. Bison were extirpated from the landscape in the early 1880s, and elk and moose were introduced in the 1930s and 1940s [49]. Permitted grazing by domestic cattle began in 1918 and four different stock associations presently hold grazing permits.

Aerial Photographs
We used historical (1976)(1977)(1978)(1979)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998) and contemporary (2018) aerial photographs to quantify land cover change. The number, resolution, and year of acquisition of aerial photo data varied between the two provinces ( Table 1). The 2018 tiles were provided in color and were already georeferenced and orthorectified. The historical photos were scanned, individual, panchromatic, raw photos from flight surveys. Georectification of those photos was done in ArcGIS (ArcGIS 10.8; Esri software) using a spline transformation with an average of 35 control points per image, evenly spaced across each photo. Orthorectification was carried out using the ArcGIS Block Adjustment toolbar and a 25-m resolution digital elevation model. The average georeferencing error in the alignment of forest edges for historical and contemporary orthomosaics was 3.7 m (range: 2.5-5.3 m). Varying greyscale and shadows in the historical imagery, which arose from complex topography and the presence of shrubland, produced several errors during attempts at both supervised pixel-based classification and object-based classification. The best land cover classification results were achieved with manual digitization, and so this was the methodology we used. Vector shapefiles with areas representing forest, grassland, and water land cover types were created from each orthomosaic ( Figure 1). We defined forests as being ≥0.5 ha in size and having ≥50% canopy cover. This size threshold is used by the Food and Agriculture Organization (FAO) of the United Nations, while a 50% cover threshold was easiest to assess visually in aerial photographs and during field work. Forest openings >0.5 ha were delineated to evaluate infilling, but smaller forest gaps were ignored. For consistency, forest areas that were subjected to harvesting or thinning were classified as forest in all images. We defined two time intervals for land cover change based on the years of aerial photo acquisition: interval 1 represented the period 1976-1998 in Alberta and1979-1991 in Saskatchewan (hereafter referred to as 1970s-1990s); interval 2 represented the period 1998-2018 in Alberta and 1991-2018 in Saskatchewan (referred to as 1990s-2018). We refer to areas that were forested at the start of interval 1 as pre-existing forests and areas that were forested at the end of interval 2 as present forests.
We rasterized the vector-based land classification at a 25-m resolution to match that of the digital elevation data available for the study area [51]. At this resolution, data were less sensitive to georeferencing errors (which were generally less than 5 m) and errors in precisely delineating the forest boundary in areas where it is diffuse. Cells with ≥50% forest cover were classified as forest, and cells with <50% forest cover as grassland. We subset the raster data before fitting statistical models of forest expansion by (1) discarding cells representing pre-existing forest or water, since these would not transition to forest cover through forest expansion; (2) discarding cells >2 km from the forest edge, since no cells were observed to transition to forest at this distance from pre-existing forest; and (3) retaining every eighth horizontal and vertical cell (i.e., 200 m separation) to better ensure spatial independence. We tested the latter assumption by fitting a preliminary binomial generalized linear model, with forest cover after interval 1 as the response and distance to present forest, elevation, and their two-way interaction as predictor variables. Residuals from this model showed very weak autocorrelation at a distance of 200 m (Moran's I = 0.02).

Airborne Laser Scanning
We collected airborne laser scanning data across approximately 75% of the study area in August 2015 ( Figure 1). We used these data to measure canopy height and cover, but did not use it as a source of topographic information because of incomplete coverage. Laser pulse returns with a point density of 0.5 pts/m 2 were acquired using an Optech Aquarius sensor operating at a wavelength of 532 nm. We used the lidR package in R [52] to classify the resulting point cloud into ground and non-ground returns, and to normalize the points to represent height above ground level. We then overlaid the point cloud with five land cover classes derived from the classified aerial imagery: (1) pre-existing forests; (2) areas that became forested during interval 1 (1970s-1990s); (3) areas that became forested during interval 2 (1990s-2018); (4) grassland areas within 50 m of the present forest edge; and (5) grassland areas that are more than 50 m from the present forest edge. Next, we grouped the point cloud into a set of spatial clusters in order to calculate canopy cover and canopy height for distinct areas and later model them against topographic predictor variables. To do this, we used a k-means clustering algorithm to assign all laser returns into one of 8357 spatial clusters based on their x-y coordinates. We chose the number of clusters in each land cover class so that each encompassed a mean area of 1 ha (n = 3219, 353, 179, 1242, and 3364 for the five classes defined above, respectively) and contained sufficient points to reliably estimate canopy cover and height. We then decimated the points to a uniform spatial density and computed the 95th percentile of non-ground return height (hereafter, upper canopy height), and the proportion of returns with height ≥2 m (hereafter, canopy cover) within each cluster.

Regeneration and Scat Transects
We measured tree regeneration density and grazing intensity outwards from the present forest edge in the summer of 2020. We selected 90 sites across the study area with balanced representation across topographic conditions and forest types (denoted by black circles in Figure 1), and established 1-3 belt transects at each site. Transect locations were recorded using a handheld GPS receiver (Garmin GPSMAP 64). Belt transects were 50 m in length, 4 m in width, originated at the present forest edge, and ran perpendicular to the forest into grassland (where forest was defined as having canopy cover ≥50% and canopy height ≥ 10 m). The number of small trees (diameter at breast height < 2.5 cm) of each tree species was counted in each of 10 contiguous 4 m × 5 m sections along the length of each transect. Aspen suckers <0.4 m in height were not counted, as they were less visible in tall grass. We also counted the total number of cow, elk, and deer scats within 0.5 m of the transect center line (50 m 2 per transect). We assumed that scat remained observable for a similar amount of time for each species and that the total density of scat reflected the overall ungulate activity around each transect.

Statistical Analysis
We developed generalized linear models and generalized linear mixed models to predict (1) the 20-year probability of transition from grassland to forest cover in each of our two time intervals; (2) upper canopy height and canopy cover of areas that had become forested in interval 1, interval 2, or that were presently adjacent (within 50 m of a forest edge) grasslands; and (3) the regeneration densities of three species (white spruce, trembling aspen, lodgepole pine) in adjacent grasslands ( Table 2). Land cover change was determined by whether a location that was grassland at the start of a given time period had become forested by the end of that time period (=1) or not (=0). The models (1-3) described how these response variables varied with elevation, slope, and aspect, as well as distance to the forest edge (1 and 3 only) and grazing activity (3 only). Table 1 provides an overview of the response and predictor variables, probability distributions, link functions, and sample sizes for each group of models. Further technical details on model formulation, selection, and Bayesian fitting can be found in Supplementary Materials.

Aerial Photographs Airborne Laser Scanning Transects
Response variables

1.
Raster cell transition to forest cover during interval 1 2.
Raster cell transition to forest cover during interval 2
Upper canopy height of forests from interval 2 3.
Upper canopy height in grassland within 50 m of the forest edge (present) 4.
Canopy cover of forests from interval 1 5.
Canopy cover of forests from interval 2 6.
Canopy cover in grassland within 50 m of the forest edge (present)
Trembling aspen regeneration density

Land Cover Change
Across the study area, observed land cover change showed that grassland area decreased by 3.7% (548 ha) during the 1970s-1990s time interval, and by another 3.6% (510 ha) during the 1990s-2018 interval. These changes correspond to 3.1% and 2.8% increases in forest cover, respectively. However, it is somewhat difficult to interpret these observed rates directly because of variation in the time interval lengths.
Both time periods had the same set of five best models (Table 3) for estimating land cover change. Model 9 performed especially well as the best model for interval 2 and second-best model for interval 1, so we retained it for further interpretation. Based on this model, which accounts for different time interval lengths, we estimate that 3.3% and 2.4% of grassland area was becoming forested every 20 years during interval 1 and 2, respectively ( Figure 2). Conversely, we estimate that forest cover was expanding by 2.8% and 1.9% every 20 years during these respective time periods. Table 3. Top five models for estimating land cover change in each time interval.

Land Cover Change Models for Time Interval 1 Land Cover Change Models for Time Interval 2
Model Notes: E = (elevation-1300)/100; S = (slope/5); EA (easterly aspect) = cos(90-aspect); NA (northerly aspect) = sin(90-aspect); G = (grazing-5)/5; ∆LOO and standard error (SE) represents the difference and standard error, respectively, in expected log predictive density estimates for a new dataset.  Models from both time periods indicate that the probability of transition to forest cover declines strongly with distance from the forest edge ( Figure 3). The probability of a grassland cell transitioning to forest in interval 1 was 10-22% (95% credible interval), 5-11%, and 1-4% at distances 25 m, 100 m, and 250 m, respectively, in areas with a 1 • slope, a 45 • aspect, and 1300 m elevation. For interval 2, the corresponding transition probabilities were 5-13%, 3-7%, and 1-4%.  Flat areas at lower elevations had the highest probability of becoming forested, whereas steep, SW-facing slopes had the lowest. For example, the models estimated that grassland cells located 50 m from the forest edge had a 9-20% and 5-11% probability of becoming forested in intervals 1 and 2, respectively, on NE-facing 1° slopes at 1250 m elevation. By contrast, the respective transition probabilities were only 2-5% and 1-4% for a SW-facing 10° slope at 1400 m elevation. Transition probabilities dropped more quickly with increasing distance from the forest edge on steeper slopes (Figure 3). Topographic effects were similar in the two time periods, although elevation had a somewhat more pronounced effect in interval 1 compared with interval 2 ( Figure 3A,B). Flat areas at lower elevations had the highest probability of becoming forested, whereas steep, SW-facing slopes had the lowest. For example, the models estimated that grassland cells located 50 m from the forest edge had a 9-20% and 5-11% probability of becoming forested in intervals 1 and 2, respectively, on NE-facing 1 • slopes at 1250 m elevation. By contrast, the respective transition probabilities were only 2-5% and 1-4% for Land 2021, 10, 1355 9 of 18 a SW-facing 10 • slope at 1400 m elevation. Transition probabilities dropped more quickly with increasing distance from the forest edge on steeper slopes (Figure 3). Topographic effects were similar in the two time periods, although elevation had a somewhat more pronounced effect in interval 1 compared with interval 2 ( Figure 3A,B).

Canopy Height and Cover
There was little difference between the upper canopy height of forest areas that established from the 1970s-1990s (mean ± s.d. = 13 ± 2 m) and those that established from the 1990s-2018 (12 ± 3 m) ( Figure 4A). The canopy height of these expanded forest areas was lower than that of pre-existing forests (18 ± 3 m). Canopy height in grassland areas within 50 m of the present forest edge was highly variable (9 ± 6 m), nearly spanning the range from present forests to open grasslands. Grasslands more than 50 m from the present forest edge had a mean canopy height of just 0.5 m (±2 m).
Areas that became forested in each interval had similar levels of canopy cover (39 ± 10% in interval 1; 40 ± 12% in interval 2) and were more open than pre-existing forests (canopy cover 50 ± 9%) ( Figure 4B). Grasslands directly adjacent to the forest edge had lower but still appreciable canopy cover (20 ± 16%), whereas grasslands more than 50 m from the forest edge had virtually no canopy (1 ± 4%). There was little difference between the upper canopy height of forest areas that established from the 1970s-1990s (mean ± s.d. = 13 ± 2 m) and those that established from the 1990s-2018 (12 ± 3 m) ( Figure 4A). The canopy height of these expanded forest areas was lower than that of pre-existing forests (18 ± 3 m). Canopy height in grassland areas within 50 m of the present forest edge was highly variable (9 ± 6 m), nearly spanning the range from present forests to open grasslands. Grasslands more than 50 m from the present forest edge had a mean canopy height of just 0.5 m (± 2 m).
Areas that became forested in each interval had similar levels of canopy cover (39 ± 10% in interval 1; 40 ± 12% in interval 2) and were more open than pre-existing forests (canopy cover 50 ± 9%) ( Figure 4B). Grasslands directly adjacent to the forest edge had lower but still appreciable canopy cover (20 ± 16%), whereas grasslands more than 50 m from the forest edge had virtually no canopy (1 ± 4%).  Table 4). The models explained almost half the variation in forests established between the 1970s-1990s but less variation was explained for the other two categories (R 2 = 0.03-0.47), and only a few consistent patterns emerged.
Areas into which forests had expanded during interval 1 were predicted to be less tall and have lower cover on steep slopes than in flat areas (12-13 m and 36-39% on a 10° slope vs. 13-15 m and 39-42% on a 1° slope). Newer forests that established during interval 2 were predicted to be taller and have greater canopy cover at lower elevations than higher ones (12-13 m and 42-48% at 1250 m vs. 9-11 m and 29-37% at 1450 m), and also had somewhat higher canopy cover on NE-facing slopes. Grassland areas within 50 m of the forest edge were predicted to be both taller and have greater canopy cover on steep, NE-facing slopes (10-12 m and 22-27% on a 10° NE slope vs. 7-9 m and 15-20% on a 2° SW slope) with canopy cover predicted to be slightly higher at lower elevations.  Table 4). The models explained almost half the variation in forests established between the 1970s-1990s but less variation was explained for the other two categories (R 2 = 0.03-0.47), and only a few consistent patterns emerged.
Areas into which forests had expanded during interval 1 were predicted to be less tall and have lower cover on steep slopes than in flat areas (12-13 m and 36-39% on a 10 • slope vs. 13-15 m and 39-42% on a 1 • slope). Newer forests that established during interval 2 were predicted to be taller and have greater canopy cover at lower elevations than higher ones (12-13 m and 42-48% at 1250 m vs. 9-11 m and 29-37% at 1450 m), and also had somewhat higher canopy cover on NE-facing slopes. Grassland areas within 50 m of the forest edge were predicted to be both taller and have greater canopy cover on steep, NE-facing slopes (10-12 m and 22-27% on a 10 • NE slope vs. 7-9 m and 15-20% on a 2 • W slope) with canopy cover predicted to be slightly higher at lower elevations.   Table 3 for symbol definitions and for ∆LOO and standard error (SE) definitions. Bayesian R 2 measures the goodness-of-fit for each model.

Regeneration Establishment
We recorded average regeneration densities of 132 stems/ha for white spruce, 1803 stems/ha for trembling aspen, and 101 stems/ha for lodgepole pine in grasslands within 50 m of the forest edge. The highest regeneration densities occurred next to the forest edge, with only 32% of all seedlings and suckers located further than 15 m into the grassland. While there was some variation in model preference between species, we retained the model that omitted grazing and a slope × aspect interaction for interpretation as its LOO values were within one standard error of the best model for each species (Table 5). Table 5. Top five models for estimating white spruce, trembling aspen, and lodgepole pine regeneration density.

White Spruce
Trembling Aspen Lodgepole Pine  Table 3 for symbol definitions and for ∆LOO and SE definitions.
The density of regeneration decreased with distance from the forest edge for all three species, but especially white spruce and trembling aspen (Figure 5). At a distance of 25 m, the predicted densities of white spruce, trembling aspen, and lodgepole pine were, on average, 6%, 3%, and 47% of their respective densities 5 m from the forest edge. With the exception of lodgepole pine, the majority of transects (76% of 180 transects for white spruce, 91% for trembling aspen, 34% for lodgepole pine) had predicted regeneration densities at 25 m that were less than 20% of the densities estimated at a 5 m distance. model that omitted grazing and a slope × aspect interaction for interpretation as its LOO values were within one standard error of the best model for each species (Table 5). The density of regeneration decreased with distance from the forest edge for all three species, but especially white spruce and trembling aspen (Figure 5). At a distance of 25 m, the predicted densities of white spruce, trembling aspen, and lodgepole pine were, on average, 6%, 3%, and 47% of their respective densities 5 m from the forest edge. With the exception of lodgepole pine, the majority of transects (76% of 180 transects for white spruce, 91% for trembling aspen, 34% for lodgepole pine) had predicted regeneration densities at 25 m that were less than 20% of the densities estimated at a 5 m distance. Table 5. Top five models for estimating white spruce, trembling aspen, and lodgepole pine regeneration density.

White Spruce
Trembling  Table 3 for symbol definitions and for ΔLOO and SE definitions. Figure 5. Effect of distance on regeneration abundance for individual transects. Note that the y-axis is on a logarithmic scale.Species-specific regeneration patterns generally followed the topographic distribution for each tree species ( Figure  6). White spruce regeneration was highest at low elevations and on steeper slopes with east and north aspects. Trembling aspen regeneration was highest at low elevations, but did not vary much with slope or aspect. Lodgepole pine regeneration was greatest at high elevations and on southerly aspects. The models showed no topographic effects on how far lodgepole pine or trembling aspen regeneration extended into the grasslands ( Figure 6). There was some indication that Figure 5. Effect of distance on regeneration abundance for individual transects. Note that the y-axis is on a logarithmic scale.Species-specific regeneration patterns generally followed the topographic distribution for each tree species (Figure 6). White spruce regeneration was highest at low elevations and on steeper slopes with east and north aspects. Trembling aspen regeneration was highest at low elevations, but did not vary much with slope or aspect. Lodgepole pine regeneration was greatest at high elevations and on southerly aspects. The models showed no topographic effects on how far lodgepole pine or trembling aspen regeneration extended into the grasslands ( Figure 6). There was some indication that white spruce regeneration did not drop off as strongly with distance at higher elevations and southerly aspects, but the predicted density of regeneration at 25-50 m from the edge was negligible.
within 50 m of the forest edge (1335 scats/ha; 76% of total scat). Wild ungulates active in these areas included deer (381 scats/ha, 22%) and elk (32 scats/ha, 2%). Regeneration models that included and excluded effects of grazing were both found among the top-ranked models for each tree species, with relatively little difference in expected predictive ability between them (ΔLOO = 0.0-0.6 for all species; Table 4). Thus, our results do not provide strong evidence for either positive or negative effects of grazing on tree regeneration in grassland areas. Based on the density of scat, cattle were the dominant grazer in grassland areas within 50 m of the forest edge (1335 scats/ha; 76% of total scat). Wild ungulates active in these areas included deer (381 scats/ha, 22%) and elk (32 scats/ha, 2%). Regeneration models that included and excluded effects of grazing were both found among the top-ranked models for each tree species, with relatively little difference in expected predictive ability between them (∆LOO = 0.0-0.6 for all species; Table 4). Thus, our results do not provide strong evidence for either positive or negative effects of grazing on tree regeneration in grassland areas.

Forest Expansion and Growth
Forest expansion has increased globally in recent decades [47], causing concern about the ecosystems into which trees are expanding and their potential impacts on ecosystem functions and services. Forests in our study area were estimated to expand by 2.8% and 1.9% over 20 years (approximately 0.12% per year), for time interval 1 and 2, respectively. This estimate of forest expansion is at the low end of woody plant expansion rates reported for North America (<0.1-2.3% per year) and for Africa, Australia, and South America (0.1-1.1%) [4,5], although our estimates do not include shrubs or isolated trees. There were corresponding 3.4% and 2.5% losses of grassland cover (approximately 0.15% per year) estimated for interval 1 and 2, respectively. Other studies in western North America have calculated rates of grassland loss to forest expansion approximately 2-6 times higher than we found [2,17,25,39].
An earlier study on the Alberta plateau of Cypress Hills found that tree cover expanded by approximately~1% per year between 1950 and 2002 [53]. Several factors could account for the greater expansion rate in that study. Their study area had a lower initial forest area (19%) compared with ours (54%), which would make relative changes in forest cover appear larger. It also spanned an earlier time period when the rate of land cover change may have been higher. Methodological differences could have played a role as well, as the authors used a lower canopy cover threshold than we did, and also included patches of trees <0.5 ha in area.
We found that forest expansion decelerated in recent decades, from 2.8% during the 1970s-1990s period to 1.9% from the 1990s-2018. This is likely part of a fluctuating pattern that has been observed in long-term studies [54,55]. The fact that the rate did not increase substantially in the second interval indicates that dramatic land cover change is unlikely in the coming years. Stevens et al. [4] found faster rates of woody plant expansion in more recent studies compared with older research and some have attributed this increasing rate in the second half of the 20th century to increasing concentrations of atmospheric CO 2 [56,57]. The decreased rate of expansion observed in this study area does not follow that pattern, however.
The mean upper canopy height of 20-year-old (established in interval 2) and 40-year-old forests (established in interval 1) only differed by about 1 m, and both were about 5.5 m shorter than pre-existing forests ( Figure 4). Pre-existing forest canopies in our study had approximately 10% more cover than the newer forests, compared with a 1% difference between the newer forests from the first and second intervals. There are a couple of possible explanations for similarity between 20-and 40-year-old forests. First, areas adjacent to the grasslands contained scattered adult trees of appreciable size. As trees fill in and a continuous forest canopy is established, the upper canopy height could still be determined by the older trees that first established in the grassland. Second, increased light availability along the edge may reduce allocation to vertical growth compared with the forest interior, where access to light may be limited by taller neighbors [31]. A third explanation is succession. As a somewhat shade-tolerant species, white spruce tends to grow slowly under an established canopy. Younger forests near the grassland edge would likely be dominated by trembling aspen or lodgepole pine, with white spruce not yet having reached the canopy or approached its maximum height of 25-30 m [58]. The upper canopy height may not change dramatically until after white spruce emerges above the surrounding canopy, a process that takes much longer than 40 years [59].

Spatial Patterns and Topographic Variation
Most forest expansion in our study area occurred within 25-50 m of the forest edge. Regeneration patterns also indicated similar spatial constraints with only 32% of all regen-Land 2021, 10, 1355 14 of 18 eration located further than 15 m into the grasslands. The strong effect of distance is likely the result of dispersal limitation and the facilitative effects of neighboring adult trees [35,60]. The increased moisture, more stable air temperatures, and decreased winds near the forest edge may be better suited to regeneration success than open grassland conditions [31,61]. Many grassland areas up to 50 m from the forest edge had appreciable tree cover (Figure 6), which suggests that forest expansion starts with sporadic tree establishment that, through the process of infilling, eventually reaches the 50% cover threshold we used to define forest.
Lodgepole pine and trembling aspen, which are shade-intolerant, pioneer species, may rely less on the facilitative effects of the forest edge. Of the three species in our study, lodgepole pine regeneration extended the farthest into grasslands. However, even though trembling aspen densities decreased more rapidly away from the forest edge, the proliferation of aspen suckers produced a similar regeneration density to lodgepole pine 25 m from the forest edge. Once established, these isolated pine trees and aspen colonies may have a nucleating effect [23], facilitating further tree regeneration in the harsher, open environment [34,35]. Isolated grassland trees also produce seed and suckers that contribute to infilling and the eventual formation of a closed forest canopy.
While dispersal from isolated trees in the grassland could lead to tree establishment even further from the forest edge, the slow rate of forest expansion observed over the last 40 years, and the low canopy height and virtual absence of tree cover in open grasslands (0.5 ± 2 m and 1 ± 4%), indicates that successful tree establishment more than 50 m from a forest edge is not common.
We observed higher probabilities of forest expansion on northeast-facing slopes, flat terrain, and at lower elevations. Flatter slopes and lower elevations were also associated with taller forest canopies and greater canopy cover for newly established forests, but topographic features did not affect how far regeneration extended into grasslands. Topographic patterns in regeneration followed the landscape-level distribution of forest types across the study area. While moisture availability at lower elevations may be a factor, the most likely explanation for topographic variation in forest expansion relates to species composition and species-specific regeneration mechanisms. For example, trembling aspen, which occupies lower elevations and shallower slopes, propagates through root suckers that produce regeneration densities an order of magnitude higher than the other species. Lodgepole pine, which reproduces by seed, had a lower regeneration density and dominates at upper elevations, where less expansion was observed. White spruce, which also reproduces by seed, had a similar overall regeneration density to lodgepole pine, but it is less likely to survive in open grasslands [29,62]. White spruce may regenerate better after pioneer species have established in the open, as the presence of other species has been shown to reduce its climatic sensitivity [11,63].

Grazing
No clear relationship between regeneration abundance and grazing intensity emerged in our results. Cattle do not eat many young trees as part of their diet [64]. Tramplingrelated mortality may sometimes occur [65], but not to an appreciable degree. It is possible that competition from grasses has only minor effects on tree regeneration and any beneficial effect from cattle grazing is small. Sankey et al. [16] and Fensham et al. [46] likewise found no overall trend of grazing on tree establishment, and Barnett and Stohlgren [66] found that greater densities of elk did not decrease aspen regeneration. An earlier study conducted on a portion of our study area also concluded that grazing did not substantially contribute to forest encroachment, despite a long-term negative correlation with tree establishment [53]. However, other studies have found effects of grazing on ecotone shifts [4,26,30].
There are several challenges to investigating long-term impacts of grazing on tree regeneration. Clear historical records on grazing management practices are often not available [47,66], and even when they are it can be difficult to separate grazing from other sources of environmental change [1,17,24]. While testing for a short-term relationship between grazer presence and regeneration density provides only a snapshot view, it could have revealed more localized effects if they were present.

Conclusions
This research quantified forest expansion across a mosaic landscape in western Canada. We estimated that rates of forest expansion were lower in the 1990s-2018 (1.9% every 20 years) compared with the 1970s-1990s (2.8% every 20 years) and that, overall, the annual rate of forest expansion (0.12%) was on the low end of woody plant expansion rates calculated for North and South America, Africa, and Australia [4,5]. Forest regeneration density decreased dramatically at distances greater than 15 m from the forest edge, and most changes in land cover over a twenty-year period occurred 25-50 m from an established forest edge. Overall, forest expansion occurred through slow, sporadic establishment of trees in grasslands adjacent to the forest edge that filled over time to establish new forest areas. These findings contribute to a better understanding of the rates and patterns of forest expansion worldwide.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/land10121355/s1, Table S1: Specifications for aerial photo data used to assess land cover change, Table S2. Formulations for the set of candidate models. Models 1-10 were considered for all response variables. Models 11-20, which included effects of grazing, were only considered for regeneration density, Table S3. Posterior estimates and credible interval values for the selected land cover change model (model 9 in Table S2), Table S4. Posterior estimates and credible interval values for selected models estimating upper canopy height and canopy cover for forests established in time interval 1 (model 3 in Table S2), forests established in time interval 2 (models 2 and 6 in Table S2), and tree cover in grasslands within 50 m of the forest edge (models 8 and 10 in Table S2),  Table S5. Posterior estimates and credible intervals for the selected regeneration density models (model 9 in Table S2) for white spruce, trembling aspen and lodgepole pine, Figure S1. Observed upper canopy height (95th percentile) and canopy cover of~1-ha areas by land cover category in relation to aspect, slope and elevation. Error bars represent standard error, Figure S2. Mean regeneration density (stems/ha) from field transects of white spruce (A-C), trembling aspen (D-F), and lodgepole pine (G-I) in relation to elevation, slope and aspect, Figure S3. Variation in activity of cattle and other ungulates, as measured by scat density, by elevation, slope, and aspect, Figure S4. Relationships between tree regeneration density in grassland areas within 50 m of the forest edge and grazing intensity. Topographic variables appeared to explain these patterns, as an independent grazing effect was not supported by our model selection proces.