Assessment of Past Decadal Dynamics of Tree Stands in Forest–Tundra Transition Zone on the Polar Ural Mountains Calibrated Using Historical and Modern Field Measurements

Abstract

Altitudinal forest limits are typically climatically dependent, such that increasing temperatures connected to global warming are causing upslope shifts in treeline ecotones worldwide. However, at the local and regional levels, the degree of such a response is dependent on differences in climate, topography and soil features. In recent decades, attempts have been undertaken to estimate tree stand dynamics with remote sensing methods, but their resolution is still too coarse for a precise assessment of stand structural changes, and requires ground-truthing, which is not possible without historical data collected on a single-tree level. We used aerial photos (1962) and satellite images (2021) in combination with historical inventory data to investigate changes in open forest positions at different spatial scales at the eastern macroslope of the Polar Urals over the past 60 years. Additionally, obtained remote sensing data were validated on a single-slope level using tree crown size estimations. Our investigations showed that since 1960 up to present day, the total crown coverage increased from 6.9 to 22.1% within the test polygon. A highly spatially variable upslope advance in an open forest boundary was identified from 1.7 up to 7.1 m in altitude per decade. We revealed that the rate of tree stand transformations was to a great extent depended on the stand density in the 1960s, soil substrate type, moisture regime, slope aspect and inclination. Our results highlighted the necessity to consider the abovementioned factors when trying to predict climate-induced tree distributional responses in subarctic mountain regions.

1. Introduction

The climate has been changing since the end of the Little Ice Age, having induced structural transformations on different biosphere levels [1]. Because mountains have a high degree of topographic heterogeneity, in turn leading to great environmental variability between sites [2], plant responses to environmental change are not uniform across mountainous landscapes [3]. Trees in regions with humid climates are typically growth-limited due to low summer temperatures at their upper range edges [4], but recent studies discovered that treeline advances are not solely determined here through warming [5,6,7,8,9], but on a landscape level constrained abiotically through features of topography (aspects or slopes) [10,11,12,13], wind load [14], microclimate and soil properties [15,16,17], and biotically due to the demands of tree species and interactions with ground vegetation [12,18,19,20,21,22,23]. Therefore, at some mountain slopes, treelines that experience warming may show no change in position [24]. The nonuniformity of elevational woody vegetation responses to climate change in mountains warrants further investigation into the factors that can influence these responses [25,26,27,28,29]. Evaluations of the parameters of individual trees and hole stands in the abovementioned treeline studies were mainly based on field time-consuming measurements [30,31]. Under short summer conditions and poor transport accessibility to mountainous areas, it is difficult to depend on a large coverage of measurements and representativeness of the data obtained. In recent decades, attempts had been undertaken to estimate the structures of tree stands using remote sensing methods, because repeated aerial photography and satellite imagery offer an excellent resource for the investigation of treeline shifts over wide territories [32,33,34]. Their resolution (0.3–1 m/pixel), however, is still too coarse for a precise assessment of stand structural changes, and requires ground-truthing or pairing with forest inventory data to provide detailed information on tree stand dynamics and changes in treeline positions [35,36,37]. However, in many cases, highly accurate estimations of former stand structures cannot be achieved without historical data collected on a single-tree level in previous periods. Thus, to determine accurate predictions about future changes in tree stand structures and treeline positions, it is crucial to increase the precision of remote sensing estimation based on former and contemporary field measurements, and to examine how different factors on local and landscape levels mediate responses to climate change.

In our study, we estimated structural changes to tree stands on their upper distributional limit in the middle part of the eastern macroslope of the Polar Ural Mountains (Mt. Chernaya, Rai–Iz massif and their surroundings) based on field observations (1960–2020) on a permanent altitudinal transect [38,39,40], interpretations of historical (1962) aerial photos and present-time (2020) satellite images of a test polygon and a comparison of open forest positions in the 1960s and 2020s in entire subregions of the Polar Urals. Additionally, we analyzed how the stand structure and dynamics couple of edaphic conditions and landscape features. Our objectives were: (1) to estimate crown closures within a test polygon (southern and eastern slopes of hills with landmarks, 312 m a.s.l.) by interpreting aerial photos and satellite images taken in 1962 and 2020; (2) to validate obtained crown coverage values using field observation data from 1960, 2011 and 2020, and to deduce the position of the upper boundary of open forests within the test polygon; (3) to investigate the vegetation patterns, detect the upper open forest boundary and assess its advances at a subregional level; (4) to identify the role of topography, sediment mechanical structures and soil wetness for determining stand structures and dynamics. We assumed that treeline positions and behaviors coupled to a great extent in regard to variations in slope aspects and ground peculiarities.

2. Materials and Methods

2.1. Study Site

Our study was performed on the eastern macroslope of the Polar Urals close to Mt. Chernaya (gabbro, 66.8377° N, 65.3836° E, 1036 m a.s.l.) and Rai–Iz massif (peridotite, 66.9544° N, 65.3304° E, 1309 m a.s.l.) (Figure 1). There was a chain of gentle-sloped hills (300–460 m a.s.l.) extending along their southeastern slopes. The east-facing slopes of the Rai–Iz massif were framed by the Slantsevaya, Yar–Keu and Pour–Keu mountains (crystalline shales, 400–880 m a.s.l.). On the plains (date of weather station Salekhard, 15 m a.s.l.), the mean air temperature in January was −22.1 °C and in July it was +14.2 °C, but at the same time, in the mountains (Rai–Iz weather station, 895 m a.s.l.) it was –17.9 °C and + 8.5 °C, respectively. The average annual precipitation was 460 mm in the valleys and 881 mm on mountain tops, one third of which would fall in the summer time. The average thickness of snow cover in the valleys was approximately 0.8–1.0 m. The stable snow cover laid approximately 240 days. In summer, western and southwestern winds prevailed, with western and northwestern winds prevailing in winter.

On the plains surrounding the Polar Urals and in the mountain valleys, the Siberian larch (Larix sibirica Ledeb.) dominated in tree stands, but, in some valleys, the Siberian spruce (Picea obovata Ledeb.) increased its abundance, in some areas even forming almost pure spruce forests. On some slopes within the treeline ecotone, birch (Betula pubescens Ehrh. ssp. tortuosa (Ledeb.) Nyman) and alder (Duschekia alnobetula subsp. fruticosa (Rupr.) Raus) prevailed. The vegetation forming the mountain tundra and dominating the open areas between groups of larches in the treeline ecotone (180–350 m a.s.l.) consisted of various shrubs, dwarf shrubs (Betula nana L., Salix sp., Vaccinium sp., Empetrum nigrum L., Arctous alpine L. and Dryas sp.) and herbs (e.g., Polygonum bistorta L., Anemone narcissifolia L. subsp. biarmiensis (Juz.) Jalas, Carex bigelowii Torr. ex Schwein., Festuca ovina L. and Solidago lapponica With.).

2.2. Trees and Stands Data Sampling and Calculation

The monitoring altitudinal (180–260 m a.s.l.) transect within the treeline ecotone was set up in 1960 on the southeastern slope (4 degrees of inclination) of a hill with a height mark of 312.8 m a.s.l., close to Mount Chernaya (Polar Urals), and had a width of 40 to 80 m and a length of 860 m [38,39,40,41] (Figure 1). It mainly consisted of rows of (2 or 4) plots 20 × 20 m. Within each plot, all adult trees and saplings taller than 20 cm were mapped and their morphometric parameters (tree height and diameter of crown horizontal projection) were measured. A repeat inventory and morphometric parameter measurements of the alive trees were carried out in 1999, 2011 and 2020. The sum of crown projection areas for each plot was calculated.

2.3. Estimation of Crown Closure by Interpreting Aerial Photos and Satellite Images

To estimate the structural changes of tree stands on the landscape level since the 1960s, the southern and eastern slopes of the hill surrounding the abovementioned altitudinal transect were taken and divided into 16,000 square plots with 20 m a side (see “Test polygon” in Figure 1).

Estimates for the crown closure of tree stands in the 1960s were obtained using monochrome aerial photographs (APs) with a survey date of 27 July 1962 (Figure 2). The spatial resolution of the digitized AP was 1 m pix⁻¹. The pixel brightness was encoded with values from 0 to 255 units. It helped to localize and identify the boundaries of tree crowns by using their shadows. As the photos were taken in the morning, shadows were directed from east to west. Using the SAGA GIS 2.3.2 software, we detected shadow pixels with a brightness value of less than 40 units in the APs. These pixel groups were vectorized into a polygon layer. For the analysis, we used contours with a size of more than 1 m², and if their width was greater than the height. The analysis of the morphometry of trees measured on the permanent transect revealed that the height of such trees were more than 2 m.

Tree crown projections (TCPs) were modeled using the contours of the tree’s shadows. The maximum values for the X coordinate of the corresponding shadow contour were taken as the coordinates of the TCP position along the X axis and the Y coordinate of the geometric center of the corresponding shadow contour were taken as the Y coordinate. The height of the shadow was taken as the diameter of the tree crown. This value was used to build the circles simulating the TCPs of the TCP center (Figure 3).

Estimates of the crown closure in the year 2020 for sparse tree stands and open forests, where shadows of neighboring trees did not overlap one another, were built using summer images from the Google Maps web service (satellite imagery supplier Airbus Maxar Technologies, true color image, spatial resolution 0.5 m pix⁻¹, see Figure 2). We also found the contours of the tree shadows and determined the width of the shadows and the length of the shadows from south to north (Figure 3). The width of the shadow, in this case, was taken as the diameter of the tree crown.

The crown closure of the modern dense stands was measured with the winter satellite WorldView-01 images from 2007 if it was not possible to measure the tree shadows (Figure 4). In the SAGA GIS 2.3.2, the Gauss smoothing procedure watershed segmentation was performed and the “Seeds points” layer was created. According to the field measurements in 2020, the average maximum diameter for the tree crown projection was determined for all 20 × 20 m plots of the monitoring transect (5.1 m). This value was used to build the TCP at the seed point, since this method did not allow for the depiction of all the trees in the dense areas (Figure 4). The share of areas where the crown closure was determined using winter images was 29%. The combined use of the two methods for modeling the TCP allowed us to assess the crown closure across the study area (16,000 square plots).

2.4. Calibration of Remote Sensing Data Using Field Observations

To increase the precision of crown closure data obtained through the use of the above described methods, we calibrated them with the TCP calculated from the crown diameter measurements of trees with heights of more than 2 m on the altitudinal monitoring transect in 1960 and 2020. At first, the crown closure (CC) was calculated for all 140 plots (20 × 20 m) based on the sum of the crown projections (obtained during field measurements or using tree shadow contours), where the double counting of crowns overlapping areas was excepted (stand crown projection—SCP).

$$CC=SCP/400,$$

(1)

where: CC—crown closure of a plot; SCP—the sum of the crown projections, except their overlapping areas.

A comparison of the CC determined from the aerial images and with the field measurements showed a strong relationship between them (R² = 0.82, Figure S1). The following formula was used to refine the final values of the closure of the entire study area.

$$CCfin1962=1.0874\times CCsh1962+0.003,$$

(2)

where: CCfin1962—the refined density of tree stands in 1962; CCsh1962—the CC determined using the shadows on the AP 1962.

The CC determined from the 2007 winter images was compared with the 2011 field measurements for squares with closely spaced trees. The data of the measurements of the crown projection areas, according to the winter images, were in good agreement with the field data (R² = 0.66) (Figure S2).

The tree crown projections obtained from the shadows on the satellite images in 2020 and with the winter images in 2007 were combined into a vector layer and used to calculate the CC for the whole study area. Then, this dataset was compared with the materials of the field measurements taken in 2020 for 140 20 × 20 m squares (Figure S3). The measurement data from the satellite images correlated well with the field measurements (R² = 0.82), and could be used to calculate the density of the forest stands under these conditions.

According to Formula (3), for the dependence of the CC determined from the field data and measurement data on the satellite images, the measured CC values were calibrated and taken as true.

$$CCfin2020=1.2661\ast CCsh2020+0.0246,$$

(3)

where: CCfin2020—calibrated forest stand density for 2020; CCsh2020—density determined using satellite images from 2020.

Field data for 140 squares for 2011 and CC data, measured with images from 2020, allowed us to find a strong relationship between them (R² = 0.81). Using the dependency, presented in Figure S4, we calculated the corrected crown closure for each square plot in 2011.

$$CCfin2011=0.9888\ast CCsh2020+0.0232,$$

(4)

where: CCfin2011—calibrated forest stand density for 2011; CCsh2020—density determined using satellite images from 2020 and 2011.

Using the calibrated CC values for 1962, 2011 and 2020, for each square plot (20 × 20 m) of the 1600-square “grid”, the 21 categories of crown closure were determined. The categories were determined by dividing the CC values from 0 to 1 into equal intervals of 0.05. Category 0 corresponded to a CC of 0 and category 20 to a CC of 0.96–1. For each crown closure category in 1962 the number of plots which came to other categories in 2021 was determined.

2.5. Detection of Position of Upper Boundaries of Open Forests

The above “grid” was used to determine the average CC along the border of low-density forests, which we carried out in the course of the visual interpretation of the satellite images taken in 2020. The average CC for the plots located along this border was 0.1 (0.088) units. This value was used to find the position of the open forest boundary for 1962, 2011 and 2020.

The generalized upper border of the open forests in 1962 was determined using plots with a CC of 0.1 and higher. Vertices were extracted from these plots and connected with lines (Figure 5).

Lines longer than 60 m were removed. The remaining lines were combined into polygons. Polygons at a distance of 60 m or closer to the river’s shoreline formed the primary line of the forest boundary; later, small polygons were connected to it, at a distance of less than 60 m (for joining such polygons, squares with a CC of more than 0 were used). Thus, the consolidated line of the upper boundary of open forests for 1962 was identified in 2 stages. The generalized upper forest border for 2011 and 2020 was determined using the same method.

2.6. Creation of Landscape Map and Detection of Upper Open Forest Boundary

The delineation of different categories of landscapes was performed with multispectral satellite imagery. We used winter, spring and summer Seintinel-2 images. All images were downloaded from the Sentinel Hub web service. For further data processing, we used spectral bands 2, 3 and 4 of the winter images, NDVI rasters and all bands of summer–autumn images. The listed dataset enabled us to perform a supervised classification in SAGA GIS. Training areas were determined using field observation data as a polygonal vector layer. Thus, we selected 8 types of landscapes: dense larch forests, open larch forests, dwarf birch (Betula nana), alder shrubs, grass, tundra, gabbro rocks and alluvial deposits. The classification procedure was performed using the method of the minimum distance and distance threshold of 0.

The grid layer we obtained as a result of the classification was vectorized into a polygonal layer. Then, we performed a manual correction of the misclassified objects. Finally, we obtained the vector landscape map to the southeastern part of the Rai–Iz range (Figure 6).

On the basis of the landscape map, we obtained a modern generalized borderline of the open forests. We applied the method that we used to construct a generalized border of open forests according to the measurement of the canopy using shadows.

2.7. Estimation of Upper Open Forest Boundary Shift

The assessment of the upward shift in the upper boundary of open forests (closure 0.1–0.2) in the Polar Urals was carried out through a comparative analysis of historical topographic maps and modern satellite images. These data were combined with a digital elevation model with a resolution of 2 × 2 m in the geographic information system ArcGIS 10.8. (ESRI Inc., Redlands, CA, USA). The boundaries of open forest distribution in 1960 were manually vectorized with topographic maps, and the boundaries in 2020, obtained from satellite images (described in Section 2.6), were added. In accordance with the categories of landscapes (Figure 6), the boundaries were divided into segments with the slight and strong influence of edaphic constrains. In addition, the studied area was divided into 5 parts with relatively similar topographical features (slope inclination and exposure). There were, mainly, gentle eastern and southern slopes in parts #1 and 2, but more steep ones in parts #3–5, in which the prevailing exposition varied from the northern and northeastern to southwestern.

An estimate of the area of the open forest expansion was obtained using the Spatial Statistics tools. The analysis did not include treeless sites within the area of the open forest expansion, where the establishment of forest stands did not occur due to the strong influence of edaphic constrains (bogging, lack of soil substrate, etc.). The calculation of the shift in the upper boundary of the open forests in a horizontal distance was carried out using the function of estimating the Euclidean distance between the borders at initial and final periods. The shift in the upper boundary of the open forests in altitude was estimated from the values (min, max, mean) of the cells of the digital elevation model corresponding to the position of the open forests.

The areas of the open forest expansion were divided into sections according to the different gradations of exposure and slope steepness. The slope steepness was split into classes in steps of 10°. Based on the enlarged raster (10 × 10 m), the exposure consisted of eight groups, 45° each: northern (N, 337.5–22.5°), northeastern (NE, 22.5–67.5°), etc.

3. Results

3.1. Tree Stand Dynamics within the Test Polygon Close to Mt. Chernaya

A comparison of crown closures reconstructed on base of aerial photographs from 1962 and satellite images from 2020 and validated using field measurements showed that there was a significant increase in the total tree crown coverage (from 6.9 to 22.1%) on the southern and eastern slopes of a hill with a height mark of 312 m a.s.l., close to Mt. Chernaya (Polar Urals) as a result of tree stand expansions in the tundra areas (Figure 7).

The share of plots without adult trees (>2 m) almost halved (from 65.4% to 34.5%) over the 60-year period as a result of transformations, mostly in sparse tree stands (crown closures of less than 0.1) and open forests (crown closures of 0.11–0.2) (see

Table 1. Share (in %) of the total number of plots with different crown densities in 1962 and 2020 on the southern and eastern slopes of a hill with a height mark of 312 m a.s.l., near Mount Chernaya (Polar Urals).

Crown closure	0.00	0.01–0.05	0.05–0.1	0.11–0.15	0.16–0.2	0.21–0.25	0.26–0.3	0.31–0.35	0.36–0.4	0.41–0.45	0.46–0.5	0.51–0.55	0.56–0.6	0.61–0.65	0.66–0.7	0.71–0.75	0.76–0.8	0.81–0.85	0.86–0.9	0.91–0.95	0.96–1.0
1962	65.5	18.6	7.9	4.4	1.9	1.0	0.5	0.2	0	0	0	0	0	0	0	0	0	0	0	0	0
2020	34.5	16.4	9.4	7.1	4.7	5.7	5.7	5.1	3.6	2.9	1.6	1.4	1.0	0.5	0.2	0.1	0.02	0.04	0.02	0.01	0.01

). The total percentage of plots with sparse tree stands remained almost unchanged (26.5 vs. 25.8%), but crown closures within a large part of such plots increased since the 1960s, and transformed in the open and closed forests by present day, while a similar quantity of tundra plots was settled with single trees. At the same time, the quantity of plots with open forests increased from 6.3 to 11.3%, and with closed forests in 17 times, from 1.6 to 27.8% (including stands with crown closures of 0.21–0.3 from 1.5 to 11.4%, 0.31–0.4 from 0.1 to 8.7% and 0.41–0.75 from 0 to 7.7%).

An analysis of changes in crown closures on individual plots since 1962 revealed that most of the tree stands with high crown closures (>0.70) at present time were formed on the sites with insignificant (0.01–0.2) initial crown closures in 1962, while on the plots with closed forests (>0.2) in the preliminary period, this increase was minimal (see

Table 2. Matrix of transition (in %) of plots with different crown closures in 1962 into plots with different crown closures in 2020 on southern and eastern slopes of a hill with a height mark of 312 m a.s.l., close to Mt. Chernaya (Polar Urals).

		2020
	Crown closure	0.00	0.01–0.05	0.05–0.1	0.11–0.15	0.16–0.2	0.21–0.25	0.26–0.3	0.31–0.35	0.36–0.4	0.41–0.45	0.46–0.5	0.51–0.55	0.56–0.6	0.61–0.65	0.66–0.7	0.71–0.75	0.76–0.8	0.81–0.85	0.86–0.9	0.91–0.95	0.96–1.0
1962	0.00	52.1	22.2	8.9	5.1	3.1	2.2	2.1	1.4	1.0	0.7	0.4	0.3	0.2	0.2	0.1	0.0	0.0	0.0	0.0	0.0	0.0
	0.01–0.05	0.0	10.1	14.4	13.5	9.0	11.3	11.1	9.5	6.9	5.2	3.2	2.0	1.9	1.1	0.5	0.1	0.1	0.0	0.0	0.1	0.0
	0.05–0.1	0.0	0.0	11.5	8.4	7.6	12.1	15.0	15.3	9.6	7.6	4.1	4.6	1.9	0.6	0.3	0.5	0.0	0.2	0.5	0.2	0.0
	0.11–0.15	0.0	0.0	0.0	12.5	5.7	16.7	13.1	14.5	11.4	9.6	3.7	4.5	4.3	1.1	1.1	1.1	0.0	0.3	0.0	0.0	0.3
	0.16–0.2	0.0	0.0	0.0	0.0	6.4	16.1	14.1	14.7	15.4	11.5	6.4	7.1	3.2	3.2	1.3	0.6	0.0	0.0	0.0	0.0	0.0
	0.21–0.25	0.0	0.0	0.0	0.0	0.0	17.1	11.0	23.2	8.6	12.2	8.5	6.1	8.5	2.4	1.2	0.0	0.0	1.2	0.0	0.0	0.0
	0.26–0.3	0.0	0.0	0.0	0.0	0.0	0.0	27.0	13.5	10.9	18.9	8.1	8.1	8.1	5.4	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	0.31–0.35	0.0	0.0	0.0	0.0	0.0	0.0	0.0	25.0	25.0	25.0	0.0	16.7	8.3	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

).

An analysis of the general dynamics of tree stands of various densities within the test polygon showed that, with an increase in the total tree crown covers of this territory from 1962 to 2020, there was upward shift in the upper boundary of the distribution of all stand categories, and, in particular, open forests (closure 0.1–0.2) by an average of 23.3 m in altitude (see Figure 8). Minimal values (10–14 m) were noted for a sector of the southern slope with moraine deposits (see S_m in Figure 8 and

Table 3. Altitudinal position of open forest boundaries (mean ± SD) on the southern and eastern slopes of a hill with height mark of 312 m a.s.l., close to Mt. Chernaya (Polar Urals), in 1962, 2011 2020, and their vertical shift from 1962 to 2020.

Sector Name	Altitudinal Position (m a.s.l.)			Altitudinal Difference (m a.s.l.)			Altitudinal Shift (m 10 y−1)
	1962	2011	2020	1962 and 2011	2011 and 2020	1962 and 2020	in 1962–2011	in 2011–2020	in 1962–2020
S	220 ± 17.8	243 ± 12.8	245 ± 10.9	23	3	25	4.6	3.0	4.4
SSE1	195 ± 6.5	234 ± 7.7	236 ± 7.3	40	1	41	8.1	1.4	7.1
SEE2	182 ± 2.2	208 ± 9.0	214 ± 10.3	26	5	32	5.4	6.0	5.5
E	158 ± 9.2	172 ± 6.9	186 ± 7.9	13	14	28	2.7	15.8	4.7
Mean for group	189 ± 9	214 ± 9	220 ± 9	26	6	32	5.2	6.6	5.4
S_m	198 ± 16.6	210 ± 19.6	213 ± 21.7	12	3	14	2.4	3.2	2.5
SE-1	176 ± 5.4	185 ± 11.2	189 ± 7.8	9	4	13	1.9	4.3	2.3
SE-2	167 ± 4.7	168 ± 6.8	177 ± 2.2	2	9	10	0.3	9.6	1.7
Mean for group	180 ± 9	188 ± 13	193 ± 11	8	5	12	1.5	5.7	2.2
Mean for all polygons	185 ± 9	203 ± 11	209 ± 10	18	6	23	3.6	6.2	4.0

) and southeastern sectors (SE1 and SE2) with waterlogged soils, whereas the maximum value (41 m) was marked for the south–southeastern sector (SSE-1). In sectors S, SSE-2 and E, the general shift was 25–32 m, although, in the first sectors, changes occurred mainly between 1962 and 2011; in the eastern sector, however, the changes were slow until 2011, and accelerated in the last decade (see Table 3).

Due to the fact that the inclination of the slopes in the sectors was different, the horizontal shifts (see

Table 4. Horizontal shifts (m) in open forest borders (mean ± SD) on the southern and eastern slopes of a hill with height mark of 312 m a.s.l., close to Mt. Chernaya (Polar Urals), in 1962–2020.

Sector Name	Absolute Horizontal Shift (m)			Horizontal Shift per 10 Years (m)
	1962–2011	2011–2020	1962–2020	1962–2011	2011–2020	1962–2020
S	152 ± 135	21 ± 30	167 ± 132	30.9	23.8	28.8
SSE1	344 ± 43	6 ± 11	356 ± 47	70.1	6.4	61.4
SEE2	260 ± 56	1 ± 4	268 ± 53	53.0	1.3	46.2
E	131 ± 112	194 ± 160	478 ± 144	26.6	215.9	82.3
S_m	62 ± 57	9 ± 19	67 ± 59	12.6	10.3	11.6
SE-1	126 ± 86	45 ± 74	193 ± 117	25.7	50.2	33.3
SE-2	36 ± 47	229 ± 188	254 ± 145	7.4	254.2	43.8

) had a slightly different sequence in their magnitudes than the vertical ones. Thus, a maximum horizontal shift was noted in the eastern (E) sector (478 m), and the second, using magnitude, was SSE-1 (356 m). The minimum values were noted in the same southern sector with a prevalence of moraine deposits (67 m).

Finally, a comparison of the shifts in the open forest boundaries in different parts of the test polygon demonstrated that such processes were 2.5 times slower on moraine deposits and waterlogged soils (see in Table 3, 4 sectors S_m, SE-1, SE-2) than on sites without them.

3.2. Characteristics of Open Forest Boundary and Its Dynamics at Studied Subregion of Polar Urals

An analysis of the distribution of moraine deposits, waterlogged soils and boulder fields at the upper location of the open forest boundary showed that their proportion varied significantly (from 0 to 0.64 of the total) at different parts of the studied subregion of the Polar Urals (see

Table 5. Characteristics of open forest boundary in the studied subregion of Polar Urals in 2020.

Part of Study Area	Portion of Total Length	Altitude in 2020 (m asl)		Distribution of Slope Inclination (Degrees)			Distribution of Slope Exposition
		Min	Max	0–10	10–20	20–30	N	NE	E	SE	S	SW	W	NW
Slopes without wide spread of moraine deposits and boulder fields
1	1	146	279	0.97	0.03	0.00	0.01	0.05	0.10	0.20	0.34	0.23	0.05	0.02
2	0.66	136	330	0.92	0.07	0.01	0.04	0.13	0.18	0.24	0.27	0.11	0.02	0.01
3	0.36	85	311	0.66	0.29	0.05	0.16	0.23	0.24	0.12	0.07	0.09	0.05	0.04
4	0.82	92	243	0.37	0.48	0.12	0.04	0.22	0.43	0.15	0.07	0.06	0.02	0.01
5	1	97	288	0.42	0.35	0.16	0.03	0.04	0.08	0.10	0.17	0.24	0.24	0.10
Slopes with a prevalence of moraine deposits and boulder fields
1	0	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.
2	0.34	145	253	0.91	0.08	0.01	0.06	0.23	0.18	0.18	0.20	0.11	0.02	0.02
3	0.64	73	235	0.76	0.20	0.03	0.15	0.30	0.25	0.10	0.05	0.05	0.04	0.06
4	0.18	140	228	0.53	0.36	0.08	0.08	0.17	0.11	0.13	0.21	0.12	0.10	0.08
5	0	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.

and Figure 9). On slopes with a prevalence of such types of edaphic conditions, the altitudinal position of the open forest boundary (min, max or mean) located lower (up to 108 m) than at slopes without their wide occurrence. The open forests ascended the maximal position (330 m a.s.l.) on the southern slope of the Rai–Iz massif and descended to minimal altitudes on the northeastern slope of the same massif. There were gentle slopes (mainly less than 10°) surrounding Mt. Chernaya (parts #1 and 2 of the subregion), but at the Rai–Iz massif and Sob’s river valley (parts #3–5), more abrupt slopes were allocated. The prevailing slope expositions at the locations of open forest boundaries changed from southern–southeastern in parts one and two to northeastern expositions at part three, eastern at part four and to southwestern expositions at part five of the studied subregion of the Polar Urals (see Table 5 and Figure 9).

An analysis of the tree stand dynamics over the 60-year period in the studied subregion of the Polar Urals showed that there was an upward shift in the upper boundary of the open forest (closure of 0.1–0.2) by an average of 33 m in altitude at slopes without waterlogged soils, moraine deposits or boulder fields (see Figure 9) and 9 m at slopes with a prevalence of them. Minimal values (0–11 m) were noted at the slopes with edaphic constrains on parts three and four (see Figure 9 and

Table 6. Altitudinal position and shift in the open forest boundary (mean ± SD) at different parts of the studied subregion of the Polar Urals in 1960 and 2020, and the area occupied by open forests over the 60-year period.

Part of Subregion	Altitudinal Position (m a.s.l.) in		Altitudinal Differences for Period 1960–2020 (m a.s.l.)	Shift in 1960–2020 (m 10 year−1)		Area Occupied by Open and Closed Forests (km2)
	1960	2020		Vertical	Horizontal
Slopes with slight influences of edaphic constrains
1	202 ± 28	244 ± 23	42	7.0	48.2	9.81
2	199 ± 39	231 ± 38	32	5.3	33.7	5.39
3	213 ± 61	238 ± 62	25	4.2	25.5	1.54
4	162 ± 39	200 ± 35	38	6.3	14.7	1.37
5	207 ± 39	238 ± 42	31	5.2	19.8	4.33
All	200 ± 42	233 ± 43	33	5.5	35.5	22.4
Slopes with strong influences of edaphic constrains
1	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.
2	194 ± 25	226 ± 29	32	5.3	41.3	3.59
3	119 ± 36	130 ± 42	11	1.8	13.8	1.71
4	169 ± 18	169 ± 20	0	0.0	3.2	0.03
5	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.
All	149 ± 47	158 ± 57	9	1.5	32.8	5.33

), whereas the maximum value (42 m) was marked in part one. A comparison of the mean altitude of the upper open forest boundaries on the slopes with different expositions within the Sob’s river valley (parts # 3–4 and 5) showed that open forests ascended higher on 33 m on SW–W than on NE–E slopes (

Table 7. Altitudinal position of the open forest boundary (mean and SD) on prevailing slope exposures at different parts of the studied subregion of Polar Urals in 1960 and 2020.

Part of Subregion	1		2		3		4		5
Year	1960	2020	1960	2020	1960	2020	1960	2020	1960	2020
Slopes slightly influence by edaphic constrains
Prevail slope exposure	S, SW		SE, S		NE, E		NE, E		SW, W
Mean	211	238	204	234	179	205	160	197	204	234
SD	43	24	36	35	42	60	36	36	42	46
Slopes strongly influence with edaphic constrains
Prevail slope exposure	n.d.		S, SW		NE, E		NE, S		n.d.
Mean	n.d.	n.d.	197	240	120	131	167	167	n.d.	n.d.
SD	n.d.	n.d.	32	34	36	42	14	16	n.d.	n.d.

). In general, the area occupied by open and closed forests increased in approximately 27.4 km² in this subregion of the Polar Urals over the 60-year period.

4. Discussion

Aerial photography and satellite images ground-truthed with plot-level forest inventory data demonstrates that at the eastern macroslope of the Polar Urals, the tree stand structures changed substantially over the 60-year period. A previous investigation [42] at that part of the Polar Urals showed that, between 1910 and 2000, the altitudinal shift in the upper boundary of open forests was 26 m (from 231 to 257 m a.s.l.) and closed forests was 35 m (from 195 to 230 m). Thus, the respective rates of altitudinal displacement over 90 years were 3 and 4 m per decade. In our research, we revealed for a coinciding part of the studied subregion of the Polar Urals (see part two on Figure 9) for the period between 1960 and 2020, upward changes in open forest positions were 32 m in the elevation (from 194–199 to 226–231 m) (see Table 5). Despite differences in our and previous investigations in time periods, methods of data requisition and definitions of forest type categories, obtained values were very close (4 or 5.3 m per decade).

We calculated altitudinal temperature lapse rates for the Polar Urals to be 0.65 °C per 100 m in altitude for the summer (see Section 2.1). Due to this temperature change with elevations, and the temperature increase in summer (0.82 °C) reported for the Polar Urals [30], the open forests line would be expected to have increased by 126 m (21 m per decade) in elevations between 1940–1960 and 2000–2020, if it was directly following the position of the isotherm. That was similar for the upward shift in isotherm lines due to climate warming estimated from different climatic records in the Alps, 10–20 m [43] or 20–40 m per decade [44]. However, the open forest position responses we provided for the Polar Urals were, therefore, significantly (four times) less than the advances predicted on the basis of climate change. Time-lags were identified between warming and tree stand structure changes, and so they may have taken some time before a response to an open forest advance could be recognized [15,45]. The disproportion between the movement of the isotherm and open forest lines indicate the complex interaction in ecological factors that could strongly influence the magnitude of the shifts on a landscape level [46]. As it was pointed out in some recent researches at the forest–tundra ecotone, this can be due to limited production, dispersal, and establishment of viable seeds [47,48,49,50,51,52,53,54], unfavorable soil properties [55,56,57] and microsite conditions [58,59,60] for seed germination and seedling survival, and competition of young trees with shrubs and heath species [61,62,63].

An analysis of the matrix of the transition of plots with different crown closures in 1962 (Table 2) supported statements about the constraint of forest advances due to insufficient numbers of viable seeds on the limit of tree distributions. Thus, our data demonstrated that tree closures increased only approximately half of the treeless plots since 1962 (see Table 2), while on plots with several single trees (crown closures of less than 0.05), it was observed in 90 % of cases. We assumed that the observed significant change in canopy closures in sparse stands and open forests (crown closures of 0.06–0.2) since 1962 was associated not only with an increase in the size of the tree crowns, but also with the intensification of seed productivity, which contributed to the active tree establishment in these stand categories. However, changes in density in closed stands (crown closures of 0.21–0.35) were not as significant since 1962, due to high competition for limited soil resources (nutrients) both between mature trees and between the undergrowth and mature trees.

Our findings also confirmed the results that the forest advance to the mountains could be significantly limited by unfavorable edaphic conditions, [55,56,57] or tree stands would be unable to occupy new territories due to unsuitable substrates [12]. Thus, the open forest line shifted 2.5–3.5 times less on moraine deposits, waterlogged soils and boulder fields on steep slopes than on sites without them on the test polygon (see Table 3 and Table 4) and throughout the studied subregion of the Polar Urals (see Table 6).

Our obtained research data demonstrated that in the studied subregion of the Polar Urals, the highest altitudinal position and the largest upward shift in the open forest line occurred on the south-oriented slopes if not restricted by edaphic constrains (Table 7). A similar situation was marked on the Khibiny massif (Kola Peninsula), where the open forests reached the highest altitude on slopes of southwestern exposure, which was 123 m higher than on the northeast-oriented slopes [64]. At this mountain massif, a maximum shift in the open forests (99–107 m) was revealed on the slopes of the southern exposition. On the Sukhie Gory massif (Putorana Plateau, northern part of Middle Siberia), the most significant altitudinal shift in the open forest boundary was discovered on the slopes of the southern (111 ± 74 m) exposures [65]. In southwestern Yukon (Canada), treelines shifted upslope higher on south- compared to north-oriented slopes [11]. In midlatitude forests of northwest British Columbia, shaded north-oriented slopes also showed lower productivity than those located on the sunnier and warmer south-oriented slopes [66]. First, variations in the altitudinal positions, responses and productivity of open forests could be explained due to differences in the incoming total solar energy, which, in cold climates, plays a significant role in the heat supply of different slopes. On the south-oriented slopes, the amount of incoming total solar energy was always greater as compared to north-oriented slopes. Thus, from June to August in the Yukspor weather station (Kola Peninsula, Khibiny massif, 910 m a.s.l.), 34.8 kcal cm⁻² of solar energy was received by the surface of the south-oriented slopes, with an inclination angle of 10°, whereas 32.9 and 30.6 kcal cm⁻² was received by west/east- and north-oriented slopes, respectively. If the angle of inclination of the surface were to increase, the difference between the northern and southern slopes would also grow up to 6.5, 8.8 and 17.7 kcal cm⁻² at 20°, 30° and 50° degrees of inclination, respectively [64]. Insolation causes the significant heating of crowns and trunks of trees, causing their temperature to become higher than the air temperature by 7–10 °C in the midday and, on average, by 0.8–2.9 °C per month on south-oriented slope [64]. In the eastern Alp treelines, the cambial activity of treeline trees lasted longer on the southern than on the northern slopes, due to the higher temperature of tree trunks at the first one [67]. Additionally, slope exposure can determine the start and end of snow melting in spring [68] and the growth period due to a greater warming up of tree trunks in the spring [67]. Differences in average altitudes in the upper forest boundary on the eastern and western slopes, which were illuminated by the sun in summer months and warmed up quiet equally, could be explained by the dominance of southwestern and western winds during the winter. This could cause a bigger snow pack accumulation on the northeastern and eastern slopes, and, as a result, could delay snow melting, thereby reducing the growing season in comparison with west-oriented slopes [65]. Thereby, it is important to keep in mind the slope exposure and its effects on tree growth and regeneration when considering models of tree stand dynamics on the upper limits of their distribution [69].

5. Conclusions

Our analysis of spatiotemporal changes of tree stand structures in the studied part of the Polar Urals showed that from 1960 to present day, total crown coverage significantly increased, and open forest lines moved up. Treeline advances originated in higher absorptions of solar radiation in newly forested areas and changes in C sequestration in plants and soils, therefore obtained by us knowledge about tree stand dynamics on the limit of their distribution could be widely applied for more precise forecasting processes of global warming. We suppose that a deeper understanding of the mechanics of tree stand dynamics based on a modelling approach that integrates climate, edaphic conditions (substrate type, thermal and moisture regime) and topographic features (slope exposition and inclination) with information on tree and tundra species ecological traits would be exceptionally helpful for improving predictions of tree stand responses to climate change in mountainous subpolar regions.