The area of interest is located in the N-W corner of the autonomous province of Bolzano – South Tyrol, Italy, and corresponds to the Senales municipality that covers the entire Senales valley, 46°42'N–10°55'E (Fig. 1). The Senales valley branches out from the Venosta valley and covers an area of approximately 210 km² with the lowest point being at the beginning of the valley, 824 m a.s.l., and the highest peak being at Similaun at 3607 m a.s.l. (Autonomous Province of Bolzano – South Tyrol 2024). The orientation of the valley is mostly N-W, with a smaller side valley, the Pfosser valley, stretching toward the NE. Based on measurements from the meteorological station in Vernagt/Vernago (Weather South Tyrol 2024), the Senales valley experienced an annual precipitation mean of 713 mm during the period 1981–2023 with a mean of 129 days of precipitation per year. Most of the precipitation occurs during the spring and summer months from April to August. The annual mean temperature registered from 1981 to 2023 for the valley was 5.5 °C. Among the main tree species that can be found at the tree line level in this valley is European larch (Larix decidua Mill.). In the Senales valley, the larch makes up the forest stands all the way up to the tree line on mountain faces with S and W exposure, where it receives more light during sun hours, while on N- and E-exposed slopes, the forest composition is more mixed. Another tree species found in the area is Swiss stone pine (Pinus cembra L.). Here, it can primarily be found on N- and E-exposed slopes towards the tree line. It is considered one of the best adapted trees to withstand cold climates (Caudullo and de Rigo 2016). Patches of Norway spruce (Picea abies (L.) H.Karst.) and broadleaved species are found at lower elevations in the montane zone.

Within the Senales valley, 16 smaller study areas (500 m × 600 m each) were selected as the basis for the analysis (Fig. 1 and Suppl. material 1). Selection of the Senales valley was guided by representing treeline-shift dynamics in the South-Tirol region, the presence of a proximal meteorological station providing reliable climatic data, and the availability of Landsat scenes with minimal cloud cover. The areas selection was then based on three topographic criteria: elevation, aspect, and steepness derived from the Digital Terrain Model (DTM) of the Senales valley (Autonomous Province of Bolzano – South Tyrol 2024).

This categorization resulted in eight possible combinations (four possible aspects, N, E, S, W, and a slope steepness below or above 30°). Based on these criteria, 16 areas were selected in total, ensuring that each set of topographic criteria was represented twice.

For this analysis, we used a Landsat 5 (TM) image from August 29, 1985, and a Landsat 8 (OLI) image from July 2, 2022, both on path 193, row 027. They were downloaded as Level-2 data, which means that they had already been atmospherically corrected to save time and avoid user errors (U.S. Geological Survey 2024). In addition, a DTM and two orthophotos (from 1982/85 and 2023) were provided by the Province of Bolzano. Description of the data (acquisition date, source, and resolution) can be found in Table 1.

Once the satellite images were collected, they were processed using the Semi-Automatic Classification plugin (SCP) (Congedo 2021) in QGIS (version 3.24.1-Tisler (QGIS Development Team 2022)). The bands were stacked into two separate rasters, one containing the Landsat 5 bands and one containing the Landsat 8 bands (Suppl. material 1). The training areas were manually labelled in each study area and assigned to one of the two classes created in our case Forest and Non-Forest (containing grasslands, bare rock, and urban areas). Then, based on the training areas, each pixel in the multispectral image was categorized into the created classes, grouping together pixels based on their spectral signatures (n° of pixel for Forest class = 2633, n° of pixel for non-Forest class = 2246). Among the different algorithms available, the RF algorithm was selected, keeping the default number of trees at 10 and enabling the balanced weight class’s function, as the number of pixels within the training areas for each class was not identical. For every pair of study areas with the same criteria (aspect + slope + year), the same training dataset was used. For each training set, we used 8 training areas for the Forest class and 8–12 areas for the non-Forest class, the latter chosen to capture the wider environmental variability of this class. This sampling strategy ensured a balance between computational efficiency and adequate representation of environmental conditions across the study area. Classification accuracy was assessed using the area-based error matrix method (Olofsson et al. 2014) implemented in the SCP plugin. This approach provides unbiased estimates of user’s, producer’s and overall accuracy. Accordingly, the accuracy of the RF algorithm was evaluated for all Landsat images and for both classes using these unbiased metrics. Furthermore, a pixel analysis of the classification layers was performed to calculate the total forest area within each study area for two years, where the total amount of pixels and the number of pixels corresponding to the forest were counted. This was repeated for the 16 areas for the years 1985 and 2022, and the results were compared. Following the classification, the tree line drawing was performed manually in QGIS separately for each year. First, the tree lines were drawn using the RF classifications for both years. Afterwards, tree lines were drawn directly on orthophotos from 1982/85 and 2023. Since the orthophotos did not have the same spatial resolution, the one from 2023 was down sampled from 20 cm to 1 m to match the resolution of the 1982/85 aerial image. The orthophoto used for 1982/85 is a composite of images from the Autonomous Province of Bolzano, collected between 1982 and 1985 and subsequently merged. Although there is no indication of which specific areas were photographed in which year, we assumed that, at tree line level within our study sites, no major biotic or abiotic disturbances occurred that would have disrupted the continuity of forest cover. After delineating the tree line, points were generated every 30 m along the lines, each carrying elevation information extracted from the underlying DTM. The mean elevation of all points was then calculated, and the results were compared between manual and ML classification. When the positions of the study areas were selected during data preparation, the polygons used did not align perfectly with the pixels of the Landsat images and therefore also contained cut pixels. This resulted in study areas of slightly different size, ranging from 304 to 340 pixels of 30 m resolution, since only complete pixels were included in the classification and following processes. For this reason, only the percentage shift in forest cover will be used for the final comparison.

The meteorological daily data from 1981 to 2023 were provided by the Autonomous Province of Bolzano – South Tyrol through the official Weather South Tyrol platform (Weather South Tyrol 2024). In the Senales valley, a meteorological station is located in Vernagt/Vernago at 1700 m a.s.l., providing both temperature and precipitation records. Although monthly precipitation data have been available since 1953 and monthly temperature data since 1967, the daily dataset was preferred in this study, as it fully covers the study period and allows for a more detailed analysis of climatic trends that could have influenced tree line dynamics. The station was located at an average distance of approximately 5 km from all study sites, with the farthest site 10 km away (coordinates: WGS84 (EPGS: 4326) 46.735658°N, 10.84928°E). Therefore, it was considered the closest and most representative source of climatic data for the study area. Climatic trends were identified using a combination of linear models (LM) and generalized additive models (GAM), performed in R (R Core Team 2014) with the mgcv package (Wood 2024). Although LM offers a robust framework for analysing dependent data, they rely on the assumption of normally distributed conditional responses (Nazeri Tahroudi 2025). GAM offers a flexible framework for modelling non-linear relationships by fitting smooth functions to the data. They are particularly well-suited for capturing seasonal patterns and long-term trends in the occurrence and intensity of daily climatic variables (Underwood 2009; Bassett et al. 2021). The use of both modelling approaches allows for a more nuanced interpretation of the data: LM captures the overall direction of change, while GAM reveals the non-linear temporal structure, offering insights into how the trends have evolved over time. This combined approach strengthens the robustness of the analysis and improves the interpretability of climate dynamics. In this study, trends were calculated for annual mean, minimum and maximum temperatures, as well as for annual precipitation. Furthermore, the analysis was conducted separately for each season to account for potential intra-annual patterns. Daily precipitation data along with daily minimum and maximum temperatures were used to estimate the solid fraction of precipitation in terms of snow water equivalent (SWE) (Carrer et al. 2019). Traditionally, the threshold for distinguishing solid from liquid precipitation was based on the mean daily temperature, typically ranging between +1 °C (Auer 1974; Feiccabrino et al. 2015) and +2 °C (Bartlett et al. 2006; Dai 2008). However, considering only the mean temperature implies that the temperatures exceed that value during the day. By incorporating both minimum and maximum daily temperatures, the estimate aligns more closely with reality, enabling a more accurate evaluation of the solid-to-liquid precipitation ratio. In this SWE estimation scheme, +2 °C is the temperature threshold chosen for both the daily maximum and minimum temperature. Therefore, for each day with a total precipitation amount equal to P and daily maximum and minimum values T_max and T_min, the SWE was established as follows:

SWE ={0 if Tmax>th and Tmin>th( i.e., no snow )P if Tmax<th and Tmin<th( i.e., all snow )P∗(th−Tmin)(Tmax−Tmin) if Tmax>th and Tmin≤th( i.e., half & and half )}

with t_h = +2 °C

Since SWE values are derived from a single meteorological station and it is not possible to apply a robust interpolation method that accounts for elevation gradients due to the lack of a sufficiently dense station network (Brunetti et al. 2014; Crespi et al. 2018; Carrer et al. 2019), they should not be interpreted as absolute values of solid precipitation. Rather, they should be considered representative of the snow trends during the time period considered. Once the SWE was defined, the potential growth season (GS) (its start, end, and length) was assessed based on both daily mean temperature and presence / absence of snow cover. Following the approach of the WATFLOOD model (http://www.civil.uwaterloo.ca/watflood/), also implemented in the TREELIM model (Paulsen and Körner 2014) to identify the minimum thermal requirement for potential tree line position limits, we calculated the snowpack variation under the assumption that snow cover can constrain the actual length of the growing period, regardless of sufficiently warm air temperatures. Subsequently, we extracted the timing and duration of the thermal growing season each year based on the mean daily air temperature. The thermal start of the season (SOS) was defined as the first occurrence of six consecutive days with mean daily temperatures above 5 °C (Frich et al. 2002). The thermal end of season (EOS) was identified as the first occurrence of six consecutive days with daily mean temperatures below 5 °C after July 1. The length of the growing season (LOS) was then calculated as the number of days between SOS and EOS. Once both parameters were integrated to define the potential “climatic” start of the growing season, a LM and a GAM were applied to identify the trend in the set time period.

The results of the Landsat image classification are presented for years 1985 and 2022, respectively (Table 2). Both producer’s and user’s accuracies were greater than 97% for both years considering the individual classes of the binary classification, while the overall accuracies were around 98% for both years.

The classified Landsat images from 1985 and 2022 provide an overview of the response of tree line dynamics in the Senales valley over a period of 37 years. In all 16 study areas, an increase in forest cover was detected between 1985 and 2022, both considering Landsat satellite images and orthophotos (Fig. 2).

The pixel analysis shows a heterogeneous spatial shift in forest cover in the 16 study areas (Table 3). Areas exposed to W and S have experienced the largest increase in forest cover in the considered time period, with an expansion ranging between ≈ +46% and ≈ +104%. Two of the N- facing areas show a mid-magnitude gain of ≈ +31% and ≈ +38% in forest cover, while the remaining two N-exposed areas showed little increases, namely ≈ +16% and ≈ +19%. The E-exposed areas show an average expansion between ≈ +18% and ≈ +25%. The study area with the lowest increase in spatial cover is S1 with ≈ +14%. The overall shift in forest cover is ≈ +44%. The slope appears to have little influence on the change in forest cover, although areas on slopes less than 30° generally have a higher percentage increase in change than steeper slopes in most cases. Study areas W1 and W2 (both W-exposed areas with steepness greater than 30°) show a substantial spatial shift, with 60.48% and 103.67% more forest cover in 2022 compared to 1985.

An overall elevational advancement of the tree lines can be seen in 15 of the 16 study areas based on the Landsat classifications (Table 3). The greatest advancements are recorded in areas exposed to W. The mean upward shifting for these four study areas is ≈ +53 m over the 37-year period. For the other cardinal directions, E, S, and N, the respective mean advancement are ≈ +31 m, ≈ +28 m and ≈ +15 m. Taking into account the areas separately, the pattern also shows that the N-facing areas generally show small upward shifts if not regression (≈-6 m to ≈ +26 m), while the areas with an E-facing slope show mid to large advancements (≈+20 m to ≈+44 m). Interestingly, S-facing areas show a relatively large advancement similar to that of E or W (≈+42 m to ≈+44 m), or a relatively small upward shifting similar to the N-facing slopes (≈+13 m to ≈+17 m).

For tree lines derived directly from orthophotos from 1982/85 and 2023, only advancements were registered (Table 4). Similarly to the Landsat-derived tree lines, the western exposure reveals the biggest advancements but with an increased magnitude. The mean tree line advances for the four aspects are ≈+77 m, ≈+38 m, ≈ +36 m and ≈ +28 m for W, N, S, and E, respectively. Since the spatial resolution of the orthophotos is 1 m, they are considered truer to the real-life values compared to the 30 m resolution of the Landsat images. Due to shadows present in the orthophotos, an accurate tree line drawing was not feasible in all areas.

A comparison between the two data sets (Table 4) reveals a considerable difference in tree line advancements in some areas. Areas N2 and W1 show the biggest mismatch between the changes registered with Landsat and orthophotos, whereas areas such as E4, S2, and S3 show very similar advancements. The tree lines from orthophotos generally show a greater elevational advance compared to those derived from Landsat, that is, an average shift of ≈+45 m compared to ≈+32 m. Only in 4 areas the difference between orthophotos and Landsat tree lines is greater than 30 meters, corresponding to 1 pixel in the Landsat images, and only two of them are greater than a two-pixel difference (marked in red in Table 4). Both data sets show a diverse magnitude of the advancements between the study areas, even with the same aspect and slope.

This section details the climatic trends for temperature, precipitation, and derived variables from 1981 to 2023, to understand the possible climatic drivers in the upward shifting of the tree line.

The annual mean temperature over the considered period shows a significant trend (p-value < 0.01), indicating that the temperature pattern is not linear (Fig. 3A). When LM is applied before and after the direction change, at the point of inflection, a significant decrease in mean temperature can be observed from 1981 to 2000 (approximately -0.07 °C per year), followed by a significant annual increase of 0.05 °C from 2000 to today. By contrast, annual precipitation did not show significant trends according to either model, highlighting its complex and irregular behaviour compared to the other climatic variable (Fig. 3B). Furthermore, considering seasonal mean temperatures (Fig. 4), spring has experienced a significant overall increase, approximately 0.04 °C each year (Fig. 4A), and together with summer (Fig. 4B) and autumn (Fig. 4C), they show a significant non-linearity in the trend. Although spring observed an increase in temperature since the early 1990s, summer only shows a change in the trend ten years later together with the mean temperatures in autumn. Maximum temperatures have decreased significantly in winter (Fig. 4D) and autumn (Fig. 4C) (approximately -0.05 °C and -0.04 °C per year, respectively), while during the warmer months, no significant increase in maximum temperatures was observed. Changes in the direction of the non-linear trend were observed earlier for autumn temperatures (early 2000) and almost 10 years later for winter temperatures. Minimum temperatures show significant trends across all seasons, with notable increases of 0.04 °C, 0.04 °C, 0.02 °C, and 0.04 °C for each year in winter, summer, autumn, and spring, respectively. In the case of spring, the increase is so pronounced that the inflection point in temperature trends is virtually absent, indicating constant increments with no setbacks during the 2000’s.

Regarding seasonal precipitation, only summer experienced a significant increase in the last 37 years, with +1.62 mm per year, resulting in a total increase of approximately 60 mm during the studied period (Fig. 5). This is also the only trend that shows two direction change points, reflecting the higher variability of the pattern considering precipitations. No significant trends were observed in the other three seasons (Suppl. material 1).

Finally, we identify the potential start and end of the growing season alongside the cumulative depth of snow, measured as SWE (mm). No significant trends were detected in any of these variables, particularly regarding the growing season length, which is intrinsically linked to both its start and end dates. In addition, the depth of the snow did not show a discernible trend during the study period, reflecting a pattern characterized by high variability and a lack of consistent directional change. This behaviour appears to be closely associated with the stochastic nature of total precipitation throughout the considered time period (Suppl. material 1).

This study demonstrates that the combined use of satellite data from the Landsat program and ML algorithms can play a key role in ecological research that addresses land-cover changes on large spatial and temporal scales. Integration of historical aerial imagery further supports and improves classification accuracy, providing valuable context for long-term landscape dynamics. The combination of Landsat data, GIS, and the RF algorithm has been successfully used for other LULC assessments (Pham et al. 2024) and exhibited a strong classification accuracy in handling complex land cover classes and acting as an efficient method to assess landscape dynamics. Although the treeline appears to follow a similar pattern when considering both the manual delineation and the classified forest boundaries, some differences in the identification of the ecotone limit were observed. Transition ecotones, such as tree lines, are likely to exhibit significant differences in the outcome classification of forest cover, often due to groups of trees showing considerable variations in appearance, even within the same area of interest (Nguyen et al. 2022). This leads to differences in the mapping continuity in the context of photointerpretation. The precision and reliability of manual interpretation are strictly influenced by the experience of the annotators, who detect and interpret geographical features in the images based on their professional knowledge (Tarko et al. 2021). However, it is labour intensive and time consuming, and the consistency of the results is also subject to operator fatigue (Borghuis et al. 2007; Svatonova 2016). Moreover, the heterogeneity of an area due to disturbances or topographic variability remains difficult to distinguish, especially for supervised classifiers, where the high variability of the spectral signal cannot be detected well in the land use/land cover (LULC) application, particularly in mountain regions (Shimizu et al. 2023). In heterogeneous landscapes, ML algorithms can be easily misled, as their reliance on spectral information alone makes it challenging to accurately separate neighbouring land-cover classes (Beaubien 1986). The observation that orthophotos better capture site complexity may also reflect their substantially higher spatial resolution compared to Landsat imagery, in addition to the greater interpretative skill of the annotators. For example, for example, in the case of area N2, the orthophotos better capture the complexity of the site and deliver results that are truer to real life. For the remaining study areas, the deviation between the tree lines derived from Landsat and orthophotos was less than 30 meters, corresponding to less than one pixel in the Landsat images. As the overall difference was 13 m, which is not considered a large deviation, RF was confirmed as a suitable, time-saving and accurate methodology also for highly heterogeneous landscape classification.

The observed temperature and precipitation trends are in line with those of Hansson et al. (2023), who investigated the link between seasonality and tree line migration in various zones of the northern hemisphere. Taking seasonal trends separately into account is fundamental, as variations in temperature and precipitation could have different effects on the phenology of trees and therefore affect their recruitment and growth. For example, opposing summer and winter trends may counteract each other when viewed annually, or the mean and summer temperature could show a divergent trend if studied separately (Zheng et al. 2021). The measurable effect of temperatures on the advancement of the tree line is opposed, as in our case, to the lower correlation with precipitation patterns. The significant increase in minimum temperatures appears to be related to the advancement of the tree line, while maximum temperatures do not contribute to migration, but rather to the growth performances of already settled trees during the growing season (Hansson et al. 2023). Warmer autumn temperatures can extend the growing season in cold-limited areas such as high-elevation tree lines, extending the duration of the necessary conditions for tree establishment and growth (Coops et al. 2013; Paulsen and Körner 2014). Thus, the timing of meeting the thermic criteria, in spring and autumn, is more important than the actual temperatures of the rest of the year. Even if in our case no significant trend has been observed, the end of the growing season seems to occur progressively later during the year. Furthermore, the pronounced increase in spring minimum temperatures appears to have contributed to the tree line advances. This trend not only extends the duration of the growing season but also reduces the probability of frost damage in plants, as higher mean temperatures are generally correlated with a lower probability of frost events (Kharuk et al. 2009; Wipf et al. 2009). However, during extreme frost events, trees are more susceptible to frost damage, particularly in the absence of a protective snow layer that would otherwise protect seedlings and newly emerged leaves early in the growing season. Despite increasing maximum and minimum temperatures in summer, which could potentially reduce the odds of tree line migration (Hansson et al. 2023) due to the higher rate of evapotranspiration (Coops et al. 2013) and the increased probability of drought stress caused by the lack of available soil water, in the Senales Valley, the significant increase in summer precipitation may have maintained soil moisture levels sufficient to prevent trees from experiencing water stress.

The aspect also seems to have a considerable impact on the different magnitudes of the upward shift of the tree line. In agreement with the study by Tourville et al. (2023), the greatest tree line advancements were observed on W-exposed slopes. As a result of higher-intensity solar radiation on southern slopes, trees in these areas face two direct consequences: drier condition due to soil water deficit and higher chance of photoinhibition (oxidative damage) from more direct sunlight, reducing their recruitment and growth potential (Holtmeier and Broll 2007; Doeweler 2021), while for exposed areas in the north, the limiting factor was assessed to be mean lower temperatures (Wang et al. 2021). Based on data from the Global Wind Atlas (Davis et al. 2023), the prevailing wind direction at 10 m height – the level that most influences tree growth and mechanical stability – exhibits a bimodal pattern along a NW-SE axis, reflecting the valley-parallel thermal circulation typical of alpine valleys. In addition, both wind frequency and intensity are dominated by flows from the N-NW sector, likely associated with föhn-type downslope events. As a result, S to SE-facing slopes are consistently exposed to the strongest and most energetic winds, and this could also be the reason why the E exposure in Senales valley seems to show limited extent in tree line advancement. Wind exposed locations are cooler than those protected from wind (Holtmeier and Broll 2010), which could influence delays in phenological development during the growing season (Telewski 1995). Leaf desiccation (Cairns 2001) and leaf abrasion of the wax layers by wind-driven ice and mineral particles during winter (Dahms 1992) could both lead to increased water loss rates in evergreen species, undermining their survival rate at tree line level. In addition, wind speeds greater than 3 m s⁻¹ reduce height growth in young trees (Kronfuss and Havranek 1999). The steepness of the slope in the Senales valley seemed to have no clear influence on the advances of the tree line. Even if some studies suggest that steeper terrain is more likely to be subjected to disturbances due to intensified geomorphological processes, which can limit further tree growth at high elevations (Holtmeier and Broll 2012; Leonelli et al. 2016), it has been proven that these factors occur mainly in areas with high slope angles (>45°), where repeated high-energy geomorphic events are likely to occur (Leonelli et al. 2016). Most of the hillsides in the Senales valley are quite steep, but also far from constant within the individual study areas. Also, our division of the areas used 30° as the threshold between steep and non-steep zones, so it might not be adequate to draw any conclusions based on this.

We are aware that working with intermediate coarse resolution satellite data in areas of complex topography and vegetation cover can lead to ambiguous results in the classification, as smaller variations are not detectable. Therefore, in a mountainous area, scarce tree cover in the tree line ecotone can cause errors or misinterpretations (Wei et al. 2020). Finer-resolution multiband satellite data such as Sentinel-2 could provide deeper insights into land cover changes, particularly in complex and mixed land-use contexts (Babitha et al. 2025). Sentinel-2 imagery has been successfully used not only for land cover and land use classification (Topaloğlu et al. 2016; Forkuor et al. 2017), but also to assess forest ecological variables such as canopy cover (CC) and leaf area index (LAI) (Korhonen et al. 2017), as well as for the detection of herbaceous and tree species (Immitzer et al. 2016). Most comparative studies between Sentinel-2 and Landsat 8 have reported equal or superior performance of Sentinel-2-based models, largely due to the inclusion of the red-edge band (B5), which has shown strong correlations with forest parameters estimation (Chrysafis et al. 2017; Astola et al. 2019). Satellite images cannot assess tree height, so it is not clear whether the tree lines derived from Landsat data follow one definition or another. Moreover, common criteria for defining forest presence or absence typically involve not only tree height, but also canopy density and structural form (e.g., the presence of shrubs). Since Landsat imagery cannot capture all these structural characteristics, the classification becomes more uncertain, especially at tree line level, where the boundaries between the two classes are inherently more fragmented. Furthermore, what can be considered as forest needs to cover a certain area to be detected as such within the Landsat image. Subsequently, this can cause a spatio-temporal delay in the assessment of the tree line as existing trees (potentially already settled in that position for a long time) are falsely registered as new, advancing individuals (Luo and Dai 2013; Garbarino et al. 2023). In recent decades, new advances in technology have led to Light Detection and Ranging (LiDAR) RS, which enables the measurement of discrete tree stand structures in three dimensions (Zhang et al. 2017; Sharma et al. 2021; Torresani et al. 2024a). However, terrestrial and airborne applications are limited by high costs, long processing times, and restricted availability of free user data (Gupta and Sharma 2022). Since 2018, the Global Ecosystem Dynamics Investigation (GEDI), a spaceborne LiDAR sensor developed by NASA has provided high-resolution 3D observations of the Earth’s forests (Dubayah et al. 2020; Potapov et al. 2021). GEDI data have been successfully used to derive information not only on forest canopy height (Marselis et al. 2022; Moudrý et al. 2024b; Kacic et al. 2025), but also on others key ecological parameters of forest ecosystems, such as forest biodiversity (Torresani et al. 2023; Moudrý et al. 2024a), above-ground biomass (Dubayah et al. 2022; Saarela et al. 2022), and topographic heterogeneity (Pracná et al. 2025). Thus, the integration of technologies such as GEDI with satellite imagery could provide more precise and near-real-time delineation, thereby improving the accuracy and reliability of treeline boundary detection. Future studies might apply larger training datasets, as well as advanced RS technologies such as GEDI to improve the robustness of the performance of the RF algorithm and detangle the complex spatio-temporal relationship between forest dynamics and climatic drivers. Another important aspect is having a certain knowledge about the vegetation cover within the area studied beforehand to avoid classifying shrubs, such as dwarf mountain pine (e.g. Pinus mugo Turra) as trees. However, dwarf pine is not widespread in the Senales valley. Finally, another important factor is the abandonment of traditional land use practices in high mountain regions. Practices such as cattle grazing and tree logging might have previously modified or pushed the tree line to lower elevations than the climate would allow and suppressed the re-establishment of tree seedlings (Bello-Rodriguez et al. 2019; Mienna et al. 2022). Therefore, advances in forest revegetation could potentially also occur independently of warming temperatures and are rather due to changes in land use. Also taking into consideration anthropogenic disturbances would have provided a broader understanding of the tree line dynamics since the 1980s and possibly indicated strong human influence in some areas. So, investigating historical data on land use in the area would add knowledge on the influence of human activities on current and historical tree lines.