Trees in the Upper Treeline Ecotone in the Polar Urals: Centuries-Old Change and Spatial Patterns

Woody vegetation at the upper limit of its growth is a sensitive indicator of climate change. The aim of this study is to provide an analysis of the centuries-old spatiotemporal dynamics of larch trees at the upper limit of their growth (mountain massif Rai-Iz, Polar Urals, Russia). We used a ground-based method of mapping the remnants of trees that grew in the study area and died during the Little Ice Age. Aerial photographs from the 1960s and high-spatial-resolution satellite images from 2015 were used as data sources to define the locations of trees. Maps of the forest–tundra phytocoenochoras (areas of the terrain that are relatively homogeneous for one or more components of vegetation and/or other indicators) were created using a modified method of boundary detection between forest parcels with different stand densities. The proposed method of boundary detection between the main types of phytocoenochoras allowed us to identify a 15% total increase in areas of closed and open forest and areas with sparse tree growth, as well as a decrease in areas of tundra with single trees over these last decades. Using our spatiotemporal analysis of forest–tundra demographics over the last 50 years, we found that the number of trees in the ecotone had doubled. However, modern trees have not yet reached the areas occupied by trees in the past.


Introduction
Plants located at their upper and/or northern range limits are sensitive to even small changes in climatic conditions; therefore, many highland and northern territories are used as monitoring sites to assess the early response of biota to climate change (Shiyatov 1993;Kullman 2002Kullman , 2007Kullman , 2016Camarero and Guti errez 2004;Dufour-Tremblay et al 2012;Hellmann et al 2016;Pellizzari et al 2017;Chhetri 2018). Many studies indicate an upward shift of the treeline in the mountains of different regions of the world (Harsch et al 2009;Kullman 2012;Hagedorn et al 2014;Sanchez-Salguero et al 2018;Mamet et al 2019). At the same time, comparative analyses of changes in the boundaries of woody vegetation and the speed of these changes show that these vary according to region, for various reasons. For example, regions may differ in (1) species composition, age, and morphological structure of woody vegetation, as well as soil type and hydrological conditions; (2) heterogeneity of climate change and quality of woody vegetation data sources, which also differ on a spatial and temporal scale; and (3) methodologies used. A search for past boundaries of woody vegetation cover is not an easy task, especially for time intervals from half a century to several centuries ago, because there is often no information about where trees have previously grown due to rapid wood decomposition.
A unification of approaches and the development of methods for processing and analyzing data would reduce subjectivity and increase the degree of automation in determining the location of vegetation boundaries.
In most cases, the vegetation units revealed by a researcher are separated by transition zones called ecotones, as well as ecolines (Hufkens et al 2009). Ecolines are used to emphasize gradual changes between plant communities due to a gradient in at least one important ecological factor (Yarrow and Marín 2007). The appearance of sharp boundaries is a special case of the ecotone concept and, as a rule, is associated with contrasting environmental conditions (Hofgaard and Rees 2008). However, even studies that adhere to the concept of discrete vegetation cover accept that the borders distinguished are rarely sharp (Grossman et al 1994).
The characteristics of an ecotone are extremely dependent on the level of detail at which it is studied (Gosz 1993;Camarero and Guti errez 2000;Fortin et al 2000;Hufkens et al 2009). This is due to the properties of vegetation, which form a complex multiscale phenomenon. Considering this feature in environmental studies, including through abstraction, is crucial to understanding the patterns and processes associated with vegetation (Bateson 1979;Grossman et al 1994;Holtmeier and Broll 2017).
Publications summarizing the terms denoting the types of forest boundaries appear regularly (Hustich 1953;Gorchakovskiy and Shiyatov 1985;K€ orner 1998;Holtmeier 2003;Smith et al 2003;Chiu et al 2014;Bryn and Potthoff 2018). However, Holtmeier (2003) found that it was difficult to provide definitions that would correspond to all aspects of their use. Therefore, currently it is advisable to use the following approach: avoid trying to unify the use of the term ''forest boundary'' but attempt, in each case, to clarify what a scientist means when using this term. The improvement of methodological approaches to the collection, processing, and analysis of woody vegetation data, and obtaining new data and results, would help to bring positions closer together or resolve differences in determining boundaries between plant communities and assessing their spatial position.
Usually, researchers use the following characteristics to separate forest communities from nonforest communities: distance between trees (Shiyatov et al 2005), tree height (Kullman 1991), crown density and tree cover (Bakker et al 2008;Kharuk et al 2010), degree of tree root density (Sannikov et al 2012), trunk and crown shape (Szeicz and Macdonald 1995), forest islands (Bryn 2008), or various combinations of these criteria Wiegand et al 2006). The threshold values of these characteristics used by different authors may vary (Harsch and Bader 2011).
Differences in terminology and semantic content of terms lead to misunderstandings among researchers and can become an obstacle in comparing the results obtained by different authors. For example, this occurs when using changes in the spatial position of the treeline as an indicator of the response of woody vegetation to climate change. Such a comparison cannot be implemented for different areas if the publications do not specify what their authors understood as the treeline (Chiu et al 2014). At the same time, many researchers believe that the term ''treeline'' or ''tree line'' is intuitive and does not need to be explained. Jobbágy and Jackson (2000) compared treeline data presented in various publications and were forced to exclude from consideration results in which there was no clear definition of this term. This approach is appropriate, because the characteristics and their threshold values necessary for identifying boundaries differ, as does understanding of what the treeline is. The term ''treeline'' can refer to a transition zone (ecotone) extending from the boundary of a closed forest to the physiological limit of tree growth (Chiu et al 2014;Holtmeier and Broll 2017), one of the boundaries of this transition zone (Smith et al 2003), or a line of intermediate position between the upper limit of the closed forest and the upper limit of the growth of tree species. It can also refer to a line that connects the highest patches of forest at a similar position in relief, for example, for a series of slopes of the same exposure (K€ orner 1998).
The presence of a boundary definition is a necessary but insufficient condition to avoid ambiguity. Guided only by a definition, 2 researchers can draw boundaries in different ways. To unambiguously determine the location of a boundary, it is necessary to describe the method of identifying this line, a specific sequence of steps that allow any researcher to obtain reproducible results. The determination of the location of ecotones and their boundaries is always associated with the identification of spatially homogeneous areas of vegetation cover and the areas that separate them.
The purpose of this paper is twofold: (1) to provide an analysis of the spatial distribution of trees on the southeast macroslope of the Rai-Iz mountain massif in the Polar Urals in the second half of the 20th and early 21st centuries, as well as the remnant trees that died in the study area during the Little Ice Age, and (2) to give a quantitative assessment of changing areas of different types of phytocoenochoras over the last 50 years. The maps used to do so were created based on an earlier method developed for the automated determination of the boundaries between vegetation units. The method allows a researcher to increase the degree of objectivity in the allocation of the boundaries, and it is applicable for the detection of treelines.

Study area
The research area (66830 0 28 00 -66847 0 42 00 N, 65849 0 28 00 -65833 0 59 00 E) is located on the southeast macroslope of the Rai-Iz mountain massif (Polar Urals, Russia). The shaded relief of the study area and the adjacent territory is shown in Figure 1. The maximum elevation of the terrain in this area is 312 masl. The Rai-Iz mountain massif is composed of peridotite and gabbro bedrock. The surface of the region was formed under the influence of glaciers during the last global glaciation. Glacial deposits have formed a large number of ridges, mounds, and depressions. There are many streams and small lakes within this region. The eastern boundary of the study area lies on the left lateral moraine. The western boundary follows the course of the Engaiu River. The right bank of this river is formed by a right lateral moraine. Pure larch forests (Larix sibirica Ledeb.) dominate the woody vegetation within the upper limit of the ecotone (Shiyatov et al 2005). In the lower part of the ecotone, there are larch forests with single trees of Picea obovata Ledeb. and Betula tortuosa Ledeb.

Landscape photography
The camera points from which ground-based landscape photographs were taken in 1961-1962 and 2016 are shown in Figure 1. Repeated landscape photographs taken in these years are shown in Figure 2. For the analysis of landscape images, we used our previously developed method, which allows us to merge a photo and its sector of visibility (Fomin et al 2015;. The sector of visibility corresponds to the area of the terrain that is visible from the camera point in a given direction. The camera point is located at the top of the sector. The right and left sides of the sector, relative to the point from which the photograph was taken, correspond to the right and left sides of the photo. Each pair of repeated photographs presented in Figure 2 corresponds to sectors of visibility created based on a digital elevation model without considering the location of large objects or woody vegetation.

Remote-sensing and field data processing and analysis
Halftone aerial photographs from 1962 and 1964 for the eastern and western parts of the study region and satellite images from 2015 (Yandex Maps 2020) were used to identify the trees. These images were georeferenced using the QGIS geographic information system (QGIS 2020). Our studies defined coordinates of trees on the ground and measured their biometric parameters in 9 forest plots (Fomin et al 2019). This verified that the spatial resolution of aerial and satellite images reliably distinguished the same trees on both types of images (Supplemental material, Figure S1: https://doi. org/10.1659/MRD-JOURNAL-D-20-00002.1.S1). According to the analysis of the tree recognition results by 4 operators, trees that are more than 4 m tall could be identified. Based on the data obtained on the plots, the average area of tree crowns, according to the class of tree heights of 4-5 to 9-10 m, was within a range of 3.6 to 11.7 m 2 . When identifying trees in the images, both the size of the crown and the height of the tree (the shadow of which may be clearly distinguished) are important.
We created vector point layers, each point of which corresponds to the location of a tree at the beginning of the 1960s and in 2015. Landscape photographs were used to illustrate the effects of the appearance of woody vegetation in the forest-tundra and tundra and helped us to identify trees in some parts of the study area when we were analyzing the halftone aerial images.

Tree remnant mapping
The climatic conditions in the Polar Urals contributed to the preservation of fragments of trunks and tree roots that grew in the study area and died during medieval cooling (Shiyatov 2003;Mazepa 2005;Hagedorn et al 2014;Shiyatov and Mazepa 2015). The duration of this cooling period, also known as the Little Ice Age, was from the late 13th to the late 19th centuries (Briffa et al 1995;Shiyatov 2003;Oosthoek 2018). The remnants of large trees that died during this period can be seen in some of the landscape photographs ( Figure 2). Dendrochronological analysis of tree remnants and analysis of their location on several elevation transects passing through the ecotone (from the closed forest to the tundra with individual trees) in the study area between Yengaiu and Kerdomanshor rivers ( Figure 1) were performed by Shiyatov and Mazepa (Shiyatov 1993(Shiyatov , 2003Mazepa 2005;Shiyatov and Mazepa 2015).
We carried out fieldwork to determine the location of tree remnants in the upper part of the ecotone in 2016 and 2017 ( Figure 3). We defined the position of 7860 tree The shaded relief of the study area (dashed rectangle) and the adjacent territory, with camera points and their sectors of visibility. FIGURE 2 Repeated landscape photographs from camera points 221 (A and B), 397 (C and D), 408 (E and F), and 461 (G and H). Images in the left and right columns were taken by SG Shiyatov in the early 1960s (C and G in 1960, A in 1961, and E in 1962 and AP Mikhailovich in 2016, respectively. remnants. The location of each remnant was determined using a Etrex 10 global positioning system receiver (Garmin). The spatial error in positioning for this device is less than 3 m. The choice of the ground-based method of mapping tree remnants instead of remote-sensing methods, for example, using a small aircraft, arose from the need to search for fragments of dead trees in areas with dense cover of dwarf birch (Betula nana L.). The remnants were usually not visible with the use of aircraft drones in such areas. Trees (and their fragments) that died relatively recently can be distinguished from ancient tree remnants by the presence of bark fragments. Small tree fragments remaining on the stump also generally indicate that the tree died in the 20th century. We excluded such tree remnants from the analysis and presentation in Figure 3B.

Mapping boundaries in the treeline ecotone
Earlier, large-scale mapping of forest-tundra communities was carried out in the study area using our previously developed method (Shiyatov et al 2005(Shiyatov et al , 2007. This was based on identifying homogeneous areas (vegetation territorial units) called phytocoenochoras (Samoilenko et al 2009;Kholod 2016). A phytocoenochora is an area of the terrain that is relatively homogeneous for one or more components of vegetation and/or other indicators, such as the density of forest stands and/or the characteristics of forest site conditions. The average distance between the trees was used to assess tree stand density. Closed forests included areas where the average distance between trees was less than 10 m.
For open forests, the values of this parameter ranged from 10 to 30 m; for areas with sparse tree growth, they ranged from 30 to 60 m; and for tundra areas with single trees, they were 60 m and above (Shiyatov et al 2005).
In this paper, we define the term ''ecotone'' as the transition zone between types of phytocoenochoras and the term ''boundary'' as the line between them. We used the following threshold values to identify the previously mentioned types of phytocoenochoras, corresponding to the middle of each of the intervals described earlier: closed forest (average distance between trees up to 8.5 m inclusive, or stand density is at least 176 trees per hectare), open forest (from 8.5 to 25 m inclusive, or at least 20 but less than 176 trees per hectare), areas with sparse tree growth (from 25 to 55 m inclusive, or at least 4 but less than 20 trees per hectare), and tundra with single trees (more than 55 m, or less than 4 trees per hectare).
Maps of the forest-tundra phytocoenochora distribution in the early 1960s and in 2015 were created using a simplified modification of the previously developed method  in QGIS. Maps of the main stages of this method are shown in Supplemental material, Figure S2 (https://doi.org/10.1659/MRD-JOURNAL-D-20-00002.1.S1). In the first stage, a layer with Voronoi polygons was created using a point layer with tree coordinate data. In the second stage, cells that belonged to one type of phytocoenochora were selected according to their area values. In the next stage, the cells were merged into larger polygons by removing the internal borders between them, with subsequent filtering of polygons using their area values. Filtering was necessary to eliminate small polygons and thus avoid degenerate cases. For example, if such a polygon (combined from cells) belonged to a closed forest, its area was too small (contained too few trees) to be considered a closed forest.

Results
A map that characterizes the spatial distribution of trees in the ecotone of the upper limit of woody vegetation in the study region in the Polar Urals in the early 1960s and 2015 is shown in Figure 3. The collected data indicate that from the beginning of the 1960s until 2015, there was an almost twofold increase in the number of trees in the study area (from 14,377 to 28,344 trees). The location of remnants of trees that previously grew in the upper part of the ecotone and died during the Little Ice Age is shown in Figure 3B. The repeated landscape photos in Figure 2 allowed us to visually assess the changes that had occurred over roughly the last 50 years, for example, in the part of the study area that is visible in the photographs in Figure 2A and B, taken from camera point number 221. The area of the terrain shown in these photos can be seen on the map in Figure 1 (sector of visibility for point 221) and Figure 3 (direction of photography is shown by a red arrow). The photos and the map shown in Figure 3B indicate that trees have not yet reoccupied some habitats that they occupied in the past. This is confirmed by the many tree remnants we found in some upper parts of the ecotone in the study area.
According to the data obtained from sample plots 1-9, the crown cover for closed forest is within a range of more than 30%; for open forests, it is 10 to 30%; for light forests, it is 5 to 10%, and for tundra with individual trees, it is less than 5%. The maps of the main types of phytocoenochoras in the study area at the beginning of the 1960s and in 2015 are shown in Figure 4A and B. The areas of closed forest increased from 23 to 54 ha, those of open forests increased from 76 to 102 ha, and those with sparse tree growth increased from 56 to 116 ha. The area of tundra with single trees decreased from 621 to 505 ha. In terms of percentage (relative to the total area of the study region), the areas of closed forest, areas of open forest, and areas with sparse tree growth increased from 3 to 7%, from 10 to 13%, and from 7 to 15%, respectively, and the areas of tundra with single trees decreased from 80 to 65%.
During the period from 1964 to 2015, the rates of change in tree density (ie trees per hectare per year) within the boundaries of closed forests, open forests, light forests, and tundra with single trees of 2015 were 2.79, 0.85, 0.17, and 0.01, respectively. The ratio of tree density in 2015 to density in 1964 for the listed phytocoenochoras make up the following series: 1.8, 2.6, 4.2, and 7.8, respectively.

Discussion
A comparative analysis of the number of trees in the early 1960s and in 2015 suggests that regional climate warming, which is observed in the Polar Urals (Briffa et al 1995;Shalaumova et al 2010;Hagedorn et al 2014;Shiyatov and Mazepa 2015), was accompanied by significant changes in the ecotone of the upper boundary of woody vegetation.
According to data from the Salekhard weather station that is nearest to the research area, the average annual air temperature from the 1960s to 2010 increased by about 48C (Supplemental material, Figure S1: https://doi.org/10.1659/MRD-JOURNAL-D-20-00002.1.S1). This is shown by the appearance of a young generation of Siberian larch trees, both in the previously unforested or weakly wooded areas of the tundra and in areas with a higher density of trees. In open areas, young trees suffer greatly from snow and wind. The average age of trees with a height exceeding 4 m in nonforested and slightly forested areas (plots 5, 6, 8, and 9) is 25 years (there are no trees with a height of more than 4 m on plot 3). This means that they have appeared on sites since the mid-1980s, during the phase of increasing average annual air temperature (Supplemental material, Figure S1: https://doi.org/10.1659/MRD-JOURNAL-D-20-00002.1.S1). The location of modern trees and tree remnants allows us to conclude that in many parts of the study region, trees have not yet occupied the places they grew in in the past, and if they do, the number of trees in these areas is still significantly lower than in the past.
Dendrochronological analysis of the tree remnants in the altitude transect (partially across the eastern boundary in the upper-right corner of the study area, indicated in Figure 1 by the dashed line) revealed that some of these trees grew in the ecotone a maximum of 1300 years ago and that the maximum stand density on this transect was between the 11th and the 13th centuries AD (during the period of the Medieval Climatic Optimum) (Mazepa 2005). A similar trend was established for Swedish Scandes: the modern upper treeline position is below its location in AD 1000-1200 (Kullman 2015).
The reasons for the lack of trees in the areas where they previously grew may be as follows: lack of seeds at a distance of more than 40-60 m from seed trees; changing soil conditions (soil loss after the death of trees); competition from grass, low-shrub, and shrub-layer plants; and negative snow effects (snow corrasion of trees in the habitats, with a small amount of snow and shortening of the growing season in snowy habitats).
The presence of areas with groups of trees without woody remnants in the upper part of the ecotone, as seen from the map in Figure 3B, does not necessarily indicate that in these parts of the study region, modern trees reached areas inaccessible to trees in the past. Some larch trees can survive under a certain combination of microclimatic conditions in creeping form, returning to multistem or single stem forms after environmental conditions improve. This means that such areas may be refugia, in which Siberian larch has successfully survived an unfavorable centuries-long cold period since the 15th century (Devi et al 2008).
The proposed method of defining the boundaries between the main types of phytocoenochoras has a high degree of universality and objectivity. If there are quantitative criteria for identifying phytocoenochora types according to tree density, the procedure for determining the location of the boundaries between them due to the high level of formalization of the proposed method should be automated. This allows a researcher to reduce the subjective features of an operator on the final result.
The maps in Figure 4 show the changes in the configuration and areas of different types of phytocoenochoras. One of the features of the expansion of trees to the tundra or forest-tundra areas is treeless areas surrounded by trees, followed by their gradual covering by trees. This process is facilitated by the provision of seeds from trees surrounding such areas.
The uneven distribution of snow in the territory studied results from the transfer of snow by wind and the presence of obstacles in the form of convex relief elements, as well as areas with dense trees. As a rule, snow is blown away by the wind from the upper parts of the hills and accumulates in the lower concave areas or areas with trees. Once tree density is high enough, tree stands can accumulate snow. As a result, snow drifts can form in the inner part of the forested area relatively close to the boundary, forming an obstacle for the wind. The effects of snow accumulation on trees may be positive or negative (H€ attenschwiler and Smith 1999;Hagedorn et al 2014). Snow cover several meters in depth can delay the start of the growing season (ie shorten the already relatively short growing season). The reduction in the duration of the growing season has a negative effect on the survival and growth of young trees with heights lower than the snow cover depth. The locations with large accumulations of snow on the leeward sides of elongated landforms formed by glacial deposits are free from trees. This is because snow in these areas may be present during June or late July (ie the ground cover is freed from snow in the second half of the growing season). The remaining warm period of the year is not enough for the normal growth of Siberian larch trees (Shiyatov et al 2005). In addition, in the spring period, the subsidence of moist snow leads to the breaking of the lower branches of large larch trees (Shiyatov et al 2005). A positive effect of the accumulation of snow is that it can protect the younger generation of trees from frost and snow corrasion. On windblown upper parts of hills and ridges, where the snow depth is small, young larch trees suffer greatly from snow corrasion and survive mainly under the protection of large single trees or small groups of trees, as well as microrelief features, such as large stones or small ridges.
Analysis of geographical patterns of species range shifts under climate change may improve our ability to assess the persistence and shift of species (Breshears et al 2008;Maggini et al 2010;Lenoir and Svenning 2015). The basic types of species range shifts grouped into some broad categories (expand, march, lean, retract, and crash), described by Lenoir and Svenning (2015), closely relate to the problem of boundary detection. This or a similar classification

Conclusions
Woody vegetation located in the upper limits of its growing range is a sensitive indicator of climate changes. We studied, at the quantitative level, the features of the spatial position of Siberian larch trees in the treeline ecotone of the upper boundary of woody vegetation in the early 1960s and in 2015, as well as remnants of trees that died in the upper part of the ecotone in the 13th to 19th centuries. The number of trees in the ecotone has doubled over the last half century.
Maps were created of the spatial distribution of closed forest, open forest, areas with sparse tree growth, and tundra with single trees in the 1960s and 2015 using a method specially developed for this purpose . This method can be universally applied. It allows a researcher to determine boundaries between groups of points, in our case, based on the location of trees and quantitative thresholds of tree density.
The development of technologies for obtaining, processing, and analyzing high-spatial-resolution data that allow us to determine the location of individual trees in the upper treeline ecotone, in combination with our approach to determining the boundary between units of vegetation, opens the possibility of mapping the upper-boundary forests in the mountains of other regions of the world. Our study contributes to understanding the particularity of climatedriven changes in the spatial position of woody vegetation at the upper treeline. In a broader sense, our studies can be used to solve problems of sustainable forest management, in particular improving existing methods and developing new ones for quantifying and predictive modeling the dynamics of forest ecosystems using surface information and remotesensing techniques.

A C K N O W L E D G M E N T S
This work was supported by the Russian Foundation for Basic Research (project 18-34-00803 mol_a). Some studies related to the assessment of the spatial distribution of snow at the beginning of the growing season were carried out within a grant from the Russian Science Foundation (project 17-14-01112) and a grant from the Russian Ministry for Education and Science (FEUG-2020-0013). We express our gratitude to Dmitriy Golikov and Yegor Agapitov for their help in setting up forest plots and data processing, and also to Denis Kapralov for help in mapping the tree remnants.

Supplemental material
FIGURE S1 Location of the forest plots in the study area, with an aerial photograph from 1962, a satellite image from 2015, and time series of mean annual temperature and annual precipitation.

FIGURE S2
The main stages of the algorithm for defining the boundaries between different types of phytocoenochoras.