Introduction

There has been an increase in forest cover in the mountain areas of many developed countries that was preceded by a period of deforestation caused by the expansion of human activity into the higher parts of those mountains (Hall et al. 2002; Meyfroidt and Lambin 2011; Price et al. 2017). This phenomenon has been recognised in many mountain areas globally and defined as the ‘forest transition’ (Mather 1992; Rudel et al. 2005). In Europe, a dominant trend of abandonment and reforestation has been found in mountain areas since 1950s (Jepsen et al. 2015). The peripheral location of most European mountain regions facilitates economic decline, followed by depopulation, a decrease in traditional farming activities, land abandonment and subsequent reforestation on abandoned agricultural land (MacDonald et al. 2000; Plieninger et al. 2016).

The quantification of reforestation on abandoned land is of high importance due to various environmental and socio-economic consequences (Kolecka et al. 2017). Contemporary and future forest expansion on abandoned land induces a significant landscape change and causes several consequences that may be perceived as either negative or positive. The positive effect concerns an increase in carbon sequestration (Kuemmerle et al. 2011) and an extension of new habitats for forest species (Ceausu et al. 2015). On the other hand, reforestation causes a decrease in biodiversity (Bennett and Saunders, 2010) and a loss of cultural landscapes (Fischer et al. 2012).

Research on changes in mountain forest cover is prevalent and most often focuses on the quantitative analysis of this phenomenon or the driving forces (Tasser et al. 2007; Rutherford et al. 2008; Munteanu et al. 2014; Kerckhof et al. 2016; Kozak et al. 2017). The results of these studies, conducted most often in large-scale areas, even covering entire mountain arcs, indicate that the places in the least convenient locations are the first to be abandoned (Lasanta et al. 2017). Little scientific research has been performed concerning reforestation factors on a local scale in mountain areas where there has been no decline in population (Ostafin 2009). What is more, there is a lack of studies based on consecutive series of aerial photographs that would make it possible to distinguish the factors determining the rate of forest succession over time from the features of plots characterised by large forest cover.

In research on the driving forces of reforestation, biophysical factors and socio-economic evolution are integrated to increase the understanding of landscape dynamics (Mottet et al. 2006). When considering changes concerning the entire mountain arc (regional scale), the influence of biophysical factors on reforestation resulting from land abandonment is visible (high slope inclination hindering mechanisation of work, poor soil quality affecting low yields, considerable remoteness and difficult accessibility for transportation) (Tasser and Tappeiner 2002; Kozak 2010; Müller et al. 2013). These factors directly determine economic viability (Flury et al. 2013). In turn, studies conducted on a small scale indicate that small enclaves of non-forest areas remain for a long time in the higher parts of the mountains, despite their unfavourable location (Sobala 2018).

This article aims to determine the role of different factors in the reforestation of abandoned plots on particular mid-forest glades in the higher parts of mountains and to detect the feature of these plots that are covered prevalently with forest. The course of reforestation of abandoned plots results from the impact of the factors which determine land abandonment and the factors which decide the course of secondary forest succession (Kolecka et al. 2017). The detailed study was conducted in the Western Carpathians where trends of forest cover changes follow the path of ‘forest transition’ and are typical of other mountain areas in developed countries (Mather 1992). Because reforestation trends vary in different parts of the Carpathians (Bucała and Starkel 2013; Kolecka and Kozak 2019), more detailed analyses are needed to distinguish the character and drivers of reforestation not only on a general but also on a local scale.

Materials and methods

Study area

The Carpathians are a prominent mountain arc in Europe. Their western part stretches along the border of Poland with Slovakia from the Low Beskids range of the Eastern Carpathians toward the Moravian region of the Czech Republic and the Austrian Weinviertel. The area of the Western Carpathians comprises about 70,000 km2 (Kondracki 1989). The highest elevation is the Gerlachovský štít (2655 m asl), located in the Tatra Mountains. However, besides the region of the Tatra Mountains, this mountain belt covers areas of mid- and low-mountains, foothills and valley bottoms (Kozak et al. 2013).

The maximal range of deforestation in the Western Carpathians occurred in the nineteenth century. A significant increase in forest cover started in the second half of the twentieth century and continues today (Kozak 2010). However, the afforestation of mountain pastures and glades located in the highest part of the Western Carpathians started in the second half of the nineteenth century (Sobala 2018).

The study was carried out on three levels: the local level, case study level and plot level. On the local level, two study areas were investigated to reconstruct changes in forest cover since the mid-nineteenth century: the part of the Barania Góra range (the Silesian Beskids) and of the Wielka Racza range (the Żywiec-Kysuce Beskids) (Fig. 1). They each occupy 45 km2, and their absolute altitude ranges from about 500 to over 1200 m asl, which is typical of the Western Carpathians. These two areas were chosen as their location in the higher parts of the Western Carpathians means that they were abandoned first (Sobala et al. 2017). The study area is situated within the lower subalpine forest and upper subalpine forest. The plant communities typical of these layers are Luzulo luzuloidis-Fagetum, Dentario glandulosae-Fagetum, Abieti-Picetum (montanum) and Plagiothecio-Piceetum. Nevertheless, the study area is mainly covered by artificial spruce monocultures characterised by a very limited species composition (Sobala and Rahmonov 2020). Non-forest patches with dominant semi-natural communities often occur within the forests. These patches mainly include various types of meadows and pastures, characterised by the following plant communities: Gladiolo-Agrostietum, Hieracio-Nardetum, Arrhenatheretum medioeuropaeum, Cirsietum rivularis and Juncetum effuse (Sobala et al. 2017).

Fig. 1
figure 1

Study area. Local level, black line; case study level: A1, Radziechowska; A2, Ostre; B1, Praszywka

On the case study level, three test areas (glades) were chosen to reconstruct a path of reforestation since 1977: the Radziechowska (A1) and Ostre (A2) in the Barania Góra range and the Praszywka (B1) in the Wielka Racza range. The plot level was used to assess the determinants of reforestation. It was connected with the fact that land fragmentation typical of the Western Carpathians influences spatial dispersion of reforestation (Kolecka et al. 2015). As the individual farmers made land-use decisions, the pattern of reforestation now differs from one plot to the other (Tasser et al. 2007).

Materials

In this study, five different types of cartographical materials were used (Table 1): archival maps and panchromatic (B/W) aerial photographs, archival and contemporary natural colour (RGB) orthophotomaps, a digital terrain model and the master map.

Table 1 Cartographical materials used in the study

The digital terrain model is a product based on airborne laser scanning (ALS) with a density of 4-8 points/m2. It was made in a flat rectangular coordinate system 1992 with a 1 × 1 m resolution and an average height error in the range of 0.2 m. LiDAR point clouds were acquired within the national ISOK (Informatyczny System Osłony Kraju) project. The master map is a basic geodetic and cartographic database covering the entire area of Poland that contains cadastral plots.

Methods

Different methods were applied on the three levels. On the local level, forest cover changes in 1848-2015 were investigated, while on the case study level, forest cover changes since 1977 were analysed; on the plot level, the Reforestation Index was determined and the determinants of its changes were analysed.

Firstly, all the historical maps and aerial photographs were georeferenced. In each case, georeferencing was specifically adjusted to the quality and type of data so as to achieve the best possible results for each series (Sobala 2012; Affek 2013). The photographs were orthorectified using spline transformation and ground control points (GCPPs) measured on orthophotomaps from 2019 or the master map. In each photograph, 28-193 GCPs allowed the result to be achieved. The spline transformation is a true rubber sheeting method which optimizes for local accuracy (https://pro.arcgis.com/en/pro-app/latest/help/data/imagery/overview-of-georeferencing.htm).

Secondly, vector maps of forest cover were created in ArcGIS® ver. 10.8. A topology construction tool was used to detect and eliminate errors (Maras et al. 2010). On the basis of archival maps, forest and non-forest areas were clearly distinguished, whereas on the basis of aerial photographs, forest, non-forest areas and groups of fairy thorns were distinguished. In the end, all cartographical materials were transformed into European Terrestrial Reference System 1989 (ETRS89)/Poland CS2000 zone 6 (EPSG:2177).

Subsequently, for each plot and for consecutive time sections, the Reforestation Index (RI) was calculated:

$$RI=\sum_{i=1}^np\ast r,$$

where: RI, Reforestation Index; p, share of i-th form of land cover within plot; r, values assigned to land cover form, n, number of land cover forms.

This index uses fuzzy logic, taking into account the following values assigned to land cover forms (r): 0 for non-forest areas, 1 for forest areas and 0.5 for groups of fairy thorns. The higher the value of the indicator, the more advanced the reforestation within particular plots. The data consisted of 389 plots, with a total of 2415 observations.

On the plot level, the data was analysed with R 4.02. statistical package. In the initial step of the investigation, exploratory analyses were conducted. First, distributions of the dependent variable and the predictors were assessed, and some variables were transformed to correct for high skewness. Second, the rate of change in the RI and the probability of reaching RI = 1 for each plot in a Bayesian multilevel zero-one inflated beta regression were estimated. Third, plot-specific model coefficients (e.g. RI rate of change) were regressed onto the set of selected predictors using multivariate regression. Significant predictors with standardised regression weights β > 0.2 were used in the main analyses (see Online Resource 1 for a list of all considered predictors).

RI was modelled with zero-one-inflated beta regression, which allows not only the conditional change of the mean as a function of predictors to be modelled, but also the conditional probability of zero and one occurrences (Nishii and Tanaka 2013). For each specific case study, 6 or 7 observations were available (depending on the glade), one for each year of observation. To control for the variation between glades in the average RI and the rate of RI change over time, random effects for the intercept and the year coefficients for the mean and zero-one inflation parameters were included in the model. The final model, a multilevel zero-one inflated beta regression model, was fitted to the data in the Bayesian framework with a brms package via MCMC sampling (Bürkner 2017, 2018).

Results

Basic characteristics

There has been an increase in forest cover in both studied mountain ranges since 1848 (Fig. 2a). The annual year increase was 3.2 ha/year in the Barania Góra range and 4.6 ha/year in the Wielka Racza range. The most significant surge of growth was noted in 1933-1960 in the Barania Góra Range, whereas in the Wielka Racza range, reforestation has accelerated since 1933.

Fig. 2
figure 2

Reforestation stages: a local level, b-d plot level: Radziechowska, Ostre, and Praszywka, respectively

Between 1977 and 2019, the forest area increased in the Radziechowska, Ostre and Praszywka glades from 49 to 66%, 45 to 73% and 34 to 82%, respectively (Fig. 2b-d). The largest surge of growth was noted between 1983 and 1993 in the Ostre and Praszywka and between 1998 and 2003 in the Ostre and Radziechowska.

The distribution of the RI clearly changed over time both in terms of the mean of the RI proportion and in terms of the decrease in the number of plots with RI = 0 and increase in the number of plots with RI = 1 (Fig. 3).

Fig. 3
figure 3

Distribution of RI proportion across years of observation. Black bars show histograms calculated excluding zeros and ones. Blue vertical lines show numbers of plots without trees (RI = 0) and with full forest cover (RI = 1). Red vertical line shows mean calculated with exclusion of zeros and ones. Light blue vertical line shows 0.5 on the x axis for reference

The results of the initial phase of the analyses indicate that the predictors mean solar, SD absolute height, area, mean exposure (coded N, W, E, S), SD exposure and mean slope may be related to the RI in general and/or to the change in the RI over the years. Descriptive statistics for these variables are presented in Online Resource 2. Due to the high skewness area, SD exposure and SD absolute height were log-transformed for further analyses.

The selected predictors were mostly unrelated to one another, with the exclusion of SD absolute height, which was strongly related to the area and weakly related to the mean slope (Online Resource 3). What is obvious, the highest levels of mean solar radiation were observed for lands with south exposure, and the lowest levels of the predictors were observed for north exposure.

Additional predictors included in the final model were test area (glade), RI1977, and planned forestation. Fourteen (4%) plots were located in A1, 81 (21%) plots in A2 and 294 (76%) plots in B1. Due to the small number of plots in A1 and the fact that the relationships between variables in this case were similar to the ones observed in B1, in the final analysis, these two glades were treated as one category ‘A1+B1’. The RI1977 represented the initial state of the plots, with ‘0’ indicating that the land had no trees (RI = 0; 12% of the plots) in the first year of observation, and ‘>0’ indicating that the RI in the first year was higher than zero (88% of the plots). The planned forestation was indexed whether the land was intentionally forested by humans (21 cases, 5.4%).

RI change modelling

The continuous predictors were scaled, and categorical predictors were coded with orthogonal contrasts. The year of observation was treated as a continuous predictor. Apart from making the estimation procedure more efficient, these transformations also resulted in comparable regression weights.

Modelling of the RI change was conducted iteratively, where in each iteration, different structures of interactions between year, glade and continuous predictors were tested. The model presented below represents the most parsimonious model that is a model which includes only interactions that were consistently statistically credible across the tested models. The Bayesian R2 value (Gelman et al. 2018) of the fixed effects for the reported model (i.e. variance explained by the population-level effects of predictors) was almost the same (Table 2) as the value of the statistic for a model which included all two-way interactions, R2 = 0.47, 95% CI: [0.45, 0.49].

Table 2 Summary of the posterior distributions of the fixed-effects of the Bayesian multilevel zero-one inflated beta regression with RI as dependent variable

Posterior distributions of the fixed-effects of the Bayesian multilevel zero-one inflated beta regression are summarised in Table 2. We observed that the average RI and proportion of plots with RI = 1 increased credibly over the years. Plots with higher mean slope, south exposure and which started out with some trees had higher average RI. Glade A2 had noticeably higher average RI than A1+B1 (Fig. 4).

Fig. 4
figure 4

Predicted average RI as a function of predictors included in the final model. a-b: blue line is predicted mean RI and the shade is 95% credible interval of the prediction. c-f: black point is predicted mean RI, error bar shows 95% credible interval of the mean. g-i: solid line is predicted mean RI and corresponding shade is 95% credible interval of the prediction

Importantly, exposure affected how fast the plots reached total forest cover, with south exposure characterised by the highest rate. This was caused by higher average RI in plots with south exposure in the first year of observation. Additionally, plots which started out with some trees reached full forest cover faster. The same tendency was observed within plots with average mean slope (15°).

Importantly, four factors in the model credibly affected the rate of change in RI over the years. Average RI increased credibly faster in A1+B1 than in A2 (Fig. 4). The change in the RI over the years was also credibly affected by the mean solar radiation, area and SD exposure (Fig. 5). The effect of mean solar radiation was the strongest: in the plots with the lowest levels of mean solar radiation, RI change was noticeably faster. Furthermore, the change in RI over the years was fastest among plots with the smallest area and low level of SD exposure. The last effect was weak but credible.

Fig. 5
figure 5

Predicted change in RI over the years as a function of mean solar radiation, mean slope, area and SD exposure

Discussion

In this paper, reforestation was analysed at a local scale, which made the detection of subtle changes possible. Small patches of forest succession, such as groups of fairy thorns, were included. In the study area located in the highest part of the mountain range, reforestation started in the second half of the nineteenth century. Nevertheless, the most significant surge of growth in the study area was noted after 1930. That was about 20 years earlier compared to the entire Polish Carpathians (Kolecka et al. 2017). The current proportion of forest cover is the consequence of the rate of secondary forest succession, the location of the plots abandoned by the owners and the time elapsed since the owner’s decision to stop cultivation. Kolecka et al. (2017) claim that spatial determinants of farmland abandonment may exhibit significant similarities to spatial determinants of increasing forest cover. This paper shows a difference between features of plots with a faster rate of reforestation and plots covered entirely with forest.

The fastest changes in the Reforestation Index (RI) are observed in plots with the lowest level of solar radiation (north exposure). In contrast, the highest RI concerns plots with south exposure, the highest mean slope and within plots partly forested in the first year of observation. This means that plots with southern exposure and high slope were abandoned in the first place. The steep slope was recognised as the most critical determinant of farmland abandonment (Tasser et al. 2007; Lieskovský et al. 2015). In turn, due to lower solar radiation throughout the summer period, northern slopes are characterised by better humidity, which predisposes them to utilisation as meadows and pasturelands (Kostuch 1976). These features could have influenced owners’ decisions as some researchers claim that land-use choices depend on spatial characteristics such as environmental suitability or accessibility (Nagendra et al. 2003; Müller et al. 2013).

By contrast, studies by Janus and Bożek (2019) show that unfavourable spatial features of parcels negatively affecting the profitability of agricultural production have no statistically significant influence on land abandonment decisions. Rutherford et al. (2008) emphasise that defining the spatial determinants of land abandonment using typical variables is difficult because the decision to abandon the land might be a result of socioeconomic factors such as changing rural culture, age and education of farmers, lack of successors, low degree of mechanisation, unfavourable subsidies, lack of motivation and financial profits and interest. Including these factors in the study could improve the results obtained. However, they are not directly available in statistical databases because they concern decisions made several dozen or even more years ago (Gellrich et al. 2007; Lieskovský et al. 2015).

Faster changes in the RI in plots with the lowest level of solar radiation could be connected with better habitat conditions for tree growth, which accelerates the rate of secondary succession. There is a domination of Picea abies preferring even humidity throughout the year in the study area. Even humidity over the year is typical of north slopes with smaller solar radiation (Puchniarski 2017).

Thanks to the use of a series of aerial photographs, the conducted study enabled the factors determining the rate of forest succession to be distinguished from the features of plots characterised by large forest cover. Using information on only current land cover does not provide such opportunities. This could be crucial for areas where land abandonment was spread over time, and the precise date of abandonment is unknown. However, the length of the analysed period is limited by the availability of aerial photographs. In the study area, the oldest aerial photographs are from 1977, so they do not show the largest range of deforestation, and because of that, it is impossible to specify the exact date of land abandonment.

Admittedly, there are archival maps available for the study area after the end of the eighteenth century (Sobala and Rahmonov, 2020). These can provide information on general land cover changes. Still, they could not be used to precisely determine reforestation factors as they were prepared at varied scales, using different map projections and for diverse purposes. Furthermore, there is also a considerable uncertainty related to the completeness and accuracy of map updates. This especially concerns marginal mountain areas where new releases of topographic maps may replicate land-use patterns from the previous map editions (Sitko and Troll 2008; Gimmi et al. 2016). All these issues can influence the applied model (Ligocki, 2001; Plit, 2006). Hence, only mutually comparable aerial photographs were used.

Diverse land use decisions made by individual farmers caused the reforestation process to be chaotic. To cope with this, in this paper, reforestation factors were determined on the level of the plot. Such a solution may be a basis for policymakers to create new policies for these plots, for example, connected with landscape conservation (Czesak et al. 2021). The necessity to study at this level results from the ability to increase the precision of tools and programs preventing the adverse effects of land abandonment (Janus and Bożek 2019). This is needed primarily due to the high probability of the further growth of abandoned areas throughout Europe (van der Zanden et al., 2017).

Some authors have used airborne laser scanning data to map reforestation using crown height model (van Ewijk et al. 2011; Szostak et al. 2014; Kolecka and Kozak 2019). Nevertheless, since various species grow at different rates, abandoned plots may show an irregular vegetation pattern or only part of it may be overgrown. As a result, the crown height model cannot provide information on the date of land abandonment.

Conclusions

This paper shows that in studies on reforestation determinants, three main elements should be considered: the rate of secondary forest succession that causes changes in forest cover, the features of abandoned plots and the time when these plots were abandoned. Distinguishing these three groups is possible thanks to the use of a series of consecutive aerial photographs that present the path of reforestation over time. Such an attitude is crucial for areas where the precise date of land abandonment is unknown. It does not allow the year of abandonment of all the plots to be precisely determined because the length of the analysed period is limited by the availability of aerial photographs that do not show the largest range of deforestation. Nevertheless, it is possible to determine what type of plots was abandoned first.

As the highest average Reforestation Index (RI) in the study area concerns plots with south exposure and the highest mean slope, these are the features of plots abandoned first. These plots also had the highest average RI in the first year of observation, so at that time, the process of forest succession was also developed. On the contrary, the fastest changes in the RI were observed in plots with the lowest level of solar radiation (north exposure).

Diverse land use decisions made by individual farmers caused the reforestation process to be chaotic. In this study, socioeconomic factors were not included. These factors are changing rural culture, age and education of farmers, lack of successors, low degree of mechanisation, unfavourable subsidies, lack of motivation, and financial profits and interest. It is because they are not directly available in statistical databases as they concern decisions made several dozen years ago.

The findings of the studies on reforestation determinants should be a basis for the proper design of legal tools concerning landscape protection and scenarios of landscape changes within marginal mountain areas.