Development and Application of Earth Observation Based Machine Learning Methods for Characterizing Forest and Land Cover Change in Dilijan National Park of Armenia between 1991 and 2019

: Dilijan National Park is one of the most important national parks of Armenia, established in 2002 to protect its rich biodiversity of ﬂora and fauna and to prevent illegal logging. The aim of this study is to provide ﬁrst, a mapping of forest degradation and deforestation, and second, of land cover/land use changes every 5 years over a 28-year monitoring cycle from 1991 to 2019, using Sentinel-2 and Landsat time series and Machine Learning methods. Very High Spatial Resolution imagery was used for calibration and validation purposes of forest density modelling and related changes. Correlation coefﬁcient R 2 between forest density map and reference values ranges from 0.70 for the earliest epoch to 0.90 for the latest one. Land cover/land use classiﬁcation yield good results with most classes showing high users’ and producers’ accuracies above 80%. Although forest degradation and deforestation which initiated about 30 years ago was restrained thanks to protection measures, anthropogenic pressure remains a threat with the increase in settlements, tourism, or agriculture. This case study can be used as a decision-support tool for the Armenian Government for sustainable forest management and policies and serve as a model for a future nationwide forest monitoring system.


Context
The Republic of Armenia (RA) is located at the junction of the biogeographic zones of the Lesser Caucasus and the Iranian and Mediterranean zones and exhibits both a great range of altitudinal variation from 375 m to the 4090 m peak of Mount Aragats, as well as a diversity of climatic zones. This combination has resulted in a variety of landscapes and ecological communities with a distinct flora and fauna, including many regionally endemic, relict, and rare species. Therefore, in the country's small territory of 2,974,300 ha, there are about 3800 species of high vascular plants, 428  The first Forest Code of Independent RA that was accepted in 1994 was rather conservative and logically connected to the former Armenian Soviet Socialist Republic (SSR)

Aim and Objectives
The objective of the present paper is to assess and report on the state of Armenia's forests with a focus on the Dilijan National Park (See Figure 1). This park is one of the four National Parks in Armenia, located in the Northern part of Armenia on the slopes of the Pambak, Miapor and Aregani mountain chains and on the basin of Aghstev and Getik rivers. Although the area became a specially protected area in 1958 in order to protect the park's flora and fauna, the National Park was only established in the year 2002. Today, the Dilijan National Park is known for its forest landscapes, rich biodiversity, medicinal mineral water springs, natural and cultural monuments, and extensive network of hiking trails. We proved the suitability of validated EO products for an independent forest monitoring and support land use planning for sustainable forest management. Two main tasks have been performed: (1) providing detailed land use and land cover classification maps for the epochs 1991-2019 as a basic planning support tool for characterizing forest ecosystems and main evolution trends and drivers of deforestation in Dilijan National Park, Armenia, and (2) delivering an in-depth examination of the forest classes of the previous service, subdividing forested areas into different forest density and type classes. The service shall focus on assessing and locating areas where deforestation or forest degradation is taking place.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 25 National Parks in Armenia, located in the Northern part of Armenia on the slopes of the Pambak, Miapor and Aregani mountain chains and on the basin of Aghstev and Getik rivers. Although the area became a specially protected area in 1958 in order to protect the park's flora and fauna, the National Park was only established in the year 2002. Today, the Dilijan National Park is known for its forest landscapes, rich biodiversity, medicinal mineral water springs, natural and cultural monuments, and extensive network of hiking trails. We proved the suitability of validated EO products for an independent forest monitoring and support land use planning for sustainable forest management. Two main tasks have been performed: (1) providing detailed land use and land cover classification maps for the epochs 1991-2019 as a basic planning support tool for characterizing forest ecosystems and main evolution trends and drivers of deforestation in Dilijan National Park, Armenia, and (2) delivering an in-depth examination of the forest classes of the previous service, subdividing forested areas into different forest density and type classes. The service shall focus on assessing and locating areas where deforestation or forest degradation is taking place.

Study Area
Armenian forests are among the most threatened ecosystems, with degradation accelerating, largely attributable to deforestation and overexploitation. This results in high rates of erosion, increasing soil salinity, lowered soil fertility, and loss of biodiversity. Thus, "expansion of forests is one of the main goals for Armenia, not only for the forests'

Study Area
Armenian forests are among the most threatened ecosystems, with degradation accelerating, largely attributable to deforestation and overexploitation. This results in high rates of erosion, increasing soil salinity, lowered soil fertility, and loss of biodiversity. Thus, "expansion of forests is one of the main goals for Armenia, not only for the forests' protective role, but also to develop forest-related businesses and to ensure fuelwood supply to the population", said Ekrem Yazici, Deputy Chief of the Joint UNECE/FAO Forestry and Timber Section. However, due to a dense population living within the protected forest areas, developed infrastructures, uncontrolled tourism, illegal logging, poaching and non-sustainable Remote Sens. 2021, 13, 2942 4 of 24 use of natural resources, environmental degradation threatens this unique biodiversity and the rich natural-historical and cultural landscapes. To become more resilient to external and internal shocks, there is a necessity for the integration of new approaches and policy instruments to rehabilitate degraded forests and increase forest cover significantly, while formulating the country's development agenda. In this context, UNDP Armenia is focusing its efforts to better understand the past forest ecosystems transformations, the land use and land cover changes, and in general, all the socio-environmental processes in the past and today that affect the sustainable management of forest resources.

In Situ Data
Armenian counterparts have provided three sets of data sources named: (a) vector data on basic cartographic layers-Dilijan NP border, hydrography, roads and railways, settlements, forest ecosystems status data on compartment level (species, density, site class, age group), red listed species (fauna and flora); (b) raster data-AlosPalsar 12 m resolution digital terrain model from 2007 and high resolution satellite images (Pléiades) for 2015, 2018-2019 years and (c) cadastral maps for the 2013-2015 period of seven communities in Dilijan NP. All the vector data, Alos Palsar raster data and cadastral maps of communities were provided by the Ministry of Environment of RA as a part of the Dilijan NP management plan for the 2007-2011 and 2017-2025 periods [6] p.215. High resolution satellite images (Pléiades) were procured within the UNDP-GEF "Mainstreaming sustainable land and forest management for the north-eastern mountain landscapes of Armenia" full size project (2016-2020) [12].

EO Data
As the analysis of land use and land cover and associated changes over the Dilijan national park spans from 1991 to 2019, covering 28 years, it was decided to use two satellite constellations: Landsat and Sentinel. As part of the Landsat family, Landsat 4, 5 and 7 times-series at 30 m spatial resolution were used for the mapping of the park in 1991, 1995, 2000, 2002, 2005 and 2010. Landsat 8 times-series at 30 m spatial resolution were used for 2015 and Sentinel-2 images at 10 m spatial resolution were used for the most recent mapping of 2019 (see Figure 2).
Two satellite VHSR data were used for calibration and validation purposes:

Data Pre-Processing
For Landsat data, the Level-2 Surface Reflectances have been processed using the ESPA on demand service. A total of 909 Landsat 4/5, 685 Landsat 7 and 323 Landsat 8 have been ordered. Using QA data (cloud detection), mosaics have been processed and cropped to the Dilijan National Park border.
For Sentinel data, the Level-2 flat Surface Reflectance using Maccs-Atcor Joint Algorithm (atmospheric correction and cloud screening software based on the MACCS processor, developed for CNES by CS-SI, from a method and a prototype developed at CESBIO) have been processed and downloaded from CNES Peps platform. After a resampling of the 20 m and 60 m bands to 10 m resolution, mosaics of Sentinel-2 imagery have been processed and cropped to the Dilijan National Park border. Remote Sens. 2021, 13, x FOR PEER REVIEW 5 of 25

Classes and Coding
Forest Density

Forest Mask
The Forest Mask is not a Forest Product, but it serves as a basis for the Land Cover/Land Use classification and associated changes, as well as for the Forest Density, the Forest Type and the Forest Degradation/Deforestation maps. The Forest Mask processing workflow can be divided into 6 main steps: (i) Pre-processing of EO data, (ii) Classification, (iii) Post-processing, (iv) Computing of the raw Forest Change Mask for all epochs, (v) Manual Enhancement of polygons of changes, (vi) Quality Check of the consistency of the Forest Masks and Forest Change Masks for all epochs.
(i) Pre-processing First, atmospheric, radiometric, and topographic corrections, using the MAJA algorithm for Sentinel-2 data, were applied on the time-series available, followed by a cloud masking. High Resolution Landsat and Sentinel-2 images for epoch n ± 1 in enlarged vegetation season from 1st of May to 31st of October were then selected in order to obtain a full mosaic coverage of the study area.
(ii) Classification of the Forest Mask Machine Learning methods based on the random forest algorithm is a common practice found in literature for the mapping of land use and land cover, and more specifically of forest areas [9,10].
The classification methodology is based on Geographic Object-Based Image Analysis (GEOBIA) using the Broceliande processing chain [13,14]. This methodology has proven its efficiency for the mapping of the Copernicus High Resolution Layers Small Woody Features and Forest for the reference years 2015 and 2018, at Pan-European scale covering an area of about 6 million km 2 [8,9]. The approach is actually a mixture of both GEOBIA and pixel-based analysis. It relies on two main components which are feature extraction based on attribute profiles, and a semi-supervised classification using a random forest algorithm at pixel level. The first step consists of computing predefined indexes such as the Normalized Difference Vegetation Index (NDVI), as well as texture characterization such as Sobel gradient and Haar-like features based on integral image representations. The Sobel gradient has proven to be more efficient than the renowned Haralick (1973) texture feature extraction based on a Gray Level Co-occurrence matrix (GLCM) which requires a time-consuming computation of a set of statistics for multiple distances and orientations. The texture information is extracted from each spectral band of the original EO multispectral image, as well as from the NDVI. This textural information will have a significant role in allowing the discrimination between Built-up and Bare soil, as their spectral response can be close. A first binary mask is produced, where all pixels flagged as no-data are set to a NoData value. Multi-scale features called Attribute Profiles (AP) and Differential Attribute Profiles (DAP) are derived through a model of morphological tree, which can be seen as a hierarchical segmentation: these trees allow to identify, depending on the way they are pruned, morphological objects from the EO data (see Figure 4). The second step is the use of the random forest classifier at pixel level.
Reference samples were collected, i.e.,~30 forest samples (examples) and 30 nonforest samples (counterexamples) to train the classifier. An in-house Machine Learning Remote Sens. 2021, 13, 2942 8 of 24 object-oriented classification algorithm was used to derive a probability map, depicting the probability for each pixel to belong to the forest class. A manual thresholding approach enabled the conversion of the latter into a binary forest/non-forest layer.
significant role in allowing the discrimination between Built-up and Bare soil, as their spectral response can be close. A first binary mask is produced, where all pixels flagged as no-data are set to a NoData value. Multi-scale features called Attribute Profiles (AP) and Differential Attribute Profiles (DAP) are derived through a model of morphological tree, which can be seen as a hierarchical segmentation: these trees allow to identify, depending on the way they are pruned, morphological objects from the EO data (see Figure  4). The second step is the use of the random forest classifier at pixel level.
Reference samples were collected, i.e., ~30 forest samples (examples) and 30 non-forest samples (counterexamples) to train the classifier. An in-house Machine Learning object-oriented classification algorithm was used to derive a probability map, depicting the probability for each pixel to belong to the forest class. A manual thresholding approach enabled the conversion of the latter into a binary forest/non-forest layer.  Considering Forest Mask 2019 as the reference and the most precise product based on Sentinel-2 data at 10 m spatial resolution, the consistency of the changes between each epoch was checked retrospectively. Based on Forest Mask products generated for the 8 requested epochs (1991,1995,2000,2002,2005,2010, 2015, 2019), forest changes were isolated and analyzed and served as a basis for the final forest degradation and deforestation map (See Figure 12).

Forest Density Mapping
The Forest Density is defined as the vertical projection of tree crowns to a horizontal earth's surface and provides information on the proportional tree canopy coverage per pixel, in the range of values from 1 to 100%. The method used for the estimation of Forest Density is similar to the one used for the Copernicus Land Monitoring Service High Resolution Forest at Pan-European level (https://land.copernicus.eu/pan-european/highresolution-layers/forests, accessed 23 July 2021) and REDD+ (Reducing Emissions from Deforestation and Forest Degradation) Copernicus project, for example in African countries (https://www.reddcopernicus.info/, accessed 23 July 2021). Forest Density modelling is based on the multiple linear regression between reference samples and vegetation indexes derived from the HR times-series available. Forest Density reference data is assessed by means of visual interpretation of VHSR data following a point grid approach. The processing workflow of the Forest Density map can be divided into 4 main steps: (i) Sample drawing and interpretation, (ii) Computing of vegetation indexes, (iii) Multiple linear regression modelling, (iv) Quality Check of the consistency of the Forest Density between all epochs.

(ii) Computing vegetation indices
The second step consists in computing several vegetation indexes within each PSU over the HR times-series, in order to establish a correlation between those indexes and the Forest Density value obtained through the interpretation of samples. Common indexes for the assessment of Forest and the characterization of vegetation and leaf properties are the Normalized Difference Vegetation Index (NDVI), Brightness Index (BI), Normalized Burnt Ratio (NBR), Green Normalized Difference Vegetation Index (GNDVI), Modified Normalized Difference Vegetation Index (GNDVI), Modified Normalized Difference Water Index (MNDWI), as well as the mean value of the spectral bands available, such as Blue, Green, Red, Near Infrared (NIR) and Short Wave Infrared (SWIR) [15].
Normalized Difference Vegetation Index (NDVI) Green Normalized Difference Vegetation Index (GNDVI) Modified Normalized Difference Water Index (MNDWI) (iii) Multiple linear regression analysis The third step lies in defining the model that will predict the Forest Density value, through an analysis of the best multiple linear regression model, according to the following equation: Y = a 1 x 1 + a 2 x 2 + a 3 x 3 + . . . . +a n x n + b (6) Y = estimated Forest Density value a n = coefficients x n = explanatory variables (vegetation indexes) b = constant Candidate models are assessed with the use of the Bravais-Pearson correlation coefficient, and the model with the higher coefficient was selected.
x = mean σ x = standard deviation Finally, the most accurate model was applied over the complete area of interest in order to derive the Forest Density values. An example of Forest Density for the most recent 2019 epoch is shown in Figure 8.
(iv) Quality Check of the consistency of the Forest Density between all epochs The Computed Forest Density value over a stable area can present change over 1991-2019 due to several factors: for instance, EO imagery from various sensors or acquisition time of EO image reflecting seasonality. To avoid artificial changes in Forest Density value, a smoothing of the Forest Density was applied in order to ensure the consistency between the different epochs. Small variations of values between epochs were then considered as artefacts unless visually interpreted as real changes on EO data. For those stable areas, Forest Density values were then averaged between two successive epochs.

Forest Types
According to the referenced literature [3], Dilijan NP mainly consists of deciduous species such as oak (Quercus iberica, Q. macranthera), oriental beech (Fagus orientalis), common and oriental hornbeam (Carpinus betulus, C. orientalis), which form homogeneous oak, beech and hornbeam forests as well as mixed forests with different combinations of the species mentioned. Coniferous forests (pine-Pinus, juniper-Juniperus and yew-Taxus) occupy a limited territory in the national park and occur in patches. Pine often makes dense forests in the basin of the River Hovajur on the slopes of the Areguni and Pambak ranges in the vicinity of serpentine Dilijan highway. There are many pine trees in Dilijan and on nearby slopes. Additionally, there are no Larix sp., neither Quercus Sempervirens at all in Armenia. Buxus is grown only artificially for greening purposes in settlements.
The Forest Type map is the distinction between broadleaved and coniferous forest. Broadleaved Density and Coniferous Density completes this map as it provides the density of broadleaved and coniferous forest at 10 m spatial resolution for the 2019 epoch. To derive those maps, the processing workflow is divided into 3 steps: (i) Identification of coniferous reference samples and analysis of their spectral signature over the EO data time-series, (ii) Classification over the whole area, (iii) Crossing with Forest Density to derive Broadleaved Density and Coniferous Density.
(i) Identification of coniferous reference samples and analysis of their spectral signature over the EO data times-series Using UNDP ancillary field data provided in the frame of this study, coniferous parcels could be identified in 2007. Visual interpretation over the 2019 Sentinel-2 scenes, using these 2007 parcels as reference, could allow the determination of training samples of both broadleaved forest area and coniferous forest area inside the 2019 forested area. NDVI was computed for 5 Sentinel-2 scenes over the 4 annual seasons, to establish "NDVI signatures" of both tree types (see Figure 6 below). (i) Identification of coniferous reference samples and analysis of their spectral signature over the EO data times-series Using UNDP ancillary field data provided in the frame of this study, coniferous parcels could be identified in 2007. Visual interpretation over the 2019 Sentinel-2 scenes, using these 2007 parcels as reference, could allow the determination of training samples of both broadleaved forest area and coniferous forest area inside the 2019 forested area. NDVI was computed for 5 Sentinel-2 scenes over the 4 annual seasons, to establish "NDVI signatures" of both tree types (see Figure 6 below). (ii) Classification over the whole area Using the NDVI signatures derived in the previous section, it was possible to derive, inside the 2019 forest mask, the discrimination between broadleaved and coniferous trees. (ii) Classification over the whole area Using the NDVI signatures derived in the previous section, it was possible to derive, inside the 2019 forest mask, the discrimination between broadleaved and coniferous trees. Using several observations allowed the reduction in artefacts that could derive from shadows.
Given the very limited area covered by coniferous species, the quality of the Forest Type map was assessed by means of visual interpretation over reference field data provided by UNDP Armenia.

Specifications
To support land use planning in the Dilijan National Park, the production of standardized land use and land cover (LULC) maps that show the current baseline and highlight past changes and trends was necessary. The production was divided into eight epochs-for each of these epochs a dedicated map was produced. Table 2 lists all of these epochs as well as the EO input data used for production. The target resolution of the final maps corresponds to the resolution of the EO data that was used as an input for the production of each year. However, for the most recent year (2019) two separate thematic maps were produced: Firstly, a map in 30 m to achieve better comparability to previous epochs and secondly, a 10 m map to facilitate a potential future Sentinel-2 based monitoring system of the national park, which would require a 10 m thematic map as a baseline. Furthermore, the minimum mapping unit (MMU) of objects on the ground was chosen based on the production year as well. The year 2019 was specified to have an MMU of 0.25 ha and all other years (based on Landsat) have an MMU of 1 ha.
Finally, the thematic content of the maps was specified in accordance with the stakeholders and covers the following classes: 1. Forest 2.

Land Use and Land Cover Classification Method
The LULC classification workflow can be divided into several parts: (i) Pre-processing, (ii) Feature Extraction, (iii) Training, (iv) Classification and (v) Post-processing. These steps are subsequently described in more detail below: (i) Pre-processing The pre-processing includes the cloud masking of all EO input data and a reprojection into the same geographic grid and projection UTM 38N. For Landsat data, the cloud masking is performed by using the provided quality layer (BQA) which provides information in clouds, cloud shadow, cirrus, snow and ice and cloud confidence. For Sentinel-2 the MAJA cloud masks were used to mask out saturated or defected pixels, dark areas, snow, thin cirrus clouds as well as clouds with high and medium probability and their shadows.

(ii) Feature Extraction
The feature extraction first includes the derivation of meaningful indexes that will help distinguish vegetated and non-vegetated classes. For this reason, for each available image, the Normalized Difference Vegetation Index (NDVI) as well as the Normalized Difference Water Index (NDWI) was calculated. After this step of feature engineering, the data were temporarily aggregated: For each spectral band and year, three percentiles (p) were derived: p20, p50 (median) and p80. This temporal aggregation is used to reduce the amount of data for the model and therefore decrease training and processing time as well as for the purpose of creating a less complex model (with less input features) with the capabilities of better generalization.

(iii) Training
The training of every supervised classification task requires reference data. These are already to some extent described in Section 2.2. However, two key datasets that were essential to the LULC classification are described here in more detail. This on one hand includes the Global Crop Extent Layer which was produced by Pittman et al. (2010) [16] and serves as the only available input dataset on agricultural activities in the area. Secondly, the Open Street map dataset which offers a very high accuracy and level of completion in the area of interest and serves as training data for most thematic classes: Settlements, Roads, Water bodies and rivers. The other classes (Bare Soil and other vegetation) were sampled manually with the help of very high-resolution data as well times-series information from Sentinel-2 imagery.
(iv) Classification The classification itself was performed with the highly popular tree based Random Forest classifier [17], the training parameters (although often neglected in case of random forest) were tuned via cross validation and overall accuracy as a target metric. The final model with 500 trees yielded an out-of-bag (OOB) error of 14.6% for the classification of the most recent year 2019 based on Sentinel-2 data.
(v) Post-Processing The subsequent post-processing included an overlay and insertion of thematic information from Open Street Map and Service 2 with the aim of enhancing the accuracy of the product. Overall, it was decided to take over the following thematic classes from Open Street Map: Settlements, roads, water bodies and rivers. This was done for the reason that a human-based manual very detailed delineation of these object was judged to be superior to the artificial intelligence-based classification based on high resolution imagery. The Forest class was entirely taken over from Service 2. Furthermore, the MMU was applied to the classification, however it has to be stressed that for Settlements and Primary Roads the MMU was lowered to 0.02 ha. This step was necessary to preserve the small-scale scattered structure of these anthropogenic objects in the area. A larger MMU would yield an underestimation of these areas.
Please note that the above elaborated workflow only describes the classification task for the baseline epoch 2019. All other epochs and their maps were produced with a subsequent change detection approach that compares the spectral characteristics of each year to the baseline year, flags noticeable changes in the area and subsequently reclassifies them into the correct class. This procedure is described in more detail in the next chapter.

Change Detection
To assess historical land cover changes, a spectral change vector analysis can be performed if the sensor and band properties are consistent over time (i.e., the Landsat Legacy). Using this method, spectral properties of the reference year 2019 can be compared to its historic spectral properties for each epoch and for each pixel, respectively. If a spectral difference between two compared epochs for a certain pixel lies above a predefined threshold, a potential change is identified and flagged.
A high change vector length value is, for example, expected during urbanization activities (i.e., the conversion of vegetation into sealed surface). Low change vector lengths may occur if there is no change at all and the spectral properties only differ due to meteorological effects, i.e., because of a delayed or shortened rain season or also during vegetation regrowth due to afforestation activities.
The applied change vector analysis is not based on a single comparison of the spectral properties of two spectral bands but considers the comparison of the entire spectrum and can therefore be called a multidimensional change vector analysis.
As already highlighted above, the algorithm considers each epoch of the reference period (1991-2019) and compares it to the reference year (2019). Ultimately, change areas are detected by dynamically thresholding each comparison individually via a high percentile value (i.e., 0.99). Figure 7 depicts this approach in a two-dimensional space where small change vector lengths are considered stable and high differences are flagged as change. By summing up the occurrences of abnormally high differences over time per pixel, a probability of change or even an exact timing of change can be derived if of interest. A high change vector length value is, for example, expected during urbanization activities (i.e., the conversion of vegetation into sealed surface). Low change vector lengths may occur if there is no change at all and the spectral properties only differ due to meteorological effects, i.e., because of a delayed or shortened rain season or also during vegetation regrowth due to afforestation activities.
The applied change vector analysis is not based on a single comparison of the spectral properties of two spectral bands but considers the comparison of the entire spectrum and can therefore be called a multidimensional change vector analysis.
As already highlighted above, the algorithm considers each epoch of the reference period (1991-2019) and compares it to the reference year (2019). Ultimately, change areas are detected by dynamically thresholding each comparison individually via a high percentile value (i.e., 0.99). Figure 7 depicts this approach in a two-dimensional space where small change vector lengths are considered stable and high differences are flagged as change. By summing up the occurrences of abnormally high differences over time per pixel, a probability of change or even an exact timing of change can be derived if of interest. The identified "change candidate areas" are then flagged and reclassified into target change classes via a supervised classification. The algorithm itself is a supervised stratified random forest classifier as described above. The identified "change candidate areas" are then flagged and reclassified into target change classes via a supervised classification. The algorithm itself is a supervised stratified random forest classifier as described above.

Product Validation
The validation of the final land cover maps is based on the good practices described by Olofsson et al. (2013) and Congalton et al. (1991) [18,19]. Therefore, the sampling design was chosen in a way that each class sample size is large enough to produce sufficiently precise estimates of the area of the class, for that reason the minimum sample size per stratum (class) was set to 30. The equation below calculates an adequate overall sample size for stratified random sampling that can then be distributed among the different strata. An illustration of a Forest Density map is given in Figure 8 below for the last epoch 2019. Dilijan National Park is composed mainly of dense mountainous forest >75%. Despite the harmonization of Forest Density values between different epochs explained in the methodology Section 3.1.3, the latter are strongly correlated to the specific multiple regression linear model used for each epoch and the corresponding EO imagery used to calculate these models. At the beginning of the monitoring cycle, in the 1990s, a strong forest density loss reveals uncontrolled logging and increasing anthropogenic pressure on the natural landscape. The middle period right after the creation of the Park 2002-2005 shows almost a complete stop of clear-cuts and the overall picture is a stable area, corresponding to protective measures taken by the Armenian Government. Real losses in forest density changes observed for the latest 2015-2019 period account for a rebound of forest degradation and deforestation. Despite this threat, and although the aim of this study did not include reforestation or afforestation, it is an encouraging sign to observe that there is also a recent increase in forest density which coincides with natural forest regeneration all across Dilijan National Park.
The validation of the Forest Density products comes into the form of a correlation between reference values (interpretation of samples on VHSR EO data) and modelized Forest Density values present in the corresponding product, as show in Figure 9 below.
shows almost a complete stop of clear-cuts and the overall picture is a stable area, corresponding to protective measures taken by the Armenian Government. Real losses in forest density changes observed for the latest 2015-2019 period account for a rebound of forest degradation and deforestation. Despite this threat, and although the aim of this study did not include reforestation or afforestation, it is an encouraging sign to observe that there is also a recent increase in forest density which coincides with natural forest regeneration all across Dilijan National Park. The validation of the Forest Density products comes into the form of a correlation between reference values (interpretation of samples on VHSR EO data) and modelized Forest Density values present in the corresponding product, as show in Figure 9 below. The correlation coefficient between reference and map Forest Density values for each epoch is presented in Table 3 below.   The correlation coefficient between reference and map Forest Density values for each epoch is presented in Table 3 below.

. Characterization of Forest Types
Analysis of 2019 EO data times-series allowed stating that Dilijan National Park forests are mainly constituted of broadleaved tree species (98.85% of Forested areas) and that coniferous species only represent a minor part of the forest included in this park (1.15% of forested areas) (see Figure 10 below).

Characterization of Deforestation and Forest Degradation between 1991 and 2019
Results concerning the forest perturbation (e.g., deforestation, degradation, or reforestation) over the study period can be found in Table 4 and Figures 11 and 12 below. Results concerning the forest perturbation (e.g., deforestation, degradation, or reforestation) over the study period can be found in Table 4 and Figures 11 and 12 below.

Characterization of Deforestation and Forest Degradation between 1991 and 2019
Results concerning the forest perturbation (e.g., deforestation, degradation, or reforestation) over the study period can be found in Table 4 and Figures 11 and 12 below. With the creation of the National Park in 2002, the forest degradation drastically dropped, and the forest regeneration increased, following the set-up of forestry management policies.
The forest degradation, which is stable over the 2002-2019 period, has been attributed to anthropogenic origin unless ancillary data can prove that the cause is natural. Therefore, it is very likely that the anthropogenic forest degradation is strongly overestimated: forest wildfire monitoring service (NASA FIRMS) only monitors wildfires since year 2000, and no field data regarding other natural causes (storm, disease, etc.) were available.

Land Use and Land Cover Change
The result of the LULC classifications as well as the mapped changes are subsequently presented below. This includes the presentation of the final LULC maps, the comparison of the 10 m and 30 m products, the highlighting of areas with prominent changes as well as the presentation of derived land cover class areas for each year. With the creation of the National Park in 2002, the forest degradation drastically dropped, and the forest regeneration increased, following the set-up of forestry management policies.
The forest degradation, which is stable over the 2002-2019 period, has been attributed to anthropogenic origin unless ancillary data can prove that the cause is natural. Therefore, it is very likely that the anthropogenic forest degradation is strongly overestimated: forest wildfire monitoring service (NASA FIRMS) only monitors wildfires since year 2000, and no field data regarding other natural causes (storm, disease, etc.) were available.

Land Use and Land Cover Change
The result of the LULC classifications as well as the mapped changes are subsequently presented below. This includes the presentation of the final LULC maps, the comparison of the 10 m and 30 m products, the highlighting of areas with prominent changes as well as the presentation of derived land cover class areas for each year.

Main Land Use and Land Cover Types
The results from the 2019 baseline LULC mapping based on Sentinel-2 (10 m) can be seen in Figure 13. In total, eight different classes were mapped, and the map shows that the area is dominated by forested areas with agricultural activities around several smaller settlements which function as a pressure regarding deforestation. The total area breakdown is presented in Figure 14 and shows that natural vegetation makes up roughly 93.7% in the area of interest, with more than half of that area covered by forest. The remaining part of the area is mostly anthropogenically altered or used to some extent: 4.62% is dedicated to agriculture, 1.19% of the area is sealed (settlements and primary roads) and 0.46% are made up of standing water and rivers. The total area breakdown is presented in Figure 14 and shows that natural vegetation makes up roughly 93.7% in the area of interest, with more than half of that area covered by forest. The remaining part of the area is mostly anthropogenically altered or used to some extent: 4.62% is dedicated to agriculture, 1.19% of the area is sealed (settlements and primary roads) and 0.46% are made up of standing water and rivers.
The total area breakdown is presented in Figure 14 and shows that natural vegetation makes up roughly 93.7% in the area of interest, with more than half of that area covered by forest. The remaining part of the area is mostly anthropogenically altered or used to some extent: 4.62% is dedicated to agriculture, 1.19% of the area is sealed (settlements and primary roads) and 0.46% are made up of standing water and rivers. As already mentioned above in the methodology section, two separate LULC maps for the year 2019 were produced at 10 m and 30 m resolutions from Sentinel-2 and Landsat 8 imagery, respectively. The Sentinel-2-based map was created to include the best available data (in the case of Sentinel-2) to produce a high-resolution map as a baseline for potential future monitoring activities. The Landsat 8 map was generated to achieve a better comparability between the preceding epochs, which are all based on Landsat 8 (because the first Sentinel-2 satellite was launched in 2015). Figure 15 shows both classifications side by side highlighting the area around the city of Dilijan inside the national park. On the left side, the Sentinel-2-based product at 10 m is shown and when directly compared As already mentioned above in the methodology section, two separate LULC maps for the year 2019 were produced at 10 m and 30 m resolutions from Sentinel-2 and Landsat 8 imagery, respectively. The Sentinel-2-based map was created to include the best available data (in the case of Sentinel-2) to produce a high-resolution map as a baseline for potential future monitoring activities. The Landsat 8 map was generated to achieve a better comparability between the preceding epochs, which are all based on Landsat 8 (because the first Sentinel-2 satellite was launched in 2015). Figure 15 shows both classifications side by side highlighting the area around the city of Dilijan inside the national park. On the left side, the Sentinel-2-based product at 10 m is shown and when directly compared to the Landsat 8 (right side) at 30 m resolution, it can be clearly seen that Sentinel-2 offers much greater spatial detail. On the left side even small (groups of) buildings are visible and due to the smaller minimum mapping unit, a more fragmented and detailed overview of the area is discernable. to the Landsat 8 (right side) at 30 m resolution, it can be clearly seen that Sentinel-2 offers much greater spatial detail. On the left side even small (groups of) buildings are visible and due to the smaller minimum mapping unit, a more fragmented and detailed overview of the area is discernable. This subsection presents the results of the LULC maps for the epochs 1991-2019. Since the area of interest is mostly inside the Dilijan National Park where conversation practices are put in place, the area is very stable in terms of land cover change. Figure 16 depicts the area breakdown for four selected epochs (1991,2002,2010 and 2019) and clearly shows that most class areas remained almost the same over the entire period of time. The biggest change can be observed in the forest class that dropped from 50.18% in 1991 to 49.71% in 2019. The main driver behind this loss of forest is either natural degradation, small scale urbanization or minor agricultural expansion.
The class "other vegetation" (than forest) increased by 0.41% over the entire period and is the main driver behind the loss in forest, this is followed by a minor increase in agricultural activity of 0.04% and lastly by the extension of settlements which grew by 0.01% of the total project area from 1991 until 2019. All other classes (water and bare soil)

Mapping of Changes from 1991 to 2019
This subsection presents the results of the LULC maps for the epochs 1991-2019. Since the area of interest is mostly inside the Dilijan National Park where conversation practices are put in place, the area is very stable in terms of land cover change. Figure 16 depicts the area breakdown for four selected epochs (1991,2002,2010 and 2019) and clearly shows that most class areas remained almost the same over the entire period of time. The biggest change can be observed in the forest class that dropped from 50.18% in 1991 to 49.71% in 2019. The main driver behind this loss of forest is either natural degradation, small scale urbanization or minor agricultural expansion.
The class "other vegetation" (than forest) increased by 0.41% over the entire period and is the main driver behind the loss in forest, this is followed by a minor increase in agricultural activity of 0.04% and lastly by the extension of settlements which grew by 0.01% of the total project area from 1991 until 2019. All other classes (water and bare soil) did not show any significant change.  Tables 5-7 below and range between 85 and 89% in terms of overall accuracy. Generally speaking, the selected method together with the input data yields good results and most classes show both high user and producers' accuracies (usually above 80%). The class with the lowest accuracies is "Other vegetation" whose producers' accuracy only ranges from 58%-63%. When looking at the validation matrices, it can be seen that this class often gets confused with agriculture. This probably arises from the fact that the training data for agriculture do not stem from a very reliable source (Global Crop Extent Layer) and the very  Tables 5-7 below and range between 85 and 89% in terms of overall accuracy. Generally speaking, the selected method together with the input data yields good results and most classes show both high user and producers' accuracies (usually above 80%). The class with the lowest accuracies is "Other vegetation" whose producers' accuracy only ranges from 58%-63%. When looking at the validation matrices, it can be seen that this class often gets confused with agriculture. This probably arises from the fact that the training data for agriculture do not stem from a very reliable source (Global Crop Extent Layer) and the very broad class definition "Other vegetation" seems to be too vague for the classifier to be able to find sharp class boundaries which yields high confusion with other classes.