The Assessment of Favourability and Spatio-Temporal Dynamics of Pinus Mugo in the Romanian Carpathians Using GIS Technology and Landsat Images

Abstract

Favourability classification for forest species represents a fundamental activity for deriving technological solutions in forestry, as specialists need detailed information about the ecological requirements of forest species from environmental factors: climate, pedological characteristics and morphometric characteristics of the study area. The purpose of the present study was the use of the qualitative data extracted from the ecological records of the Pinus mugo species and the generation of a complex geospatial database for the entire territory of Romania. The results were represented by a collection of thematic maps generated on favourability classes for the Romanian Carpathians, as well as for the major landform subunits which had been the basis for the statistical analysis of the results. The validation of the results was performed by comparing the results obtained through the application of the model which used the frequency points reported in the European Atlas of the Forest Tree Species from Europe, 2016. In order to identify the spatio-temporal dynamics, LANDSAT satellite images from 30 years were used, which enabled the identification of the expansion and the reduction in size of the Pinus mugo area at a zonal level, a process which is dependent on natural factors, like climatic variations, or anthropic factors (overgrazing or works of cleaning the montain pastures).

1. Introduction

The dwarf mountain pine is the most important wood species from the alpine vertical zone of the Romanian Carpathians, due to its stabilising effect on the ecosystems in which it is present and to its supporting role in the shaping of the alpine landscapes.

In the Romanian Carpathians there are 50,000 ha of terrain where the dwarf mountain pine is spatially distributed [1], in the mountain massifs from the subalpine vertical zone from the Oriental Carpathians: Maramureș Mountains, Rodna Mountains, Călimani Mountains, Giumalău Mountain, Ceahlău Mountain, Gurghiu Mountains, Bistrița Mountains, Nemira Mountains, Ciucaș Mountains, Piatra Mare Mountain, Piatra Craiului Mountains, and Postăvarul Mountain; in the Meridional Carpathians: Bucegi Mountains, Făgăraș, Iezer-Păpușa, Căpățânii, Parâng, Retezat, Țarcu-Godeanu, Cernei and Semenic Mountains; and in the Occidental Carpathians: Vlădeasa Mountains (Figure 1) [1].

When classifying habitats, Pinus mugo is included in habitat number 4070 (bushes with Pinus Mugo and Rhodendron Myrtifolium in the sites of the NATURA 2000 network from Romania) [2], which is included in the category of alpine and subalpine habitats with a community interest included in the project LIFE05 NAT/RI/000176—Priority alpine, subalpine and forestry habitats in Romania.

Habitat number R3105—Southeast Carpathian dwarf pine bushes (Pinus mugo) with rhododendron (Rhododendron myrtifolium)—also includes boreal alpine species like Rhododendron myrtifolium (syn. R. kotschyi) and Calamagrostis villosa, which represents a local characteristic for the Romanian Carpathians as they are oligotherm, hygrophilic, oligotrophic, and acidophilic species [2].

Over the time, different studies have been made of the different types of vegetation developed in Pinus mugo associations in the mountain areas, the focus being laid on identifying the geographical and ecological conditions which influence the different vegetation associations [3,4,5] and not on the spatio-temporal dynamics of the previously mentioned habitats.

As a consequence of its adapting to environmental conditions, Pinus mugo is considered a pioneer species [6,7]. It can be found at high altitudes in the east and southeast of the Alps, in the central and northern part of the Dinaric Mountains, in the Rila, the Pyrenees and the Carpathian Mountains where it grows on large areas and as isolated species in the Abruzzy Mountains from the Apennines, in the Jura and the Vosges Mountains [8,9,10,11,12]. It is a species which is used to delineate the mountain forestry vegetation zone from the mountain grass zone [13,14,15,16,17,18].

The modification of this species’ zone was also highlighted by a series of studies made at the level of the mountain area from Western and Central Europe [19], which highlighted the fact that the inferior limit of the dwarf mountain pine in the Apennines had increased from 35 to 65 m and there is a prognosis of a migration with another 30 m until 2060 as a consequence of summer grazing as a main factor [20] and of the increase of average annual temperature as a secondary factor.

The study area is represented by the Romanian Carpathians which are extended on 66,286 km², representing 27.8% of the Romanian territory which has a total of 238,391 km². The Romanian Carpathians represent the main territorial structure of the country as they are positioned on both sides of the Transylvanian Depression and the plateau and plain landform units are disposed concentrically around them. Due to their morphometric characteristics which refer to altitude, slope angle, slope exposure, etc., the mountains determine a vertical zonation of the Romanian vegetation from the steppe vegetation zone, which develops up to 200 m, to the oak and sessile oak zone (200–600 m), the beech zone (600–1200 m), the fir tree zone (1200–1800 m), the dwarf mountain pine zone (1600/1800–2200 m) and the alpine grass vegetation zone (at over 2200 m) [21,22,23,24].

The changes occurring in the spatial extension of the mountain areas with Pinus mugo can be identified with the help of LANDSAT satellite images which have a sufficiently large temporal extension to cover a representative temporal interval: 1988–2018 and a spatial resolution which offers enough details to make a statistical evaluation of the surfaces in order to identify their expansion or reduction. Thus, the remote sensing techniques become extremely useful instruments considering the objectives of the present study which enable the indirect monitoring of the territories with Pinus mugo. The purpose of this study is to classify the Carpathian area from Romania according to the favourability for the growth of Pinus mugo as well as to identify the pretability of using GIS technology, spatial analysis and remote sensing to attain this purpose.

The relationships between climatic conditions and stationary conditions (forest sites) are synthetically represented in the ecological chart of each species [25], which attempts to give the optimum-suboptimum-limitative area according to the determining ecological factors (

Table 1. Determining ecological factors for the favourability of Pinus mugo.

	Ecological Factor	Favourabily		Interval for the Study Area
		Very Low	Low
1.	Average annual temperature (°C)	<2 °C; >3 °C	−1–3 °C	−2.3–10.4 °C
2.	Average annual precipitation (mm/year)	<900 >1400	901–1400	557–1620
3.	Length of bioactive period (months)	>4	1–4	0.6–6 months
4.	Relative atmospheric humidity -July	<75	>75	66–91
5.	Wind (Beaufort degrees)	1,2,3,4,5	6,7	3.9–36.7 km
6.	Altitude (m)	<1300	>1300	56–2544 m
7.	Slope exposure	Flat surfaces	All exposures	All exposures
8.	Saturation degree of bases in soil (V%)	>45	<45	0.002–82
9.	Soil acidity (pH in water)	5.4–7.2	5.4–3.6	4.3–8.2
10.	Azote supply (V = 0.01)	>4	<4	1.06–7.11
11.	Edaphic volume (m3/m2)	>0.45	<0.45	0.17–1.67
12.	Water supply capacity (%)	<12; >20	10.1–20	1.8–20.9
13.	Soil compactness	Very compact	Very loose, loose, moderatey loose, compact	All compactness classes
14.	Soil texture	Loamy-clayey, Clayey-loamy, clayey	Sandy, sandy-loamy, Loamy-sandy, loamy	All texture classes

).

Morphologically, Pinus mugo is a shrub species with stems that may lay close to the ground or be erect and branched (Figure 2), which is adapted to the harsh climatic conditions from above 1600 m, forming pure associations in areas with favourable atmospheric and soil humidity (shady valleys, snow accummulations).

It is not a requiring species in relation to soil trophicity, however, a good development of the dwarf mountain pine is dependent on the soil characteristics determined by the underlying rock, its pH values or its humidity degree [26]. It represents the best fixing factor of eroded slopes and debris cones, this species having the best hydrologic control from the alpine catchments. It is characterised by a slow growth, it makes fruits annually but at various intensities, and the saplings are sensitive to diseases in the first years while their aerial part has a very slow growth.

The ecological formula has the following form: Alpine-subalpine, oligotherm, nivophil, heliophil, oligobasic, oligotrophe, acidophil, mezohygrophit, subsaxicol [27].

At present, in favourability studies for forest species of economic and ecological importance, the literature in Romania offers only general information, ranking by favourability classes only by synthetic indicators, tables that, at the time of realization, have not been correlated with the spatial distribution at level subgroups of the Carpathians. In the present study we propose to use modern GIS technology for spatial analysis so that the classification on favourability classes is achieved at the level of the entire Carpathian chain, allowing the statistical analysis of the results, the evaluation of the correlations between the factors and the satellite images from time to time to analyse the temporal distribution of the studied species.

2. Database and Methodology

In order to identify the favourability classes from the Romanian sector of the Carpathians, to characterise the specific habitat of Pinus mugo, it is necessary to start from a spatial digital database. This must include the most important ecological, climatic, pedological and morphometric factors, which individually and cumulatively influence the various degrees of favourability or restrictivity for the analysed species (

Table 2. The database used.

Nr.	Database	Database Type	Format	Source/Obtaining Process
1	Weather stations	primary	Vector (point)	Vectorisation
2	Average annual precipitation (mm)	modelled	Raster (GRID) 30 m	Implementation in GIS environment of the equation of the regression function Y = 0.435 X + 533.51, where: Y—average annual precipitation (mm) X—relative altitude (m)
3	Average annual temperature (°C)	modelled	Raster (GRID) 30 m	Implementation in GIS environment of the equation of the regression function Y = 0.395 X + 363.76, where: Y—average annual precipitation (mm) X—average annual temperature (degrees C)
4	Number of days with average temperature >10 °C	modelled	Raster (GRID) 30 m	Implementation in GIS environment of the equation of the regression function Y = 205.29 – 0.075 X, where: Y—average temperature >10 °C (days) X—altitude (m)
5	Relative atmospheric humidity in July (%)	derived	Raster (GRID) 30 m	Interpolation of punctual values from the weather stations
6	Wind regime	derived	Raster (GRID) 30 m	Interpolation of punctual values from the weather stations
7	Saturation degree of bases in soil (V%)	derived	Raster (GRID) 30 m	Interpolation of punctual values from the field
8	Soil acidity	derived	Raster (GRID) 30 m	Interpolation of punctual values from the field
9	Azote supply	derived	Raster (GRID) 30 m	Interpolation of punctual values from the field
10	Edaphic volume	derived	Raster (GRID) 30 m	Interpolation of punctual values from the field
11	Capacity of water supply	derived	Raster (GRID) 30 m	Interpolation of punctual values from the field
12	Soil compactness	derived	Raster (GRID) 30 m	Interpolation of punctual values from the field
13	Soil	primary	Vector (polygon)	Vectorisation
14	Soil texture	primary	Vector (polygon)	Vectorisation
15	Digital Elevation Model (DEM)	primary	Raster (GRID) 30 m	Direct acquisition
16	Slope exposure	derived	Raster (GRID) 30 m	Derived based on the DEM

).

The digital database consists of primary data represented by vector data in the case of weather stations, which were generated by digitisation, as well as soil classes, soil types and texture, which were also digitised from the soil map of Romania, 1:200,000.

In the category of primary raster databases we included the digital elevation model with a resolution of 20 m, which was used to model various other databases (the average annual precipitation grid, average annual temperature, number of days with average temperature >10 °C) and derive other databases (slope exposure).

The category of derived grid databases includes climatic characteristics like wind regime and relative atmospheric humidity, while from the pedological characteristics it includes the base saturation degree (V%), soil acidity, etc. (Table 2), obtained by interpolating the punctual values obtained in the field.

The satellite images belong to the Landsat 5, 7 and 8 satellite missions and were downloaded from the US Geological Survey (http://glovis.usgs.gov), so as to cover the time interval 1987–2018, and were selected from the autumn period (August–October) for a better differentiation of the coniferous species in the mountain area (

Table 3. Acquisition data of the Landsat images.

Acquisition Date/ID	Mission	Band
		RED	NIR
2018 12 X LC08_L1TP_184029_20181012_20181030_01_T1	LANDSAT 8 OLI	Band 4	Band 5
2013 5 X LC08_L1TP_185029_20131005_20170429_01_T1	LANDSAT 8 OLI	Band 4	Band 5
2008 15 X LE07_L1TP_185029_20081015_20161224_01_T1	LANDSAT 7 ETM+	Band 3	Band 4
2003 31 VIII LE07_L1TP_185028_20030831_20170124_01_T1	LANDSAT 7 ETM+	Band 3	Band 4
1998 28 X LT05_L1TP_185029_19981028_20180211_01_T1	LANDSAT 5 TM	Band 3	Band 4
1993 28 IX LT05_L1TP_185029_19930928_20170116_01_T1	LANDSAT 5 TM	Band 3	Band 4
1987 30 X LT05_L1TP_185029_19871030_20180211_01_T1	LANDSAT 5 TM	Band 3	Band 4

).

The monitoring of the coverage degree of mountain surfaces with Pinus mugo was performed through the NDVI index, which had previously produced very good results on territories of different extension [19,28,29,30,31,32,33].

NDVI values is tightly correlated with the variations in average annual temperature of the analysed territories as well as with the general trend of temperature increase at a global level, due to global warming [28,29,30,31,34,35,36,37,38].

This index was used as a mapping tool in areas affected by forest fires where the tree density was unsatisfactory after the destructive event [39,40,41,42,43,44,45]. It was also successfully used in regional studies at the level of coniferous forests in order to identify their territorial extension [46,47,48], even at the level of phenophase. Very tight connections have been identified between the values of the NDVI index and the air temperature, as well as a moderate correlation between the NDVI values and the average precipitation amount (in this case influencing the water reserve from the soil) [48].

The database preparation stage represented by the spectral bands of downloaded satellite images on https://earthexplorer.usgs.gov assumed their radiometric correction. The radiometric correction of satellite imagery was done using the Erdas 2013 software based on two Hase and noise reduction models to reduce the impact of atmospheric aerosols and water molecules that lead to solar radiation absorption and scattering, also leading to an improvement in its recording by remote sensors.

Using the extension Map Algebra, Raster Calculator of the ArcMap 10.6 software (ESRI, Romania), the value of the NDVI index was calculated for each year from the analysed period, with a time interval of five years, by applying the formula:

NDVI = (NIR – RED)/(NIR + RED)

where NDVI—Normalized Difference Vegetation Index; NIR—near-infrared; RED—Red band.

The NDVI index (Normalized Difference Vegetation Index) was obtained as a ratio between the difference and the sum of the spectral signatures from the limit of the visible spectrum (RED) with the infrared (NIR—near infrared). The values thus obtained vary between −1 and 1. Based on these values one can identify the areas with trees and shrub vegetation from the mountain area, knowing the fact that healthy vegetation (clorophyll) reflects closer to infrared (NIR) and the green light, in comparison to other wavelengths, and absorbs more red and blue light.

The research takes the form of a complex model of GIS spatial analysis which follows the stages of deriving a prognosis and spatial prediction model based on databases, expert knowledge and dissemination through analysis of the final results.

Thus, the model was structured in three main stages: database acquisition, spatial analysis, model validation and presentation of the final results (Figure 3). In the database acquisition stage, two methods were used. The first was based on the direct acquisition of numerical data (temperature, precipitation) and of the vector databases (soils) and raster databases (the digital elevation model and the Landsat images downloaded from specialised sources—http://glovis.usgs.gov). The indirect acquisition of the main databases which are included in the spatial analysis model is represented by the development of submodels as a result of integrating directly acquired databases through spatial analysis equations (regression function applied on the digital elevation model to interpolate climatic values) using specific techniques (vector-raster conversion, classification), according to the purpose and the necessary database type.

The spatial analysis stage has at its centre the integration of the submodels’ results as a general model based on the general equation of spatial analysis, which defines the conceptual model, with the main purpose of identifying the surfaces which are favourable to dwarf mountain pines. Additionally, in this stage, special attention was given to the NDVI estimation and the supervised classification of this indicator, with the main purpose being the indirect validation of the model, as well as the identification of the spatio-temporal dynamics of the analysed species, in case the final favourability result is valid.

The use of the proposed equation to integrate all the variables of favourability into the spatial analysis equation based on the ArcGIS geoinformation software highlights the the great influence of the two factors (the average annual temperature and the amount of rainfall), as evidenced by the development of Pinus Mugo in these sectors.

The validation stage is performed on two main directions: the direct validation, based on samples from the field and the indirect validation, based on the areas obtained by deriving the NDVI index for the year 2018, considered a characteristic year when the areas with Pinus murgo were identified in the field. The validation of the present model is highly important because, based on the results of the validation, one can determine the final results, the identification of the changes occurring in this territory over the seven years for which the NDVI values were derived. Only if the validation is a positive one can a final result be obtained, otherwise one must return to the model’s structure and change it, as well as the quality and the consistency of the databases.

3. Results

3.1. Spatio-Temporal Variability Assessment of Pinus Mugo Using Landsat Imagines

The monitoring of forested areas from Romania, and especially of the changes occurring due to human intervention on mountain territory (tree cutting, pests, forestry-cultural works), has been performed indirectly by using Landsat OLI8 satellite images [49,50]. Two main categories of interest have hence been identified: forest gain and forest loss. However, the research was limited to areas of approximately five hectares where disturbances occurred.

On larger territories, this process requires considerable effort (a large number of specialists, numerous resources, reduced precision), which imposes the use of indirect monitoring methods of the areas with tree vegetation, by using geospatial techniques, offering an overall view and highlighting, at the same time, the vulnerable areas.

The spatial analysis from the present study was performed based on the derivation of the NDVI index for two years which were considered representative, 1987 and 2018. High-quality LANDSAT 5, 7 and 8 images were identified for these years, without cloud coverage in the analysed area.

The NDVI index derived for 1987 highlights the large extension of the areas covered with vegetation, approximately 84% of the whole territory, while almost the same percentage (86%) was calculated for 2018. The homogeneity of the alpine area due to the associations of species makes the differentiation between grass and wood species a difficult task, a problem which also affects the identification of the studied species.

For the study area, the spectral response of the bands being integrated in the equation of deriving the NDVI index is represented by numerical values between 0.98 and 0.60 for vegetation. For the year 1987, the value is 0.92, while for 2018, the value is 0.51, a result which highlights the reduction of the compact vegetation surface for the studied habitat (Figure 4). The analysis of the spatial distribution of the vegetation in the study area indicates a very heterogenous distribution of vegetation classes for both years, due to the intercalation of areas without vegetation (bare lands) among those with grass vegetation and with dwarf mountain pines (Pinus mugo shrub) from the upper forest tree line, section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation as well as the experimental conclusions that can be drawn.

The analysis of the wooden component, with the purpose of highlighting the dwarf mountain pine associations, faces great difficulties due to the fact that most of the study area is covered with different coniferous associations (spruce, juniper, Swiss pine) which have almost the same spectral response as the dwarf mountain pine.

In order to highlight only the areas classified as being covered with Pinus mugo, one used the technique of reclassifying the raster databases to integer values (value 1 was attributed to all the areas with Pinus mugo and value 2 for the rest of the areas where other types of vegetation were identified), which were used in the analysis process of the temporal changes affecting the studied species in this territory.

The analysis of the final results indicates the spread of dwarf mountain pines on 74% of the study area in 1987 and on 63% in 2018, the large surface being motivated mainly by the topography, the favourable altitude as well as by the favourable climatic and pedologic conditions.

A reduction in the areas covered with Pinus mugo is also visible in the time interval under analysis, which is motivated by some studies through massive cutting of dwarf mountain pines for the expansion of alpine pastures and the destructuring of the Pinus mugo habitats by extensive harvesting of its buds for alternative medicine.

The analysis of the changes in vegetation types indicates (Figure 5 and Figure 6) firstly the areas without modifications in the structure of the analysed vegetation (approximately 57% of the entire study area), a fact which is explained by the large percentage of areas without vegetation and short grass vegetation, which is characteristic to Carpathian areas above 2000 m. The spatial extension of Pinus mugo is reduced (approximately 15% of the total area of Piule-Iorgovanu Massif), which is explained by the migration of the dwarf mountain pines towards the areas with loose vegetation and their growth on the debris close to rocks, as well as in the areas with windthrow where the dwarf mountain pine settled as an invasive species. The ratio of expansion-reduction in surface covered by Pinus mugo is dominated by reductions (with approximately 27%), which is very well correlated with the results of the spatial analysis based on the supervised classification, which validates the fact that there has been a reduction in the extension of the Carpathian dwarf mountain pines between the two representative years.

3.2. Identification of Favourability for Pinus Mugo in The Romanian Carpathian Area

Among the environmental factors, the climate is highly important through its influence on forest vegetation. The most important climatic factors are light, heat, humidity and wind [51]. The dependency between the environmental factors and the requirements of Pinus mugo was highlighted by the resulting cartographic database which classifies the territory in favourability classes according to 14 main ecological factors under analysis (Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19 and Figure 20).

The variations in air temperature in Romania between the warm and cold seasons influence the vegetation period of tree species. Thus, the vegetation periods alternate with resting periods during winter, the variations in temperature being influenced by latitude, altitude, slope exposure, etc.

The average annual temperature in the study area ranges between −2.3 and 10.4 °C and manifests its influence by inducing a very low favourability for Pinus mugo in 50% of the Oriental Carpathian territory, 26.3% of the Occidental Carpathians and, due to their higher altitude, 17.1% of the Meridional Carpathians (Figure 7). The most extended territories where this influence is present occupy 2870 km² in Făgăraș, Parâng and Retezat-Godeanu mountain divisions from the Meridional Carpathians, 1267 km² in the Oriental Carpathians at high altitude and smaller areas, of 221 km², in the Apuseni Mountains from the Occidental Carpathians. Thus, Pinus mugo (oligotherm species), must adapt its rhythm of biological activity to the variations in thermal regime from the Carpathian area, as it is necessary to resist to frost in the cold season and to massive snow fall.

The influence of humidity on dwarf mountain pine development manifests itself on two directions, with favourable, but also restrictive effects. From the category of negative effects induced by a large amount of precipitation, one notices the high soil erosion on already skeletal soils and the influence on the pollinating period (which is short). As a positive aspect, however, high humidity allows this species to form pure vegetation associations on large surfaces, creating areas with a specific microclimate.

In the study area there is a variation of the average annual precipitation amount between 557 and 1620 mm/year. In the precipitation range of 901–1400 mm/year, the favourability for Pinus mugo is low and in the areas where the precipitation amount is below 900 mm/year and above 1400 mm/year, the favourability is very low (Table 2).

The Oriental Carpathians, through their specific average annual precipitation, offer acceptable growth conditions for Pinus mugo on 30% of their territory (19,939 km²), being followed by the Meridional Carpathians, on 13.3% of the analysed area (8810 km²) (Figure 9). The least represented division of the Carpathians in what concerns the favourability for Pinus mugo induced by the average annual precipitation amount is the Occidental Carpathians. Here, the lower altitude and the effect of orographic barriers in the way of the western air masses that bring precipitation determine conditions for low favourability on 20% of the territory (13,261 km²) and restrictive conditions on 6.7% (4450 km²) (Table 3).

The length of the bioactive time period (Figure 8) is determined by the number of days with temperature above 10 °C. This factor varies in the study area between 0.6 and 6 months, the longest duration being specific to the area of low mountains and depressions inside the mountain range. The favourability for Pinus mugo is positively influenced on 23.1% of the analysed territory where the bioactive time period is between 1 and 4 months/year, which helps the development of the dwarf mountain pine, while on the remaining 76.9% of the territory (51,035 km²) this indicator limits the spatial and temporal development of the species (Table 3).

The relative atmospheric humidity in July (%) is important mainly because, during this month, there is generally the lowest amount of precipitation. This indicator varies in the Romanian Carpathians between 66% and 91% (Figure 10), the most favourable territories for Pinus mugo being characterised by values above 75%, mostly found in the Oriental Carpathians (on 77.4% of their territory) and on 23.9% of the Meridional Carpathians (Table 3). For this indicator, 66.8% of the mapped territory is included in the medium favourability class, being the only one out of the 14 analysed factors which reaches this level of favourability.

In the Carpathian mountain range from Romania, the wind regime ranges between 3.9–36.7 km/h (Table 2). The negative influence of the wind regime on dwarf mountain pines is felt on 66,071 km² (which represents 99.5% of the whole study area), according to the analysis performed (Figure 11). Thus, according to the Beauford classification, the wind speed in the range of 3.9–5 km/h is correspondent to a light breath of wind, 5.1–11 km/h corresponds to a light breeze, 11.1–19 km/h corresponds to little wind, 19.1–28 km/h corresponds to moderate wind, while wind speeds between 28.1 and 37 km/h correspond to strong winds. Although the wind speed is a limitative factor in the development of the dwarf mountain pine, this species’ resilience to wind is due to its growth pattern (high elasticity of the stems, crawling and ascending shrub aspect) and its plastic rooting system.

Altitude is the determining ecologic factor in the development of the dwarf mountain pine due to its effect on the previously presented climate factors (temperature, average annual precipitation, wind speed, atmospheric humidity) as well as on the soil characteristics which work together in supplying the necessary factors for the development of Pinus mugo. Specialised literature mentions the altitude of 1350–1450 m as a minimum for the growth of the dwarf mountain pine in the Western Carpathians and a maximum of 1850 m in the Tatra Mountains, and even above 2100 m in the southern Romanian Carpathians (Retezat, Parâng, Bucegi) [16,52].

The altitude in the study area is in the range of 56–2544 m (Figure 12). The mountain areas with lower altitudes than 1300 m offer a very low favourability for the growth of Pinus mugo, while those above 1300 m offer better conditions (Table 2). Hence, the Meridional Carpathian range ensures good development conditions on 4439 km² (6.7% of the whole Carpathian range), while the most restrictive Carpathian unit (30,939 km²) is represented by the Occidental Carpathians, characterised by the lowest altitude in the whole Carpathian range (its maximum being only of 1849 m) (Table 3).

The slope aspect represents a morphometric parameter which influences indirectly the favourability and the restrictiveness for dwarf mountain pines. The mountain surfaces with slopes exposed towards the south or the southeast receive a higher caloric energy than the slopes facing north and northwest. The northern and northwestern slopes are also more exposed to the air masses which are rich in precipitation, followed by the southern and southeastern slopes.

From this point of view, 92.1% (61,160 km²) of the Carpathian territory offers good conditions for the development of Pinus mugo on quasi-horizontal surfaces, which receive a higher quantity of caloric energy and determine for 7.9% (5220 km²) of the analysed area a medium favourability.

In what concerns the influence of soil on the development of forest vegetation, it is known that trees and bushes set their roots in the soil and absorb through them the mineral salts needed for their development. Thus, the thicker soils offer the possibility of growing deeper roots. In the case of shallow soils which are specific to the mountain area above the forest boundary, it is important to know the physical and chemical characteristics of soils as the trees need to adapt their roots and are more vulnerable to windthrow and massive snowfall.

The saturation of the soil with bases in the Romanian Carpathians varies between 0.002 and 82%, the areas with larger values than 45% (624 km² of the analysed territory) offer a very low favourability, while the areas with values below 45% offer a better favourability (Figure 14, Table 3).

Soil acidity, determined as the value of water pH, varies in the study area in the range of 4.3–8.2, the most favourable areas for the growth and development of dwarf mountain pines are those characterised by values between 5.4 and 3.6, which can be found on 16.3% of the territory (10,801 km²). The Oriental Carpathians are the best represented mountain range from this point of view, on 8.2% of their territory (5895 km²) (Figure 15, Table 3).

The azote from soil was determined procedurally using the Kjeldahl procedure and the punctually value was interpolated using GIS technologies. In the Romanian Carpathians ranges between the values of 1.06 and 7.11 (Figure 15). As in this case we are interested in the influence of this indicator on Pinus mugo, the focus is laid on the territories with values above 4, which represent 97% of the study area, mostly extended in the Oriental Carpathians (on 32,632 km²) (Table 3). A reduced quantity of the azote from soil will lead to the yellowing of the needles and a decrease in the plant’s growth. That is why an analysis of the azote trophicity is also required as wooden plants cannot take in atmospheric azote directly, thus, a thick layer of humus is required.

Edaphic volume (Ve, fractions of unity), particle size distribution, as well as waterstable macroaggregates (structural macrohydrostability) have been determined for all monitoring plots at national level and this values was interolated using GIS technologies.

The edaphic volume (m³/m²) determines a very low favourability for Pinus mugo when the values are higher than 0.45. In the study area, this indicator ranges between 0.17–1.67 (Table 1). The best values of the edaphic volume for the development of dwarf mountain pines are found on 5956 km² from the Oriental Carpathians (Figure 17,

Table 4. Favourability classes for Pinus mugo.

Favourability Classes	Eastern Carpathian	Southern Carpathian	Western Carpathian	Area (km2)	Factors
				%
Very Low	33,206	11,325	17,491	62,023	Average annual temperature
	50.00	17.10	26.30	93.40
Low	1267	2870	221	4357
	1.90	4.30	0.3	6.60
Very Low	14,534	5385	13,261	33,180	Average annual precipitation
	21.90	8.10	20.00	50.00
Low	19,939	8810	4450	33,200
	30.00	13.30	6.70	50.00
Very Low	26,669	8154	16,212	51,035	Bioactive period
	40.20	12.30	24.40	76.90
Low	7804	6041	1500	15,345
	11.80	9.10	2.30	23.10
Low	2402	6264	13,345	22,011	Relative atmospheric humidity
	3.60	9.40	20.1	33.20
Medium	51,381	15,885	14,285	81,551
	77.40	23.90	21.50	66.80
Low	34,348	14,040	17,683	66,071	Wind regime
	51.70	21.20	26.60	99.50
Very Low	75	135	0	210
	0.10	0.20	0.00	0.30
Very Low	30,939	9756	17,001	57697	Altitude
	46.6	14.7	25.6	86.90
Low	3534	4439	710	8683
	5.30	6.70	1.10	13.1
Very Low	2891	1089	1239	5220	Aspect (Slope exposure)
	4.40	1.60	1.90	7.90
Low	31,582	13,106	16472	61,160
	47.60	19.70	24.80	92.10
Very Low	34,147	14,131	17378	65,656	Degree of soil saturation with bases
	51.40	21.30	26.20	98.90
Low	276	44	305	624
	0.40	0.1	0.50	0.90
Very Low	28,527	11,287	15,666	55,480	Soil acidity
	43.00	17.0	23.60	83.60
Low	5895	2888	2017	10,801
	8.90	4.40	3.00	16.30
Very Low	1790	0	127	1917	Azote from soil
	2.70	0.00	0.20	2.90
Low	32,632	14,175	17,556	64,364
	49.20	21.40	26.40	97.00
Very Low	28,467	13,970	17,377	59,814	Edaphic volume
	42.90	21.00	26.20	90.10
Low	5956	205	306	6467
	9.00	0.30	0.50	9.70
Very Low	3189	1010	1743	5942	Usable soil water
	4.80	1.50	2.60	9.00
Low	31,233	13,166	15,940	60,339
	47.10	19.80	24.00	90.90
Very Low	0	42	576	618	Soil compactness
	0.00	0.10	0.90	0.90
Low	34,425	14,133	17,109	65,667
	51.90	21.30	25.80	98.90
Very Low	3306	857	2191	6354	Soil texture
	5.00	1.30	3.30	9.60
Low	31,110	13,319	15,494	59,922
	46.90	20.10	23.30	90.30

).

The usable water from soil (%) ranges between 1.8 and 20.9% for the study area, but the best conditions for Pinus mugo development are found in the range of 10.1–20% (for 90.9% of the study area, the largest territory included in this category being located in the Oriental Carpathians −31,233 km²) (Figure 18).

Soil compactness represents an important pedologic factor for the development of trees and bushes, thus, loose soils which are rich in humus and have a low acidity are favourable to plant growth, while on compact soils, the favourability is very low (this being the case of the dwarf mountain pine as well). According to pedologic studies, in the study area, surfaces with a limitative compacteness degree were identified on only 618 km² (Figure 19), these territories being located in the Banat Mountains and in the Retezat Mountain unit from the Meridional Carpathians.

Soil texture conditions the regime of soil humidity, soil permeability and consistency, as well as the value of the edaphic volume. As a consequence, it has an important role in the growth and development of the trees and bushes from the study area. The soils from the Romanian Carpathians offer favourable conditions for the development of the dwarf mountain pines on 90.3% of the whole study area, where there are soils with sandy, sandy-loamy, loamy-sandy and loamy texture. The rest of the territory is characterised by loamy-clayey, clayey-loamy and clayey texture on 6354 km² which imposes restrictive conditions for the development of Pinus mugo (the most extended surfaces with these types of textures being found in the Oriental and Occidental Carpathians (Figure 19).

In order to identify the favourability classes for the development of Pinus mugo in the Romanian Carpathians and the vegetal associations in which it is included, a digital database was used, including ecological, climatic, pedologic and morphometric factors which determine through their cumulated influence various degrees of favourability or restrictivity for the analysed species. The final map of favourability for Pinus mugo was created using the techniques of GIS spatial analysis and the Raster Calculator extension, using the formula:

$$\mathrm{P}\mathrm{M}\mathrm{F}\mathrm{a}\mathrm{v}=\sqrt{\frac{\mathrm{T}\mathrm{m}\mathrm{a}\ast \mathrm{P}\mathrm{P}\mathrm{a}\mathrm{v}}{4}\ast \left({\mathrm{Q}}_{\mathrm{W}}+{\mathrm{Q}}_{\mathrm{H}}+{\mathrm{Q}}_{\mathrm{T}10}+{\mathrm{Q}}_{\mathrm{A}}+{\mathrm{Q}}_{\mathrm{B}}+{\mathrm{Q}}_{\mathrm{V}}+{\mathrm{Q}}_{\mathrm{A}\mathrm{C}}+{\mathrm{Q}}_{\mathrm{U}}+{\mathrm{Q}}_{\mathrm{C}}+{\mathrm{Q}}_{\mathrm{N}}+{\mathrm{Q}}_{\mathrm{S}}\right),}$$

(1)

where: PM Fav—Favourabily for Pinus mugo, Tma—average annual temperature (°C); PPav—average annual precipitation (mm/an); Q_W—quantitative index for wind; Q_H—quantitative index for relative atmospheric humidity in July; Q_T10—quantitative index for temperature above 10 °C; Q_A—quantitative index for Aspect; Q_B—quantitative index for degree of soil saturation with bases; Q_V—quantitative index for edaphic volume; Q_AC—quantitative index for soil acidity; Q_U—quantitative index for usable soil water; Q_C-quantitative index for soil compactness; Q_N—quantitative index for azote grid; Q_S—quantitative index for soil texture.

By analysing the final result, classified according to the favourability classes for Pinus mugo, one identifies three classes in the Romanian Carpathians (low, medium and high) (Figure 21).

In the Oriental Carpathians there is a high favourability for Pinus mugo on 266.1 km² (0.77% of the mountains’ surface) and medium favourability on 949.5 km² (2.76%), while in the Meridional Carpathians the high favourability class is found on 525.5 km² (3.71%) and the medium one on 1957.3 km² (13.81% of their surface). The most reduced values which offer good development conditions for Pinus mugo are found in the Occidental Carpathians (only 40 km² with high favourability and only 180 km² with medium favourability).

4. Discussion

The identification of the dwarf mountain pine shrubs was performed using the method of supervised classification [53,54,55,56], by classifying the result of deriving the NDVI index for the study area. In order to implement the supervised classification in the geoinformatic software, two different digital databases were required: vector databases represented by samples identified in the territory as areas with Pinus mugo (575 samples and a total area of 4.82 km²), samples of areas with rocks without vegetation (125 samples and an area of 0.07 km²), samples with forest vegetation, different from dwarf mountain pine (43 samples and an area of 12.53 km²), samples representing grass vegetation (36 samples and an area of approximately 3.03 km²) and raster databases represented by the NDVI index, derived for the two representative years, 1987 and 2018. After performing the supervised classification with the previously mentioned databases, several raster databases were obtained which highlight by extrapolation the areas with the three types of samples used in the analysis.

The validation of the spatial analysis based on the supervised classification was performed by comparing the final result with the samples from the field where Pinus mugo was present. Ten samples representing 0.0015 km² were compared, others than the samples used as input data for the supervised classification, the 92.3% overlap indicating a very good validation of the classification process.

The monitoring of the temporal changes (between the representative years) of the areas with dwarf mountain pines was performed by implementing equations of spatial analysis in geoinformational software, the result being a reduction of the two databases with the attributes 1 and 2. The final result is represented by a digital database with integer attribute values: −1 (areas characterised by an expansion of the dwarf mountain pine between 1987 and 2018), 0 (areas with no changes in the types of vegetation from the study area) and 1 (areas with a reduction from the year 1987 in the surface of the territory with Pinus mugo).

As a first method of validating the model, the areas with a higher frequency than 25% of Pinus mugo were digitised using the Atlas of Forest Tree Species [57], where the data referring to the frequency of Pinus mugo in the field observations from Romania correspond to national forestry inventories. Hence, out of the 19 areas reported as having a frequency of the dwarf mountain pine higher than 25% of the total vegetation, 6 areas are located in the proximity of the high favourability for Pinus mugo according to the GIS model of analysis developed in the present study.

The second method of validation was based on the identification of the overlaying degree for the 585 polygones representing the vegetal associations with Pinus mugo with the theoretical favourability classes identified by applying the model (Figure 22). After this validation, for an area of 4.23 km² with Pinus mugo out of the 4.84 km² identified for the test sector there is a very well overlay with the high favourability class. In addition to this, there is an overlay of 0.26 km² with the medium favourability class for which the validation rate is of 92.82%. This validation rate indicates a high predictability degree for the quantitative model applied.

In order to clearly show the dependence of the final result according to each factor we made a transverse profile on the NE-SV directives (Figure 23), thus exemplifying the factor’s dependence on the favourability classes that the profile intersects. This makes it clear that most factors have the same variation in getting the final result. Analysing the figure we can see that both the rainfall and the average annual temperature have a small favourability class for small final favourability at Pinus mugo (visible on the cross-sectional profile between 3400 and 4200 m), and for the precipitations in the range 8100–9600.

Additionally, in the same trend is the slope orientation factor that integrates medium favourability for the same low-favour surface for Pinus mugo. From the point of view of the orientation, obviously its influence in the final result is evident, so we can identify favourable values from medium to very low in sectors intersected by the transversal profile on which the final average values for Pinus mugo are recorded. The decrease of orientation as influence in terms of final result is given by soil texture, soil compactness, useful water capacity, factors that make the slope orientation as a factor that largely influence the quantity of useful water in the soil to be significantly diminished identified by profile on the 1400–2300 m sector and the 5500–5700 sector).

Noteworthy is that the only parameter that does not correlate with the final result is the relative atmospheric humidity in July, a factor that maintains the highest maximum value for all the Pinus Mugo favourability classes intersected by the cross-sectional profile.

On the opposite side, it is found that all the factors taken into account fall within the low favourability class for all sectors intersected by the cross-sectional profile. Positive diminution or compensation of their influx in the final result, that of identifying areas of favourability to Pinus mugo, is due to the significant variation of the other factors taken in the analysis and exemplified earlier.

5. Conclusions

The purpose of the study was the classification of the area with the Pinus mugo species, found in various vegetation associations, according to its favourability, and the identification of the dynamics of the spatial extension of Pinus mugo for the 1987–2018 time period, using the technique of GIS spatial analysis and LANDSAT satellite images.

The implementation of the spatial analysis model based on databases, expert knowledge and modern GIS techniques of investigation has led to a highly accurate identification of the spatial dynamics of the vegetal associations with Pinus mugo, highlighting the fact that this spatial analysis can be successfully implemented and can support environmental practices based on postdiction and prediction, having as the main purpose a sustainable management of the forestry resource, depending on an integrated analysis of the area/a zonation of the final result.

A quantitative GIS model was firstly implemented in order to classify the study area in favourability and restrictivity classes for the Pinus mugo species, based on the transfer of a pre-existing quantitative model designed as a table to a GIS spatial model, which enabled the correlation of the factors included in the model and the generation of analysis data for this correlation between ecological factors. This made possible the extension of the study area at the level of the whole Romanian Carpathian mountain range by which an overview was produced of the areas where this species requires special conservation measures.

Considering that the literature does not frequently deal with this subject and there are no models of spatial analyses derived for the identification of Pinus Mugo species or other forest shrubs, we have developed our own model Landsat images based on the fact that Landsat’s images are free to use and a good resolution for Pinus Mugo’s spatial distribution specificity in Romania leaving the possibility to scale-up the model on test areas by using higher resolution images and not based on the UAV images.

From the analysis of the spatial and temporal dynamics of Pinus mugo on Piule- Iorgovanu massif (selected as a test area), no changes in the structure of vegetal associations were identified for 57% of the area, for the previously mentioned time period. The extension of dwarf mountain pines was approximated for 15% of the test area, as these plant associations migrate towards the areas normally covered with rare forests at the limit of the alpine zone, where they behave as an invasive species, and on debris and rocks.

For the analysed time period, the reduction of the area is visible for 27% of the territory, this reduction being caused mainly by anthropic activities which had various purposes: the extension of alpine pastures, the use of dwarf mountain pines as fire wood, the harvesting of this species for the pharmaceutic and cosmetic industries, all these being correlated with a low rate of growth.

The results are useful in the process of territorial management for the alpine areas. This management is implemented by local public authorities with the purpose of valorising the ecologic and economical functions of alpine zones and of identifying specific conservation measures, considering the fact that most of the areas where Pinus mugo grows are integrated in the protected areas of the NATURA 2000 sites, requiring special attention.

By the way in which the spatial analysis model is structured (in clearly presented stages), the type and structure of the input databases (databases acquired freely, structures of common databases—vector, raster, numerical), the model presented in this study can be repeated and implemented successfully on areas with similar environmental components.