Spatiotemporal dynamics of encroaching tall vegetation in timberline ecotone of the Polar Urals Region, Russia

The study region of the Rai-Iz massif is located on the eastern slope of the Polar Urals range (200–350 m asl.), in the Sob' River watershed (ca. 66°46' N, 65°49' E) underlain by the continuous permafrost (figures 1(a) and (c)). The bedrock is mainly comprised of Paleozoic amphibolite and crystal granodiorite (Mazepa et al 2011). The most common soil in the Polar Urals is mountain-tundra gleysol, and rare inclusion of mountain-tundra turf soils. On the ultramafic massif of Rai-Iz, soils are characterized by neutral and near-neutral pH (Kataeva et al 2004). According to data collected at Salekhard meteorological station (55 km southeast of the study area, 35 m asl.) between 1883 and 2005., the mean annual, mean monthly minimum (January) and maximum (July) air temperatures are −6.7°, −24.4°, and 13.8 °C, respectively. The mean annual precipitation is 600 mm with 50% as snow and sleet. The growing season lasts from mid-June to earlyAugust with a mean frost-free period of 94 d. Westerly wind prevails, while northeasterly wind is also present in early summer (see S1 available online at stacks.iop.org/ERL/17/014017/mmedia), with an average wind speed of 1–3 m s⁻¹ and maximum speed observed during the winter season (Mazepa 2005).

Zoom In Zoom Out Reset image size

Starting in the early 1960s, six permanent altitudinal transects of 300–1100 m long and 20–80 m wide were established in the eastern foothills of the mountain range for long-term monitoring of spatiotemporal dynamics of alpine forest-tundra and forest-meadow communities. Transects typically start at the upper timberline boundary, where live trees are not yet present, and extend down into the upper limit of closed forest (Mazepa et al 2011). This study focuses on the earliest ecotone transect within which complete data of tree characteristics, terrain elevation, and snow distribution have been collected for multiple years (figure 1(d)). This transect also exhibits varying topography, as compared to the other transects. The transect is 860 m long and 40–80 m wide (the total area is 5.6 ha, with the upper left corner at 66°48'57'' N, 65°34'09'' E) and oriented along the direction of the predominant westerly wind (Mazepa 2005, Shiyatov and Mazepa 2011).

This transect is located close to the southern boundary of the bioclimate subzone E (Walker et al 2005) as shown in figure 1(b). Tree stand composition within the timberline ecotone is relatively simple: larch (Larix sibirica Ledeb.) forest-tundra communities dominate in the upper part, while open larch forest with Siberian spruce (Picea obovata Ledeb.) and downy birch (Betula tortuosa Ledeb.) occur in the lower part of the area (Shiyatov et al 2005). Dwarf birch (Betula nana L.) is abundant in the lower and middle parts of the transect. Names of plant species correspond to the Integrated Taxonomic Information System (ITIS 2021). Historically, boreal trees grow at the fringe of the forest-tundra ecotones in various growth forms (creeping, prostrate, single-stemmed, and multi-stemmed), which reflect adaptations to environmental changes (Mazepa and Devi 2007, Devi et al 2008). Over 90% of the young trees appearing in the study transect after 1950 are single-stemmed (Mazepa et al 2011). The studied forest stands, and recent vegetation growth have emerged without significant impacts from reindeer grazing, fire, and human activities (such as logging) over at least the last millennium (Mazepa 2005, Devi et al 2020).

The records of three census campaigns in 1960, 1999, and 2011 are used in this study to quantitatively assess the changes in the composition, structure, and spatial distribution of the forest-tundra communities. Specifically, detailed mappings of all alive and dead tree locations, measurements of their allometric and geometric characteristics such as height, crown size, branch position, and diameter at multiple heights were produced during these campaigns. All three censuses predominantly recorded larch with the occasional presence of spruce along the study transect. The numbers of trees for the 1960, 1999, and 2011 censuses are 1853, 1700, and 1398, respectively, after pre-screening (see supplementary materials (S2)).

After reconstructing the missing tree diameter data using an allometric function (see S2), we compute the change of tree diameters ΔD between censuses as the proxies for the increase of biomass and tree productivity, and then relate them to environmental features. Specifically, to compare changes of tree diameters ΔD for different periods (i.e. from 1960 to 1999, and from 1999 to 2011), the relative diameter change ΔD_rel = ΔD/D₀, where D₀ is the tree diameter at the first year of the corresponding period, is computed first. ΔD_rel is then normalized by the period of time T (years) to get the annual relative change rate of stem diameter ${\dot D_{{\text{rel}}}} = \sqrt[T]{{1 + \Delta {D_{{\text{rel}}}}}} - 1$ assuming a constant annual change rate for a given period. This rate was computed only for trees recorded both at the start and end year of the time interval between the two censuses. The changes in tree diameter may also be used as a proxy of tree biomass change through its relationship with sapwood (see S2).

Snow data were collected along the study transect in 1961, 2006, and 2013 in the spring, when snow depth (Z_S) was at its peak by direct measurement, tree marking, and visual assessment (see S3). The relative snow depth, defined as ${\hat Z_{\text{S}}} = \left( {{Z_{\text{S}}} - {Z_{{\text{Smin}}}}} \right)/\left( {{Z_{{\text{Smax}}}} - {Z_{{\text{Smin}}}}} \right)$ where ${Z_{{\text{Smax}}}}$ and ${Z_{{\text{Smin}}}}$ are the maximum and minimum snow depths within the transect for a given survey year from the three surveys, is shown in figure 2.

Zoom In Zoom Out Reset image size

Digital terrain models (DTMs) are needed to capture topography at fine scales, but there are limited areas in the Arctic for which high-resolution DTMs are available. There is a range of digital elevation model (DEM) products that can provide high-resolution elevation maps. However, these DEM products are affected by the presence of surface biomass masking the location of actual topography. To reconstruct the DTM for the transect, we combined a resampled 2 × 2 m resolution grid from a 20 × 20 m regularly spaced field terrain measurements of the transect, and a 2 m resolution DEM product called ArcticDEM (Porter et al 2018). After co-registering two ground sources (i.e. field measurement and ArcticDEM), we manually selected bare ground grids from ArcticDEM together with the field measurement to extract a high-resolution DTM for the study transect by natural neighbor interpolation (Tily and Brace 2006).

Terrain curvature is one of the fundamental topographic properties and an important driver of hydrologic processes (Moore et al 1991). Two types of curvature indices, i.e. profile curvature and planform curvature defined by the Environmental Systems Research Institute (ESRI 2021), are used to represent micro-topographic features in this study. Specifically, profile curvature is the curvature of the surface in the direction of the steepest slope, while planform curvature is the curvature in the plane perpendicular to the direction of the steepest slope. For profile curvature index, a negative value indicates an upwardly convex or dome-shaped surface, while a positive value represents a concave surface or depression; for planform curvature index, negative and positive values indicate a laterally concave and convex surface, respectively. Note that the opposite sign is used in the definition of planform curvature concavity/convexity to that of profile curvature. In conventional terrain analysis, the curvature indices are calculated from DTM as the second-order derivative of elevation map using various geomorphometric methods (e.g. Zevenbergen and Thorne 1987). The curvature indices calculated from the resampled DTM with various resolutions are used to capture topographical features at different scales (see S4).

To attribute topographic locations to niches of higher or lower favorability for tree encroachment, we computed the weighted average curvature index using stem diameter at locations where trees are present as a weight:

$$\begin{equation}{\bar \kappa _p} = \frac{{\sum {D_i}{\kappa _{i,p}}}}{{\sum {D_i}}},\end{equation} \tag{ 1 }$$

where ${D_i}$ is the diameter of tree i, ${\kappa _{i,p}}$ is the curvature index of a terrain grid cell where the tree is located, and the summation is for all trees along the transect. The index p represents either planform or profile curvature index. Noted that ${\kappa _{i,p}}$ as well as the number of trees within a given grid cell change with resolutions, so the summation is not carried out over the same set of trees for different grid resolutions. The above formulation emphasizes the contributions of curvature indices of terrain locations with higher D. The rationale is that stem diameter is argued to be an indicator of biomass and long-term productivity (section 2.2). If there is an enhanced tree growth at topographic locations of specific curvature, it can be elucidated by comparing ${\bar \kappa _p}$ in (1) with the average curvature of transect topography as the arithmetic mean of curvature over all grid cells along the transect at different grid resolutions.

Spatial point pattern analysis is commonly used to characterize the patterns of the spatial distribution of tree species (e.g. Condit et al 2000). Most of these analyses use a family of statistics derived from the Ripley's K-function (Ripley 1976) to quantify the species clumping characteristics by comparing them with complete spatial randomness (CSR). In this study, following the approach proposed by Plotkin et al (2002), we apply a density-based spatial clustering method with the presence of noise (i.e. points located too far from other points) (Ester et al 1996) to capture spatial pattern of tree distribution. The advantage of this method is that it does not require a priori number of clusters, such as in k-means analysis (Macqueen 1967), and can identify clusters of arbitrary shape. Clusters of tree locations (referred to as 'points') are derived at different spatial resolutions using the records of the three tree censuses (section 2.2). The algorithm requires two parameters: neighborhood search radius (d) and the minimum number (q) of points required to form a cluster located within the search radius. By setting the parameter q to 1 in all analyses, we make the clustering process depend only on the parameter d, which thus represents the maximum distance between any pair of connected trees belonging to the same cluster. Therefore, for any given value of d all stem locations are partitioned into a unique set of clusters.

Tree clusters were formed at different spatial scales by changing the parameter d from 0 to the maximum distance between locations of any two connected trees in the transect (65 m). For a given d, the total number of obtained clusters is denoted as m and the size of an ith cluster as c_i . The sample size, i.e. the total number of trees in the transect, is denoted as $n = \mathop \sum \limits_{i = 1}^m {c_i}$. The unique arrangement of clusters for a specific value of d can be summarized with the normalized mean cluster size:

$$\begin{equation}\hat c = \,\frac{1}{{{n^2}}}\mathop \sum \limits_{i = 1}^m {c_i}^2.\end{equation} \tag{ 2 }$$

This variable represents the probability that two randomly chosen tree locations belong to the same cluster (Plotkin et al 2002). For example, $\hat c = 1/n$ when $d = 0$, corresponds to a clustering pattern that each point forms its own cluster; $\hat c = 1$ when d equals the maximum pairwise distance between any two connected tree locations, i.e. all trees form a single cluster. The clustering pattern varies for intermediate values of d and can be obtained by recording the concomitant change of $\hat c$, known as the mean cluster size curve $\hat c\left( d \right)$. To make the results for censuses with different sample sizes comparable in the mean cluster size curves, we use the normalized distance $\hat d$ defined as

$$\begin{equation}\hat d = d\sqrt {n/A} ,\end{equation} \tag{ 3 }$$

where A is the transect area. Spatial patterns of clustering can then be compared with CSR (see S5).

The sampling distribution of the annual relative change rate of stem diameter ${\dot D_{{\text{rel}}}}$ for the two periods between the three censuses, (1960–1999 and 1999–2011) is shown in figure 3(a). The mode of the distribution increased from 2%/year to 3%/year, indicating accelerated growth in the later period. The right-skewed distribution of the second period implies that the proportion of trees with higher growth rates is larger in the recent period. Figure 3(b) shows the relationship between ${\dot D_{{\text{rel}}}}$ and D₀ for the two periods between the three field surveys. The maxima of ${\dot D_{{\text{rel}}}}$ form a clear upper boundary with respect to D₀ for both periods, indicating the upper limit of the annual growth rate and the results illustrate an increase of these growth rates from the earlier to the later period (red and blue lines). The highest growth rate decreasing with D suggests that larger, older trees grow slower than smaller, younger trees. Growth rates can vary significantly, between just above zero and to a few percent of their stem diameter per year.

Zoom In Zoom Out Reset image size

Figure 4 compares the average, D-weighted terrain curvature index for grid cells with live trees (equation (1)) and the transect average index (i.e. uses all grid cells), as a function of DTM grid cell resolution. The results indicate scale-dependent differences with respect to the transect average for both planform and profile curvatures, and for all censuses. The difference is significant for DTM grid sizes ranging from 2 to 14 m, with the largest difference occurring at the grid size of 10 m. This implies the preference of tree growth in terrain curvature is characterized most significantly at this scale (i.e. 25–30 m length scale of landforms, which is the size of the moving window used to compute the curvatures, see S4). At locations with live trees, the profile curvature index is more negative, and the planform curvature index is more positive than the overall transect averages, implying that trees tend to cluster at sites with convex topography. The curvatures of the locations with trees and the transect average become indistinguishable for grid cells of sizes greater than 14 m. Since the number of grids along the transect decreases with grid size, this leads to the convergence of the curvature indices to the same value.

Zoom In Zoom Out Reset image size

The resultant cluster size curves for the three censuses are shown in figure 5. The sharp transitions in the three curves at various distances $\hat d$ (equation (3)) correspond to the 'continuum percolation' phenomenon (Meester and Roy 1996, Plotkin et al 2002), i.e. the aggregation of isolated smaller clusters into a larger cluster. The corresponding distance is known as the 'critical distance' or 'percolation threshold'. When $\hat d$ is below the threshold, elements in the system of interest are disconnected and form many smaller-sized clusters; above the percolation threshold, the elements are aggregated to form larger clusters. In theory, for true CSR (i.e. if the sample size is infinitely large) the cluster size curve will show a discontinuous transition corresponding to the step function at a critical distance of $\hat d \approx 1.2$ (an abrupt transition for the developed CSR curve in figure 5 is not discontinuous due to the finite sample size and transect edge effects on the computation of $\hat d$).

Zoom In Zoom Out Reset image size

The cluster size curves constructed using the actual tree locations from all censuses have three transitions at varying critical distance $\hat d$. Their differences with respect to the CSR curve imply a non-random spatial tree distribution at different spatial scales. The plateaus of the curves between the abrupt transition in cluster sizes represent the ranges of d within which the change of cluster pattern is insignificant. The number of plateaus (two for this case) indicates the number of scales of non-random aggregation (Plotkin et al 2002). Therefore, one may perform a clustering analysis using d corresponding to a plateau region to represent such a stable distribution at a selected scale. Here we choose the values of normalized distance corresponding to the start of plateaus immediately after the occurrence of 'continuum percolation'. The following clustering analysis uses $\hat d$ corresponding to the first plateaus (${\hat d_1}$) in figure 5: 1.24, 1.4, and 1.51 (i.e. the actual distances are d₁ = 6.7, 7.9, and 9.4 m) for the 1960, 1999, and 2011 censuses, respectively. For the second plateau (${\hat d_2}$) in figure 5, the normalized distances are ${\hat d_2}$= 2.37, 3.63, and 5.58 (d₂ = 12.8, 20.5, and 34.7 m). When a cluster pattern corresponding to ${\hat d_1}$ is defined as the 'primary cluster set', then the value of ${\hat d_1}$ is a representative proxy of the distance among trees within any cluster of this set, while ${\hat d_2}$ represents a proxy of the distance between borders of clusters formed by the primary set. As figure 5 shows, over time relatively smaller changes of the distribution pattern occurred within the primary cluster set, as reflected by the comparable magnitudes of ${\hat d_1}$ across the three censuses. The distances between the cluster borders however become larger, as illustrated by the increase of ${\hat d_2}$, particularly for the census of 2011 (mainly due to the mortality of trees at cluster borders).

Figure 6 shows the spatial distributions of tree clustering for the two critical distances $\hat d$. Transitions of cluster patterns can be detected by comparing plots in the same row. For example, when the distance changes from d₁ = 6.7 m to d₂ = 12.8 m (the census of 1960), the number of clusters decreases significantly from 102 to 21 due to cluster aggregation. The same 'continuum percolation' phenomenon is observed for the other census data of 1999 and 2011. Comparing plots in the same column reveals that the number of primary clusters decreases over time. This phenomenon can be explained by the disappearance of cluster 'noise', i.e. those trees that form their own cluster. Notice that for the primary cluster set (the left panel of figure 6), the sum of cluster sizes of the three largest clusters is over 80% of the total number of surveyed trees (not shown).

Zoom In Zoom Out Reset image size

Figure 7 illustrates the temporal change of the distribution of the largest cluster from 1960 to 2011. In 1960, there were many small trees scattered in this cluster and larger trees were mainly located in the middle part. A large fraction of small trees previously located in the lower part disappeared in later years, while those in the upper part grew bigger and likely seeded new trees that appeared in the transect in 2011.

Zoom In Zoom Out Reset image size

Concurrently with the change of the tree distribution pattern, the area with high snow depth had shifted from the lower part to the middle part of this cluster. Similar transitions of areas with high snow depth are also apparent in the regions outside of this cluster (not shown). The area of high snow accumulation shifted from the lower part to the middle-upper part of the transect, coinciding with the growth and densification of trees in clusters at the top.