Land Use Modeling Predicts Divergent Patterns of Change Between Upper and Lower Elevations in a Subalpine Watershed of the Alps

The synergic influence of land use and climate change on future forest dynamics is hard to disentangle, especially in human-dominated forest ecosystems. Forest gain in mountain ecosystems often creates different spatial–temporal patterns between upper and lower elevation belts. We analyzed land cover dynamics over the past 50 years and predicted Business as Usual future changes on an inner subalpine watershed by using land cover maps, derived from five aerial images, and several topographic, ecological, and anthropogenic predictors. We analyzed historical landscape patterns through transition matrices and landscape metrics and predicted future forest ecosystem change by integrating multi-layer perceptron and Markov chain models for short-term (2050) and long-term (2100) timespans. Below the maximum timberline elevation of the year 1965, the dominant forest dynamic was a gap-filling process through secondary succession at the expense of open areas leading to an increase of landscape homogeneity. At upper elevations, the main observed dynamic was the colonization of unvegetated soil through primary succession and timberline upward shift, with an increasing speed over the last years. Future predictions suggest a saturation of open areas in the lower part of the watershed and stronger forest gain at upper elevations. Our research suggests an increasing role of climate change over the last years and on future forest dynamics at a landscape scale.

(2050) and long-term (2100) timespans. Below the maximum timberline elevation of the year 1965, the dominant forest dynamic was a gap-filling process through secondary succession at the expense of open areas leading to an increase of landscape homogeneity. At upper elevations, the main observed dynamic was the colonization of unvegetated soil through primary succession and timberline upward shift, with an increasing speed over the last years. Future predictions suggest a saturation of open areas in the lower part of the watershed and stronger forest gain at upper elevations. Our research suggests an increasing role of climate change over the last years and on future forest dynamics at a landscape scale.

INTRODUCTION
Mountain ecosystems and populations around the world are increasingly affected by the combination of changes in climate and land use (Bugmann and others 2007;Lasanta and others 2017). Tropical mountains, for instance, are facing expansion and intensification of agriculture and an over-exploitation of natural resources (Peters and others 2019). In contrast, temperate mountain forest dynamics, especially in Europe, are controlled by land-abandonment processes (Chauchard and others 2007) that, along with gradual effects of climate change, are expected to modify landscape structures and ecosystem services supply in the future (Van der Sluis and others 2019). Nevertheless, it is hard to disentangle the role of climate and land use changes (hereafter, CC and LUC), since they influence in synergy, especially in highly exploited landscapes (Clavero and others 2011).
At lower elevations, LUC plays the most relevant role through secondary successions represented by in-filling of open areas and forest gaps at the expense of abandoned meadows, grasslands and arable lands (Gautam and others 2004;Garbarino and others 2014;Malandra and others 2019). The loss of these areas is strongly related to the abandonment of marginal areas during the twentieth century and to the decline of traditional land uses and practices (Chauchard and others 2007;Tattoni and others 2017).
In the upper part of mountain watersheds, the limit of tree distribution, that is, the treeline, is one of the most studied ecotones where the role of long-term CC can be assessed since it is mainly limited by heat availability (Kö rner 2015;Fajardo and others 2019). However, other environmental variables such as soil and topography also proved to be important treeline drivers (Holtmeier and Broll 2020). Several studies have shown evidence of wood encroachment on upper grasslands (Barros and others 2017;Malfasi and Cannone 2020) and of treeline upward shift in different ecosystems around the world (Fang and others 2009;Elliott 2011;Ameztegui and others 2016). Land abandonment effect is often highlighted as the main driver of the treeline upward shift, but it is expected to be equalized by the CC, which will become the most important in future years (Gehrig-Fasel and others 2007).
Vegetation maps and aerial photographs are fundamental historical data sources for the comprehension of past dynamics and future trajectories of forests. These spatially explicit datasets are useful for LUC detection and can be used to assess the contribution of climate and land use on forest gain (Cousins and others 2015;Filippa and others 2019;Ridding and others 2020). A long time-series dataset with a fine temporal resolution is a crucial aspect for setting observed current changes in a comprehensive historical context in order to produce more realistic future predictions of a process that is highly variable and nonlinear (Becker and others 2007;Tattoni and others 2017).
Many approaches and classification methods can be adopted in LUC forecast modeling. Six types of models exist (sensu Lantman and others 2011): (i) agent-based; (ii) artificial neuron networks (ANNs); (iii) cellular automata (CA); (iv) economics-based; (v) Markov chains (MCs); and (vi) statistical. Since there is not a 'right model' for all ecological purposes, different tools can be integrated in a single framework to boost their strengths and minimize their weaknesses (Mas and others 2014). Another important decision that must be made in forecast modeling is the one between scenario-based simulation models and Business as Usual (BaU) ones. BaU models base their predictions according to the continuation of current land use practices and policies. BaU models are data-driven tools that permit to: (i) forecast landscape trajectories according to past situations and (ii) stress past dynamics. There are two model types that can be considered adequate to carry out this task: Markov chains (MCs) and artificial neuron networks (ANNs). MCs are widely used to model both anthropogenic and natural or semi-natural LUC at different spatiotemporal scales (Muller and Middleton 1994;Tattoni and others 2011;Al-Shaar and others 2021). MCs forecasts are based on previous changes and on the assumption of the persistence of historical dynamics, thus they are reliable tools for BaU scenarios. Their main strengths are simplicity and flexibility and the ability to describe complex and lengthy processes of land use dynamics as simple transition probabilities (Lantman and others 2011; Iacono and others 2015). However, these models are not spatially explicit (Noszczyk 2019), and because of these aspects they are often combined with other models such as ANNs and CAs (Tattoni and others 2011;Ozturk 2015).
ANNs are machine learning algorithms built around several layers of interconnected neurons, similar to human brains (Noszczyk 2019). Therefore, they can recognize patterns and facilitate the development of irregular relationships between past and future LUC (Lantman and others 2011). ANNs are usually connected to suitability maps, expressions of transition potential from one state of the system to another. Their principal limitations are the 'black box' approach and the time required to build them. They have different applications on LUC modeling like urban growth, natural trends, habitat loss (for example, Gontier and others 2010;Tattoni and others 2011;Fattah and others 2021). Though the interaction between MCs and ANNs has been thoroughly analyzed, few studies use more than two or three land cover maps (that is, time steps) to predict future dynamics, oversimplifying the complexity of LUC over time.
In this study, we analyzed past dynamics and future trajectories of LUC using five land cover maps and a nonlinear BaU modeling approach to test our three hypotheses on an inner subalpine watershed: (i) post-abandonment forest gain at lower and upper elevations shows different landscape patterns; (ii) primary successions at higher elevations are becoming more important in the last decades due to increasing climate change and decreasing land abandonment effects; (iii) forest canopy closure at lower elevations and the availability of open areas at the treeline ecotone will favor a stronger forest gain at higher elevations. We had the opportunity to test our hypotheses at the Mont Avic Natural Park, Aosta valley (hereafter, MA), a human-dominated forest ecosystem with a history of intense forest exploitation and a relatively low agro-pastoral impact.

Study Area
MA is located on the southwest part of Aosta Valley, a mountainous autonomous region in northwestern Italy. The Natural Park covers more than 5800 ha and contains Chalamy and Champorcher valleys. It was created in 1989 to protect natural resources in the upper part of the valleys, where high value landscapes had not been strongly modified by humans because of their rough terrain. For this reason, the main human disturbance is not livestock grazing as in many other alpine valleys, but forest cutting for mining activities (Mont Avic 2018). The climate is Alpine, with annual average temperatures ranging from 1 to 3°C and precipitations (800-1200 mm year -1 ) mainly concentrated in autumn and spring (Tiberti and others 2019). The valleys belong to the 'Mont Avic Ophiolitic Complex', and the prevailing lithology is serpentine (D'Amico and others 2008). Deciduous trees such as beech Fagus sylvatica L., chestnut Castanea sativa Mill., downy oak Quercus pubescens Willd. and birch Betula pendula Roth prevail below 1100 m a.s.l. and are sporadic over 1500 m a.s.l. Three coniferous forest species (Larix decidua Mill., Pinus sylvestris L. and Pinus mugo Turra subsp. uncinata) dominate the landscape in the upper montane and subalpine zones (1100-2000 m a.s.l.).
Elevations above 2600 m a.s.l. were excluded from analyses because conditions are considered unfavorable for plant species; in this way, the analyses were limited to the vegetated areas, reducing the weight of unvegetated soil (Garbarino and others 2020). For the analysis of trends and successions, three areas were chosen: (i) the full study area (5478 ha), (ii) the upper elevations (1230 ha) and (iii) lower elevations (4248 ha). The threshold between upper and lower areas was determined based on the maximum elevation of the dense forests in the oldest available aerial image (that is, 1965), corresponding to a medium elevation of 2182 m a.s.l.

Image Analysis and Environmental Predictors
Historical forest dynamics were evaluated by using land cover maps obtained from the classification of five aerial images spanning 52 years (Table 1 Supplementary Materials). Historical aerial photographs were orthorectified and segmented. A supervised classification based on an initial set of polygons was then performed, and lastly a manual classification of residual unclassified polygons (Garbarino and others 2020). Five land cover classes (LCCs) were considered: dense forest (FO), sparse forest (SF), grassland (GR), urban surface (UR) and unvegetated (UV). The land cover class 'dense forest' includes high ( ‡ 80%) canopy cover stands; the 'sparse forest' class represents lower (< 80%) canopy cover stands and shrublands, which have a similar spectral signature and pixel texture and are therefore hard to separate in this landscape. The unvegetated class is a residual class that includes rocks, bare soil, gravel, sand and water. Urban surfaces were not considered in further analyses because our main goal was to explore natural forest dynamics, independent of the marginal expansion of human settlements in the lower part of the valley. Among the maps, overall accuracy (OA) ranged from 78% (1988) to 92% (2017) with a Cohen's Kappa coefficient between 0.71 and 0.89, respectively (Table 1 Supplementary Materials).
Several environmental predictors of land use change (Table 1) were produced such as topographic variables (elevation, aspect, slope, heat load index (HLI), topographic wetness index (TWI), topographic position index (TPI), terrain ruggedness index (TRI), roughness, curvature), anthropogenic variables (cost of movement, Euclidean distance to buildings and roads), the distance from preexisting forest and sparse forest edges and the likelihood of class transitions. The choice of predictors was based on a preliminary literature search and expert knowledge (for example, Rutherford and others 2008; Dubovyk and others 2011; Garbarino and others 2020). A neighborhood size of 8 was used for all the GIS variables computed at the neighborhood level. We derived the accumulated cost of movement through Tobler's hiking on-path function using slope and buildings from Open Street Map as starting points. The cost was expressed as hours required for the movement. Al-most all the predictors were produced in R environment: Cost of movement was implemented with the movecost package (Alberti 2019), the heat load index was calculated with the spatialEco package (Evans 2020), and the topographic wetness index was calculated according to Beven's classical topmodel with the dynatopmodel package (Metcalfe and others 2018). The likelihood of transitions was calculated in the TerrSet environment (Eastman 2016).

Historical Landscape Pattern Analysis
The historical patterns were evaluated looking at both quantity and allocation of changes. The quantification was assessed in the R environment by producing 10 transition matrices regarding all the possible combinations between land cover maps. This procedure was conducted for the full landscape and the upper and lower parts of the landscape (Table 2 Supplementary Materials). The qualitative and spatial dynamics were assessed with landscape metrics, produced with the R package landscapemetrics (Hesselbarth and others 2019). These two analyses allowed a more comprehensive insight into the nature of ecological successions. We considered four landscape metrics for the comparison between higher and lower elevations: edge density, patch density, contagion and Shannon's evenness index. The Shannon's evenness index was preferred over diversity indices because it does not include the richness, which is a meaningless variable in our landscape since the number of classes does not change over the years.

Landscape Pattern Forecast
Future changes were predicted for both short (that is, 2050) and long-term (that is, 2100) predictions using a GIS modeling approach in TerrSet environment within the panel 'Land Change Modeler' (LCM), conveniently implemented with other TerrSet tools. This framework joins two separate models: MC and multi-layer perceptron (MLP) ( Figure 2).
First, a 'transition probability matrix' was produced through the use of MC model that expresses the likelihood of change from one state to another in a consecutive way, by providing two distinct land cover maps (Jokar Arsanjani and others 2013). A MC is a system of elements that pass through one state to another over a discrete time space (Balzter 2000) consistent with the Markov property: For all times t = 1,2,3,… and for all states s = s 0 , s 1 , …, s t , s. Accordingly, X t+1 depends upon X t , but it does not depend upon X t-1 , …, X 1 , X 0 . The transition probability matrix P reports the probability that each land cover type would be found after a certain number of time units for k states: P ¼ p 11 ::: p 1k ::: ::: ::: p k1 :::  Takada and others  2010 for the implementation of yearly matrices).
Second, a spatially explicit model was trained to produce a suitability map that describes the likelihood of change among the different cells. For this study, we applied a MLP model, one of the most common types of ANN (Sangermano and others 2012) that consists of an input layer, one or more hidden layers and an output layer, where every neuron in each hidden layer is connected to other neighboring layers' neurons (Ozturk 2015). MLP is a non-parametric algorithm; thus, it allows multicollinearity and meaningless variables to be excluded.
The model consists of a series of submodels that represents specific transitions. For this study, only the main dense forest dynamics, that is, the transitions from SF, GR and UV to FO, were considered when training MLP. Model parameterization is a crucial part of the process and requires adjusting especially the start learning rate and the hidden layer nodes in order to produce the highest accuracy (Eastman and others 2005). The explanatory variables in the submodels can be both static and dynamic. Static variables do not change over time, while dynamic variables are recalculated during the prediction. In this study, we considered all the variables as static except for the distance from forest and sparse forest areas. It is common in landscape ecology studies (Mishra 2016;Ozturk 2015) to consider distance from buildings and roads as dynamic variables, but since the area did not show great urbanization, these variables were counted as fixed. The MLP model also applies a Jackknife test to measure the relevance of each variable in driving the changes.
The joining of the two models (MLP-MC) allows for the change allocation prediction according to the amount of change (MC model) and the potential for change (MLP model). A common experimental approach is to predict future cover classes based on the oldest and newest land cover maps (for example, Al-Shaar and others 2021; Mirici 2018). In this research, we compared this standard model (hereafter, SM; 1965 as a starting date for the calibration period) with an optimized model (OM) that considered the ability in predicting future land cover for the year 2017. To do this, all possible combinations between maps before 2017 were tested to extract all the change detection matrices (Markov chain's probability transition matrices and area matrices) of LU change (Table 3 Supplementary Materials). Note the numbers in this table refer to probabilities of change according to MC, and consequent areas, while the numbers in Table 2 Supplementary Materials are observed transitions in the past. By using the 'Markov' tool of TerrSet, the probability transition matrices were calculated considering the classification accuracy as a measure of uncertainty called proportional error, which in this study was based on the mean of the OAs of the calibration periods. We then validated these MCs by comparing observed and predicted land cover values for 2017 using Pearson's Chisquare test according to the Eq. (3): where n = 4 (FO, SF, GR and UV classes were considered); O i is the observed value of land cover; E i is the expected one obtained from the MC model. The initial state chosen for this OM was the one with the lowest Pearson's Chi-square in respect of the 2017 prediction. We predicted the 2017 land cover map according to the SM and the OM with reference to 2006 to validate the results (Verburg and others 2004). The process was conducted with the TerrSet tool 'Validate', which provides an assessment of the spatial prediction according to four components of the kappa index of agreement (KIA): Kstandard, Kno, Klocation and Klocationstrata (Table 4 Supplementary Materials). The combination of these four kappa indices allows the overall success rate to be assessed and to understand the strength factors of the prediction (that is, allocation and quantity).
Once this MLP-MC validation had been performed, two BaU predictions were assessed for both the SM (1965-2017 calibration period) and OM (initial state selected according to Pearson's Chi-Square results and final state as 2017), one for the short-term (2050) and one for the long-term timespan (2100). The two models were also averaged to evaluate mean trends. The transition potential matrices obtained with the MC were applied to all the LCCs (Table 3 Supplementary Materials). The future allocation (MLP-MC model) considered only the major transition (all-to-FO) to obtain superior results (Ozturk 2015) and to gather just future forest dynamics. Three class metrics (edge density, patch density, mean core area) based on FO class were performed to compare past and future dynamics.

Historical Landscape Pattern
Historical landscape dynamics were measured with quantitative and spatial metrics, according to transition matrices and landscape metrics. Figure 3 shows the evolution of the five LCCs area from 1965 to 2017 on the total landscape. The most important trend is the expansion of dense forest, from 1778 ha in 1965 to 2978 ha in 2017, corresponding to an increase at a rate of + 1.3% year -1 . The biggest loss was experienced by unvegetated areas, more than 36% of their surface and a decrease at a rate of -0.7% year -1 (from 1876 to 1194 ha). According to the Producer's Accuracy, the most uncertain class is SF, while other classes are relatively accurate. Temporal ( Figure 3) and spatial ( Figure 4) LUC patterns at MA emerged as being not linear. FO increase between 1975 and 2006 was constant (+ 2.0% year -1 in the period 1975-1988 and + 1.5% year -1 between 1988 and 2006), while it decreased in the last 11 years (+ 0.6% year -1 ). SF decreased between 1975 and 2006, but showed an improvement in the last 11 years (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), with an increase rate of + 2.1% year -1 .
At lower elevations, the decreasing trend of edge density, patch density and Shannon evenness and the increase of the contagion index suggest a simplification of spatial pattern due to the forest gain (Table 2). In the upper part of the landscape, the most relevant class over the period remains UV, thus facing a decrease at a rate of -0.4% year-1 in the 50 years (from 1082 ha in 1965 to 891 ha in 2017). SF and GR increased over the period with a rate of + 7.8% year-1 and + 1.6% year-1, respectively (Figure 4b). Looking at the transition matrices (Table 2 Supplementary Materials), UV contributed to increasing especially the GR (168 ha exchanged between 1965 and 2017) and then SF (46 ha). Landscape metrics (Table 2) indicate an opposite trend with respect to the lower elevations, corresponding to fragmentation and plant colonization of UV areas. Biggest changes-represented Past and Future Dynamics of a Subalpine Forest Ecosystem in Figure 4b and Table 2-were experienced in the last years, especially from 2006 to 2017; UV surface decreased at a rate of -1.1% year -1 in this period, more than nine time faster than the previous 40 years -0.14% year -1 , Shannon evenness increase in this period (+ 3.4% year -1 ) was more than ten times higher with respect to the previous period (+ 0.3% year -1 ).

Model Outcomes
The Pearson's Chi-square analyses show that the highest predictive power for 2017 projection was the one obtained from the 1975-2006 period, with a value of 40.1. Thus, we selected 1975 as the initial state for the OM. The mean accuracy of MLP models was about 60% (Table 4 Supplementary Materials). All the kappa values for the 2017 validation (1965-2006 and 1975-2006 calibration periods) are greater than 0.80, reaching 0.88 peaks for Klocation and Klocationstrata in the SM (Table 4 Supplementary Materials). Predictions can thus be considered strong, and both the optimized and standard MLP-MC models can be applied to predict future dynamics. Figure 5 shows the relative importance of the most relevant driving variables common to all the MLP models (distance from FO, the likelihood of class transitions, the cost of movement and the TWI) according to the Jackknife test of MLP models. These four driving factors contribute to more than 99.8% of the total accuracy. There are little  (1965,1975,1988,2006,2017) and area of each land cover class over time. Urban class was removed from the analysis. The error bars represent the Producer's Accuracy of the different classes (Table 1 Supplementary Materials). differences between the models regarding the driving variables explanatory power, and the relevance order is the same.

Landscape Pattern Forecast
Predicted forest gain is higher for the OM (1975-2017 calibration period) than the SM (1965-2017) ( Figure 6). Initial FO expansion (2017-2050 shows an increase at a rate of + 0.29% y -1 and + 0.16% y -1 , respectively, for the OM (3258 ha in 2050) and SM (3139 ha). Between 2050 and 2100, the OM's trend remains quite constant (+ 0.33% y -1 , 3802 ha), while it increases for the SM (+ 0.33% y -1 , 3674 ha). SF class shows an initial growth until 2050 for the two models, then a decrease related to FO transitions, confirmed by the mirrored behavior of FO and SF transitions and by the potential matrices (Table 3 Supplementary Materials), which reflects a great inclination of SF classes toward FO. UV, the second biggest class in 2017, will be outclassed by SF before 2050. GR cover was predicted as decreasing in a similar way according to the two models: about -0.30% y -1 for the period 2017-2050, about -0.55% y -1 for the period 2050-2100.
Dense forest gain was forecast toward upper elevations, especially in the long-term prediction (Figure 7). Class metrics relative to the dense forest stands reflect different future forest dynamics at lower and upper elevations (Table 3). Below the timberline, edge density and patch density were predicted to decrease between 1965 and 2100 (-50% and -80%, respectively), while mean core area was forecast to increase tenfold over the years. Edge density and patch density showed a divergent trend in the upper part of the landscape, while mean core area showed an increase where in 1965 there were no dense forests.

Historical Landscape Pattern
Forest gain following land abandonment is a common dynamic in temperate mountain areas (Kamada and Nakagoshi 1997;Gautam and others 2004;Benayas and others 2007). Considering dense (FO) and sparse (SF) forest classes together, we observed an overall increase at a rate of + 0.56% year -1 (1965-2017 period) that is similar to the mean value recently observed by Garbarino and others (2020) for other landscapes of the Alps and Apennines (+ 0.60% year -1 ) and to the reforestation rate of temperate areas reported by Sitzia and others (2010) of 0.59 ± 0.30% year -1 . The FO class alone shows a higher past overall increase (+ 1.30% year -1 ). The length of time since abandonment affects forest density; thus, the value we observed highlights an important role of land abandonment (Tasser and others 2007;Orlandi and others 2016).
We observed a divergent pattern of change between lower and upper elevation in our study area. Landscape metrics and transition matrices display an increase of homogeneity at lower elevations due to the expansion of dense forests led especially by gap-filling processes to the detriment of sparse forests (SF) and grasslands (GR) (Tables 2, 2 Supplementary Materials). The increase of dense forests started in 1975, while between 1965 and 1975 their amount remained quite stable, and the SF class increased from 1164 to 1353 ha (Figure 3). Looking at the transition matrices, these starting dynamics seem to be linked to a mutual exchange of surface between the two land classes (SF/FO). The observed pattern may be explained by the slowness of initial forestation processes and the occurrence of the final stage of forest use and rural     1965-2006 and 1965-2017) and optimized (1975-2006 and 1975-2017) models. The error bars represent the standard deviation. (that is,, treeline advance and tree density increase), but the upward shift rate was lower (Trant and others 2020). The transition from unvegetated to vegetated classes (primary succession) increased in magnitude in the last decade (+ 15.9 ha year -1 ), compared to the previous period (+ 2.7 ha year -1 ). The acceleration of primary succession processes at MA may reflect the leading role of the CC, which is probably equaling the relevance of LUC and is expected to become the dominant driver of change, especially at upper elevations (Gehrig-Fasel and others 2007; Filippa and others 2019). A particular case of transition in the upper part is the expansion of sparse forests and shrubs, which highlights a fast woody intrusion in the subalpine and alpine zone. The SF class increase is due to the transition from unvegetated soil (46 ha), which can be considered a primary succession, and from grassland (37 ha), a secondary succession (Table 2 Supplementary Materials). Mountain areas with past heavy exploitation followed by abandonment typically show a similar pattern, with gap-filling processes at lower elevations (Malandra and others 2019), and an upward shift of treeline and woody encroachment at upper elevations (Harsch and others 2009;Malfasi and Cannone 2020). Historical land use legacies on forest structure and landscape patterns is a well-known process (Garbarino and Weisberg 2020). New forests that developed in the last 50 years through encroachment on grasslands and meadows show a higher proportion of larch among seedlings, young, and dominant trees layers, whereas mountain pines dominate forest patches developed on unvegetated areas (primary succession).

Landscape Pattern Forecast
Our land cover change models predict a future change that is consistent with the historical one. Gap-filling dynamics will be dominant at lower elevations and a fragmentation will continue in the upper part of the catchment (Table 3). However, the overall forest gain and the expansion of FO class alone are predicted to slow down in the future (2017-2100) according to the two models at a mean rate of + 0.31% year -1 (Figure 6). These predictions are considered realistic because of a reduction of space availability for forest gain at lower elevations and the slow processes at the timberline ecotone. Moreover, future predictions tend to allocate forest gain at upper elevations more than to the detriment of grasslands located in the lower portion; this may be explained by the present utilization of these meadows that do not persuade the MLP-MC model to forecast future forest gap-filling in these areas (Figure 7). A decrease in forest edge and patch densities is predicted by our models. Mean core area, instead, shows different behaviors, with a regular increasing trend for the standard model and a higher increase for the optimized model that suggests strong future gap-filling dynamics causing a further saturation of available open habitats.
At upper elevations, the density of forest patches is predicted to follow two different future patterns: a linear expansion of trees (standard model) and a rapid expansion of small patches followed by a decrease in the number of patches (optimized model), maybe due to a closure of the patches. Since the two models differed only for the starting period, the amount of dense forests of 1965 and 1975 influences the future forest gain prediction. According to long-term future predictions, sparse forests will be replaced by dense ones in our study area (Figure 7, Table 3 Supplementary Materials). However, it is important to consider that the prevailing serpentine lithology of the catchment limits soil fertility and subsequently slows down the growth and establishment of trees (Kim and Shim 2008). Nevertheless, the increase in extreme drought frequency due to CC may counteract the expensive dynamics of forest, leading to different and complex treeline dynamics, ignored by models like MLP-MC (Allen and others 2010). Other extreme events, such as wildfires, heavy rains with potential of landslides, insect disturbances and their interaction, are expected to increase in frequency and intensity according to climate and LU changes and may dampen future predictions reliability (Barros and others 2017;Seidl and others 2017).

Models Discussion: MC, MLP-MC and Driving Factors
LUC historical patterns show a nonlinear trend in our study area, so we adopted a non-linear modeling approach to make predictions. Assessment of future dynamics based on past land use data proved to be a useful tool not only for predicting future trends, but also for amplifying and better understanding past dynamics, especially with a fine temporal resolution. Still, forest landscape ecological forecasting studies often lack a common design and rely on a few maps or low time frequency. Our methodological approach based on the selection of the initial stage with the Pearson's Chi-square highlights that land use dynamics between 1965 and 1975 were different from the recent ones. Furthermore, by looking at general trends it is possible to notice that the slope of transition remains constant from 1975 to 2006 for almost all LCCs. This means that the OM focuses on the most relevant degree of change.
The distance from preexisting forest edges is the most important driver for the forest cover change ( Figure 5). The proximity to closed canopy affects seed recruitment and creates favorable microsites for forest gain in different forest ecosystems (for example, Gü nter and others 2007; Garbarino and others 2020).
The second driving variable according to its relevance is the evidence likelihood, which represents the future transition probability based on past observations. Another important driver appears to be the cost of movement according to the Tobler hiking function that joins anthropogenic and topographic aspects to determine the accessibility of Table 3. land for human activities and is thus a proxy for site remoteness that is strictly related to historical human legacies such as harvesting and pasturing. Prediction of forest landscape changes through historical aerial images and environmental driving variables can be important to satisfy the increasing request of future LUC scenarios to guide decision making (Guan and others 2008;Tattoni and others 2011;Stü rck and Verburg 2017). With this study, park managers can understand future trajectories of forests under BaU scenarios. Our predictions confirm the hypothesis under which climate change and land abandonment promote gap-filling dynamics at lower elevations and woody encroachment at upper elevations, leading to an overall loss of open habitats (Barros and others 2017). This may cause different cascading effects such as loss of biodiversity, increased risk of fire ignition and propagation due to landscape uniformity and fuel buildup, and the loss of cultural landscapes and management techniques (Lasanta and others 2017; Mantero and others 2020). In European mountains, the loss of a-diversity caused by the decline of open habitats is mostly related to least concern species with low vagility, while the past forest gain led to an enrichment in forest specialist birds and mammals (Guliherme and Pereira 2013; Martínez-Abraín and others 2020).

CONCLUSIONS
Our data provide evidence of divergent land use change dynamics between lower and upper elevations in a subalpine watershed of the Alps over the last 50 years. Secondary succession gap-filling processes dominate at lower elevations, whereas primary successions and treeline advancement are stronger at upper elevations and have accelerated in recent years. Climate and land use change effects are difficult to disentangle, but our predictions suggest that the influence of climate will be stronger than land use legacies on future upper elevation forest dynamics. Our results are in line with other studies in mountain watersheds and represent a possibility of a BaU forecasting approach in alpine ecosystems. One caveat to this approach is the assumption that socioeconomic conditions will endure in the future. It is important to remark that different land use scenarios might arise unpredictably because of changes in agricultural and forestry policies (for example, the European Union's Common Agricultural Policy, CAP). In some mountain regions, recent pastoral promotion is leading to an increase in domestic density and traditional land uses (Lasanta and others 2016).
However, it is rather difficult that an abrupt change in socioeconomic patterns will occur in a remote and marginal area such as MA, where the soil conditions have always been a limitation for pastoral practices as a consequence of the toxic and unfertile serpentine lithology (D'Amico and others 2008). Long fine-resolution time series allowed a good evaluation of past dynamics on a subalpine watershed and led to the evidence of the role of climate change at upper elevations. We believe that this approach should be applied to reconstruct time series of land use change over a wide range of mountain forest ecosystems. As shown in this study, a higher number of LC maps over the years can be crucial to assess the role of climate change and land abandonment in alpine landscapes and to predict future dynamics. Moreover, the application of a nonlinear model along with different LC maps allowed a thorough model calibration, avoiding the traditional approach in which the adopted calibration period is given by the oldest and newest land cover classification. The relatively low agro-pastoral activity at Mont Avic and the high classification accuracy of the unvegetated class are important to emphasize the role of climate-driven primary successions. The main limitation of our work is the low accuracy of the SF class, which includes both rare woody areas and shrubs, which can lead to some mistakes in the accuracy assessment and future projection. Climate change is expected to be the most relevant process in the near future, overcoming LUC in alpine forests and our results suggest that it is likely that forest gain at upper elevations will increase in the future, with consequences on habitat composition, biodiversity, natural disturbance regime and forest management. nd_cover_raster_data_of_Mont_Avic_NW_Alps_/1 6718737.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/li censes/by/4.0/.