Identifying the Spatiotemporal Transitions and Future Development of a Grazed Mediterranean Landscape of South Greece

Abstract

Spatiotemporal changes over previous decades in grazed Mediterranean landscapes have taken the form of woody plant encroachment in open areas (e.g., grasslands, open shrublands, silvopastoral areas), altering its structure and diversity. Demographic and socioeconomic changes have played a significant role in landscape transformations, mainly by causing the abandonment of traditional management practices such as pastoral activities, wood harvesting, and agricultural practices in marginal lands. This study aimed to quantify and evaluate the spatiotemporal changes in a typical grazed Mediterranean landscape of Mount Zireia during 1945–2020, and to investigate the effect of these changes on the future development (2020–2040) of land use/land cover (LULC) types. Cartographic materials such as aerial orthophotos from 1945, land use maps of 1960, Corine Land Cover of 2018, and recent satellite images were processed with ArcGIS software. To estimate the future projection trends of LULC types, logistic regression analyses were considered in the framework of CLUE modeling. The results indicated that the strongest trend of spatiotemporal changes were forest expansion in open areas, and grasslands reduction, suggesting that the LULC types that were mainly affected were forest, grasslands, and silvopastoral areas. Future development prediction showed that forests will most probably continue to expand over grassland and silvopastoral areas, holding a high dynamic of expansion into abandoned areas. The reduction in grasslands and silvopastoral areas, independent of environment and biodiversity implications, represents a major threat to sustainable livestock husbandry based on natural grazing resources.

1. Introduction

Mediterranean landscapes are considered highly diverse areas in terms of history, geography and land uses. Several civilizations from ancient times have left a rich cultural heritage promoting this variety [1]. The Mediterranean landscapes, as a result of their long history of human activities, with a unique combination of topographic and climatic variability, have generated a rare combination of unique, but fragile, diverse species-rich ecosystems [2,3]. The Mediterranean basin is the second largest biodiversity hotspot in the world, holding more than 25,000 plant species [4]. The long history of human intervention in this area has formed plant communities that are considered as “man-made” and composed of natural components, a fact that has a significant value in setting goals and methodology for sound conservation interventions [5]. The last 75 years of technological advances, such as the introduction of heavy machinery in farming activities [6], trade globalization, the creation of the European Economic Community [2], and the Common Agricultural Policy (CAP) [7,8], have driven dramatic changes in these ecosystems unlike those experienced in the past [9]. In recent years, climate changes, along with unbalanced land use activities (e.g., coastalization, undergrazing, and land abandonment), have facilitated Mediterranean ecosystems change [10]. Two opposite trends of landscape evolution have occurred in the Mediterranean region in recent decades. Forest cover increased around the northern edge of the Mediterranean region (south European countries) and decreased around its southern edge (mainly in the Maghreb countries). This increase in forest cover in northern Mediterranean landscapes is mainly attributed to the abandonment of marginal agricultural lands [11,12], while the decrease in forest in the south is attributed to the expansion of cropland in marginal areas initially dominated by woodlands [13]. The above changes have followed the socioeconomic trends of land abandonment in rural areas in the north versus the increased population pressure in rural areas in the south [9,14].

One of the main land use activities in Mediterranean landscapes is pastoral activities [2,9,13,15]. Approximately one-fifth of European agricultural lands are dedicated to extensive livestock grazing, with the majority being situated in southern Mediterranean Europe, including the Balkans. Furthermore, 80% of Europe’s sheep and goat flocks are located in Spain, Italy, Greece, and southern France [16]. Grazing is considered a major landscape-changing factor directly related to human activities, especially in Mediterranean areas [17,18]. Greek landscapes have historically been grazed by livestock in quite a similar way as modern practices, and are highly influenced by the changes in traditional pastoral activities [19,20]. Recently, significant changes emerged in the traditional extensive livestock production systems of Greece, mainly related to the reduction in the number of local and transhuman flocks of free-grazing animals (sheep, goats, and cattle) [8,20,21,22,23]. These changes follow the land abandonment trend already mentioned for the European part of the Mediterranean region, and they highly contribute to the spatiotemporal transitions occurring in grazed areas. These transitions are taking the form of woody plant expansion in open areas, transforming grasslands, open shrublands, silvopastoral areas and abandoned agricultural areas, into forest or dense shrublands [21,24,25,26,27,28].

The study of land use/land cover (LULC) change provides an important aspect in understanding the history of spatiotemporal transition patterns, derived from landscape changes. Spatiotemporal transition patterns produce useful data for studying the effect of physical and socioeconomic interactions, land use conflicts, and influences on landscape changes [29,30,31]. Analysis of spatiotemporal changes and transitions is typically conducted within the geographic information systems (GIS) environment [30,32], with visual photointerpretation of a time series set of aerial photographs [33,34] through digital processing of multispectral satellite images [32], or more recently through object-based recognition technics [35]. The development of transition matrices has become an important part of landscape history analysis [36]. New tools and indicators of LULC changes derived from the matrices have emerged, addressing issues related to the annual rate of changes [30], persistence and net changes as quantity difference and swap as allocation difference [37], and identifying systematic or main transitions [29,36].

Predicting the future development of LULC types and transitions is an effective and reliable technique for evaluating both the causes and the significance of past and present conditions, usually under future scenarios [27,38,39]. Several spatiotemporal models for LULC future projections have been proposed over the years [31], including the adoption of empirical models for LULC prediction such as logistic regression approaches (e.g., CLUE modeling framework, LCM, MaxEnt) [39,40,41]. The use of regression analysis in landscape prediction studies contributes to understanding and describing the change mechanisms and processes of LULC types, provides an advanced statistical environment for analyzing multivariate components, and finally, predicts the LULC changes [14,42]. The above prediction models can also produce accurate results to support policy makers, land managers, and scientists in reaching sustainable landscape management decisions [41].

Spatiotemporal changes have a significant effect on altering landscape structure in terms of landscape composition and configuration [43]. These changes can be easily evaluated with the use of landscape metrics [44,45], applied in spatiotemporal studies of landscape changes [22,24,46,47], or in the future projections of land use changes [48].

Overall, there is a limited amount of published information regarding the spatiotemporal changes in grazed landscapes, especially for the eastern part of the Mediterranean region, and particularly about the influences of land abandonment in the future development of land uses that are related to pastoral activities. Therefore, the present research aims to: (a) quantify and evaluate the spatiotemporal changes of a typical, grazed Mediterranean landscape of south Greece (Mt Zireia landscape), (b) investigate the effect of these changes on the future development of the most significant LULC types, and (c) identify their correlation to a set of landscape driving factors. Finally, the overall effect and interactions of socioeconomic changes are explored, focusing on pastoral activities in LULC transitions and future development.

2. Study Area

Mount (Mt) Zireia (or Kyllini), located in the Peloponnese peninsula (South Greece), was selected for the study. Mount Zireia is the second highest mountain in the Peloponnese, located in the Korinthos prefecture 115 km west of Athens (Figure 1).

The study area covers 39,762 ha of land inhabiting 3777 people living in 19 village communities–municipalities subdistricts. Elevations in the study range from 310 m to 2374 m a.s.l. A large gorge, called Flampouritsa, divides the mountain into two areas, “Mikri” (small) Zireia and “Megali” (big) Zireia. Mt Zireia, apart from the highest point of 2374 m, has other seven peaks above 2000 m (four in Megali and three in Mikri Zireia). The multiple ridges created by the mountain tops, in combination with valleys and plateaus, create a particularly diverse relief of hills, plains, cliffs, and canyons. More than two-thirds of the study area is part of the network of Natura 2000 protected areas (pSCI, SCI or SAC, SPA) [49]. Two main hydrological basins are found in the area, creating the natural lake Stymfalia (area 15,285 ha) to the south, and the artificial lake Doxa (area 48 ha) in the west (Figure 1), which greatly affect the microclimatic conditions and facilitate the touristic development of the area. Lake Stymfalia is closely connected to Greek mythology, and especially with Heracles’ legendary labors. According to mythology, the lake was full of aggressive man-eating Stymphalian birds, and Heracles’ sixth labor was to exterminate them [50].

The climate, according to Köppen–Geiger climate classification, is a hot summer Mediterranean climate (coded as “Csa”) [51]. The mean annual precipitation has varied over the last 60 years, from 418.62 mm (in 1993) to 1056 mm (in 2005), while the mean annual temperature varied from 12.59 °C (coldest year in 1976) to 15.55 °C (warmest year in 2010) [52].

The main land uses of the area are forests, rangelands, and agricultural areas. Rangelands include grasslands, shrublands, and silvopastoral areas with less than 40% tree cover and grazed by sheep and goats. Agricultural areas are cultivated mainly with annual crops such as beans, corn, barley, and wheat [52].

According to the official census report derived from the Hellenic Statistical Authority [53], the temporal evolution of socioeconomic data from 1961 (oldest available data) until the most recent available data of 2011, showed that in the last 50 years, the total population, active workforce, and employees in the primary economic sector has rapidly been reducing (

Table 1. Temporal evolution (1961–2011) of the total population, active workforce, and primary sector employees of the nineteen village communities in the study area.

	Total Population	Total Working Force	Employees in the Primary Sector	% of Primary Sector Employees per Total Working Force
1961	7420	3632	3169	87.25
2011	3777	1354	682	50.37

(Source: Hellenic Statistical Authority.)

), following the general trend of land abandonment that many researchers have reported for the Mediterranean region [8,44,54,55]. Age structure analysis indicates that the human population is becoming older. Indeed, 62% of the population was under 44 in 1961, versus 47% in 2011 (

Table 2. Temporal evolution (1961–2011) of age structure (as percentages) of the local population of the nineteen village communities in the study area.

	1961		2011
	% (0–44)	% (45–)	% (0–44)	% (45–)
Total	62.34	37.66	46.93	53.07

(Source: Hellenic Statistical Authority.)

).

In contrast, the local population over 45 years old increased from 38% to 53%, for the same period. The above data are in line with demonstrated demographic change in the Mediterranean region and the movement of the mainly younger population from rural areas to urban centers [8,11,44].

Census data from the Hellenic Statistical Authority and the Payment and Control Agency for Guidance and Guarantee Community Aid [56], regarding the historical data of transhumans [57], revealed that the number of grazing animals (mainly sheep and goats) and their farms have significantly reduced in the last 50 years (

Table 3. Temporal evolution (1961–2020) of grazing animals (sheep and goats) in the study area.

	Number of Animals		Percentage of Change
	1961	2020
Sedentary	28,750	27,595	−4.02
Transhumant	38,230	13,717	−64.12
Total	66,980	41,312	−38.32

) [20].

The total number of grazing animals decreased from 1961 to 2011 by 38% (Table 3). This reduction was more intensive for transhuman animals (more than 64%) and less for sedentary animals (almost 4%). According to the available inventory data, the number of sedentary animal farms significantly reduced by 80% during a similar period (Hellenic Statistical Authority, 1961 to 2000). This reduction follows the similar trend of change as the number of people that are employed in the primary sector of the economy (Table 1).

3. Materials and Methods

3.1. Land Use/Land Cover Changes. Spatiotemporal Transitions

The following cartographic materials (Figure 2) were considered: (a) digital aerial orthophotographs of 1945 with a spatial resolution of 1 m obtained from the National Cadastre of Greece (georeferenced to the Hellenic Geodetic Reference System 1987-HGRS87); (b) satellite images obtained from the Google Earth Pro program for the years 2017, 2019 and 2020 (georeferenced to HGRS87); (c) maps of forest vegetation and land cover for 1960 (scale 1:20,000), obtained from the Ministry of Agriculture in digital format (shapefile in HGRS87); and, (d) digital maps of Corine Land Cover 2018 (shapefile reprojected in HGRS87).

Aerial orthophotographs from 1945, as well as the recent Google Earth satellite pictures, were digitally processed using the software ArcGIS v.10.8.1, to produce LULC maps for 1945 and 2020. To identify the distinct LULC types, on-screen visual photointerpretation and manual delineation of LULC polygons in shapefile format were performed within the ArcGIS environment (Figure 2). The chosen analysis used a classification scheme consisting of eight categories of LULC types and was based on the Greek Forest Service’s LULC classification system (

Table 4. Classification scheme of land use/land cover types in the study area.

Land Use/Land Cover Types	Description	Codes
Agricultural areas	Areas covered with annual or permanent crops	AG
Grasslands	Areas dominated by herbaceous plants, with woody vegetation cover of less than 10%	GR
Open shrublands	Areas dominated by sparse woody shrubs with less than 40% cover	OS
Dense shrublands	Areas dominated by dense woody shrubs higher than 40% cover	DS
Silvopastoral areas	Open grazed forest with tree cover between 10 and 40%	SI
Forest	Forest areas with tree cover higher than 40%	FO
Barren areas	Mainly bare lands with little or no vegetation	BA
Urban areas	Areas with man-made features, mainly villages	UR
Lakes	Large area of water surrounded by land	LA

). According to the chosen classification system, numerous elements on aerial orthophotos and Google Earth images were recognized by using common photographic keys (tone, texture, pattern, shade, form, and size) and feature association [15,21,29,33,34]. Special attention was placed on identifying tree and shrub cover density patterns with the use of crown density scales [58]. The 1960 forest vegetation and land cover maps in shapefile format were a valuable resource for the 1945 LULC mapping, since they served as a reference map and guided the photointerpretation. The minimum mapping unit of the reference map was one hectare, and the same unit was chosen for the 1945 and 2020 mappings. For the 2020 LULC mapping, additional supporting materials were considered from the 2018 Corine Land Cover digital map, and from several elements of the Google Earth application software, such as 3D views and street view images available from many narrow-paved roads between villages of the study area. The visual interpretation was also supported by field sampling verifications from well-experienced human image interpreters with good knowledge of the area. The above cartographic materials were further processed using ArcGIS and Excel to create tables and digital maps of the temporal evolution of LULC types. This approach produced two digital maps of LULC types for 1945 and 2020, as well as a temporal evolution table (Figure 2).

According to Puyravad [59], the annual rate of change of LULC types was calculated for the overall study period (1945–2020). The annual rate calculation Equation (1) was based on the formula developed from compound interest law, and offers a better assessment and biological meaning to the LULC change comparisons because it is insensitive to the different time periods between observation dates [29]:

$$r=\left(\frac{1}{t_2-t_1}\right)\times \ln \left(\frac{A_2}{A_1}\right)$$

(1)

where r is the annual rate of change, and A₁ and A₂ are the LULC class areas at time t₁ and t₂, respectively.

The next phase in the process was to estimate the spatiotemporal transformations of the study area for 2020 as a result of the diachronic transitions of all LULC types from their original surfaces in 1945. This was accomplished by employing a common post-classification comparison (PCC) change detection method across the study’s periods of various dates [30]. The PCC method produced a LULC change transition matrix, which was calculated using ArcGIS overlay functions for all time periods. In addition, a map showing the spatiotemporal transition of LULC types was constructed. Additional components of land changes, such as gains and losses, net changes, total changes, and swap [60], were included in the LULC changes study due to transition matrices. The proportion of the landscape that underwent gross gain or loss of LULC type j between times 1 and 2 was represented by the letters P_+j and P_j+, respectively. The proportion of the landscape that demonstrated the persistence of category j was indicated by the diagonal elements (denoted as P_jj) of LULC types [60]. The difference between gain and loss is called net change and was denoted as D_j. Swap is the simultaneous gain and loss of LULC type j, and was calculated as two times the minimum gain and loss (S_j). The sum of the net change and the swap, or the sum of gains and losses for each LULC type j, abbreviated as C_j, is the total change [29,60]. In order to calculate net changes, swaps, and total changes, Equations (2)–(4) were applied:

$$D_j=P_{+j}-P_{j+}$$

(2)

$$S_j=2\times MIN \left(P_{j+}-P_{jj}, P_{+j }-P_{jj}\right)$$

(3)

$$C_{j }=D_j+S_j$$

(4)

Recent scientific views have defined net change as quantity difference (or quantity disagreement), and swap as allocation difference (or allocation disagreement) [37].

Identification of the most systematic transitions or dominant signals of change is another critical component in evaluating LULC alterations [29,61]. The most important form of transition can be determined using the transition matrix data by adding the total area of change for each LULC type over the time periods. This technique cannot consider the random process of LULC changes caused by the dominant LULC types and, therefore, interpreting LULC transitions based on their sizes is the correct way to evaluate them [29]. The predicted gains (denoted as G_ij) and expected losses (denoted as L_ij) that will occur if random changes among the LULC types occur, were computed using a process that was first proposed by Pontius [61] (Equations (5) and (6)):

$$G_{ij}=\left(P_{+j}-P_{jj}\right)\left(\frac{P_{i+}}{100-P_{j+}}\right)$$

(5)

$$L_{ij}=\left(P_{i+}-P_{ii}\right)\left(\frac{P_{+j}}{100-P_{+i}}\right)$$

(6)

The difference between the observed (P_ij) and expected (G_ij or L_ij) transitions in a random process of gain (P_ij − G_ij) or loss (P_ij − L_ij) is indicated as D_ij, and the ratios meaning (P_ij − G_ij)/G_ij or (P_ij − L_ij)/L_ij are denoted as R_ij. D_ij and R_ij values show the tendency of a LULC type j to gain from type i (focus on gains) and the tendency of LULC type i to lose from type j (focus on losses) [60]. Systematic transitions or dominant signals of change are indicated by values having a considerable positive or negative deviation from zero [29]. R_ij ratios are equivalent to the (observed value − expected value)/expected value ratios that are used in chi-square tests [61].

The results from the annual rate of change, absolute values of net change, and the main systematic transitions of all LULC types were used to identify the main LULC types that had undergone significant spatiotemporal changes.

3.2. Logistic Regression, Probability Maps, and Future Projection

The future projection trends of the main LULC types that were identified as experiencing the most significant spatiotemporal changes, were determined by logistic regression analyses, under the methodological approach of the CLUE modeling framework [21,62,63,64,65]. According to CLUE modeling, a set of landscape driving factors (LDF were used as independent variables in the regression analysis. In this research, 20 LDF variables were identified and collected based on the physiographic, accessibility, and socioeconomic conditions of the study area (

Table 5. Type, units, and data sources of the independent variables used in the logistic regression analyses for future projection modeling of land use changes in the study area.

a/a	Independent Variables *	Type/Unit	Data Source
1	Elevation	Continuous/m	DEM Aster 2
2	Slope	Continuous/%	DEM Aster 2
3	Alluvial deposits/very deep soils	Binary/0–1	Soil map of Greece (Nakos 1991)
4	Hard limestone/ahallow to bare soils	>>	>>
5	Limestone colluvium/seep to moderately deep soils	>>	>>
6	Doline-deposition cones/seep soils	>>	>>
7	Tertiary deposits/seep to moderately deep soils	>>	>>
8	Tertiary deposits/shallow soils	>>	>>
9	Schist/shallow soils	>>	>>
10	Schist/deep soils	>>	>>
11	Erosion potential	Continuous/t × ha−1 × year−1	Soil erosion by water (RUSLE 2015)/ESDAC **
12	Distance from unpaved roads	Continuous/m	Digital files from state Cadastre, Google Earth
13	Distance from paved roads	>>	>>
14	Distance from water courses	>>	Hydrological model from DEM, topographic maps, Google Earth
15	Distance to settlements	>>	Land use map 2020
16	Population density	Continuous/number × ha−1 of the total area	Hellenic Statistical Authority [53]
17	Sheep/cow density	>>	Hellenic Statistical Authority [53], PCAGGCA *** [56],
18	Goat density	>>	Hellenic Statistical Authority [53], PCAGGCA *** [56],
19	Annual mean temperature	Continuous/°C	https://worldclim.org/ (accessed 20 April 2022)
20	Annual precipitation	Continuous/mm	>>

* Landscape driving factors (LDF), ** European Soil Data Cent, *** Payment and Control Agency for Guidance and Guarantee Community Aid.

).

In addition to the above independent variables, the identified main LULC types were selected as dependent variables. According to the spatial module of the CLUE model, both the LULC types and the independent variables were transformed into digital raster files (ArcGrid format) with a pixel size of 100 m. The raster files of all dependent variables and 8 out of 20 independent variables (Table 5) were binarily rendered. As a result, each pixel of a given main LULC type and the eight independent variables received a value of 1, and those without received a value of 0. All the other independent variables received a continuous value according to their definition. The raster data sets were then transformed into ArcGIS ASCII grids, and with the use of the “File Convert v2” application of Dyna-CLUE modeling version [64], were further transformed into tabular format necessary for entry in the SPSS statistical package v. 27.0 (IBM Corp., Armonk, NY, U.S.A.).

In SPSS, the data were analyzed by the method of binary logistic regression of absence/presence, using forward conditional analysis as a step-by-step regression method. The input and output probabilities of the independent variables in the equation were set to not surpass the significance levels of input = 0.01 and output = 0.02, respectively, during the process [64]. Variables that did not meet the above criteria were rejected as exhibiting a low degree of correlation to the model. This procedure resulted with fewer independent variables from the original selection, and resolved problems due to multicollinearity [66,67,68]. The regression coefficients (b_i) of the remaining independent variables in the logistic equation were tabulated. Furthermore, the area under the ROC (relative operating characteristic) curve (AUC) was estimated, as a measure of controlling the goodness of fit of the data to the logistic regression model [64,69], and was used to validate the model (Figure 2) [65,70].

The AUC number indicates how well the model can differentiate across classes [71,72]. The greater the value, the better the model’s ability to distinguish between classes. AUC values range from 0 to 1, with 0.5 indicating that the model is unable to distinguish between the classes, and 1 indicating that the model is perfectly fitted [73,74,75]. AUC levels of 0.7 to 0.8 are considered acceptable, 0.8 to 0.9 are considered excellent, and values beyond 0.9 are considered exceptional [76,77]. All the available data concerning the dependent and independent variables and the logistic regression results from SPSS were introduced into the Dyna-CLUE version of the model, to produce a set of land use probability maps. Probability maps represent the distribution of the results of the logistic regression equations in the landscape [64].

The land use demands for 2040 of the identified main LULC types that have undergone significant spatiotemporal changes were computed using linear interpolation of their historical trend (Figure 2). This technique is often used to construct “Business as Usual” model scenarios (BAU scenario) [27,62]. The BAU scenario for a 20-year prediction period (2020–2040) was calculated by adding one-third of the total positive or negative trend of change from the most recent available historic trend, which were the years 1960 and 2020 (60-year trend). The areas of identified main LULC types for 1960 were obtained from the available forest vegetation and land cover map. As Mamanis and coworkers [27] suggested, the one-third ratio was used because the 20-year prediction period is equal to a third of the historical trend (20 years/60 years = 1:3). The projected land use demands under the BAU scenario were spatially allocated into the probability maps based on the higher probability of occurrences, which resulted in the creation of the predicted potential map of the future distribution of the main LULC types. The 2040 prediction maps did not consider LULC interactions.

Finally, the projected results were examined using the ArcGIS Patch Analyst program [21,45,78,79,80]. The number of patches (NumP) and mean patch size (MPS, ha) as overall measures of landscape fragmentation, edge density (ED, m/ha) as a measure of the number of ecotones [44], and mean nearest neighbor (MNN, m) as a measure of patch isolation, were calculated as indicators of spatial heterogeneity in landscape and class levels. The mathematical formulas for the specified indices can be found in the user manuals for Patch Analyst and Arc Fragstats [45,78,79].

4. Results

4.1. Land Use/Land Cover Changes. Spatiotemporal Transitions

The results of photo interpretation and LULC changes over the 75-year periods in terms of area, percentage, and annual rate of changes are shown in

Table 6. Temporal evolution and the annual rate of changes (1945–2020) of land use/land cover types in the study area (ha).

Land Use/Land Cover Types	1945	2020	Area Change from 1945 to 2020 (ha)	Percentage Change from 1945 to 2020 (%)	Annual Rate of Change (% per Year)
Agricultural areas	9581.59	8092.84	−1488.75	−15.54	−0.23
Grasslands	8052.24	5893.51	−2158.73	−26.81	−0.42
Open Shrublands	3756.35	3571.97	−184.38	−4.91	−0.07
Dense Shrublands	4390.13	4849.15	459.02	10.46	0.13
Silvopastoral areas	4671.17	2850.80	−1820.37	−38.97	−0.66
Forest	8000.30	13,430.98	5430.68	67.88	0.69
Barren areas	788.44	432.18	−356.26	−45.19	−0.80
Urban areas	317.46	447.05	129.59	40.82	0.46
Lakes	203.87	193.18	−10.69	−5.24	−0.07
Total	39,761.55	39,761.66

. The LULC types that increased in the study area were forest (68%), dense shrublands (10%), and urban areas (41%) (Table 6).

All the other LULC types decreased, with the more important changes being the reduction in silvopastoral areas (39%), grasslands (27%), and agricultural areas (16%) (Table 6). Open shrublands were reduced in area to a limited extent (5%). Barren areas and lakes were also reduced by 45% and 5%, respectively, but they covered only a limited part of the study area. Forest and dense shrubland expansion, at the expense of silvopastoral areas, grasslands, and agricultural areas, demonstrated that in the last 75 years, woody vegetation in the study area had significantly increased. Analyzing in more detail the annual rate of changes of LULC types (Table 6) during the study period 1945–2020, suggested that the most significant changes were the declining trend of barren areas, silvopastoral areas, and grasslands, and the increasing trend of forest and urban areas. Agricultural areas presented a considerable declining trend, but were less severe compared with grasslands and silvopastoral areas.

Gradual conversion of silvopastoral areas and grasslands into forest was observed in the northern areas between the villages of Feneos, Tarsos, Karya, and Trikala (Figure 3). Forest also seemed to have expanded in the southern area near the village of Drosopigi at the expense of shrubland areas. Grasslands decreased over time in the study areas, except for the central area east of the village Goura, where a higher elevation of landscape occurs (>1200 m). Changes in agricultural areas did not have a strong spatiotemporal orientation, suggesting that they covered a broader range of landscape territories.

The LULC transition matrix of the study area showed that between 1945 and 2020, 65.71% of the total landscape remained unchanged, while 34.29% was transformed into a different LULC type (

Table 7. Land use/land cover (LULC) change transition matrix between 1945 and 2020 (%) in the study area.

2020 LULC
1945 LULC	AG	GR	OS	DS	SI	FO	BA	UR	LA	Total 1945	Loss
AG	19.16	1.12	1.05	1.10	0.35	0.97	0.02	0.26	0.08	24.11	4.95
GR	0.41	10.94	2.21	0.67	2.57	3.27	0.15	0.03	0.00	20.25	9.31
OS	0.06	0.68	4.65	1.81	0.41	1.78	0.01	0.04	0.00	9.44	4.79
DS	0.20	0.17	0.51	7.80	0.20	2.12	0.02	0.02	0.00	11.04	3.24
SI	0.09	0.90	0.38	0.35	2.85	7.16	0.00	0.01	0.00	11.74	8.89
FO	0.04	0.78	0.04	0.24	0.68	18.30	0.00	0.00	0.04	20.12	1.82
BA	0.24	0.24	0.14	0.20	0.10	0.17	0.89	0.01	0.00	1.99	1.10
UR	0.03	0.00	0.00	0.00	0.00	0.01	0.00	0.75	0.00	0.79	0.04
LA	0.13	0.00	0.00	0.02	0.00	0.00	0.00	0.00	0.37	0.52	0.15
Total 2020	20.36	14.83	8.98	12.19	7.16	33.78	1.09	1.12	0.49	100.00	34.29
Gain	1.20	3.89	4.33	4.39	4.31	15.48	0.20	0.37	0.12	34.29

Note: The values in the shaded box (diagonally) indicate the unchanged LULC types from 1945 to 2020. The underlined values indicate the most important land use/land cover transitions (>1%). AG: agricultural areas, GR: grasslands, OS: open shrublands, DS: dense shrublands, SI: silvopastoral areas, FO: forest, BA: barren areas, UR: urban areas, LA: lakes.

). According to the matrix, the most important LULC transitions (>2%) were those of silvopastoral areas (SI) into forest (FO); of grasslands (GR) into open shrublands (OS), silvopastoral areas (SI) and forest (FO); and finally, of dense shrublands (DS) into forest (FO). Additional significant changes (1–2%) were presented in agricultural area transition into grasslands (GR) and shrublands (open, OS; or dense, DS); and open shrublands (OS) transition into dense shrublands (DS) and forest (FO). The results of the matrix suggested that the most important changes in the Mt Zireia landscape were woody plant encroachment into open areas such as grasslands, open shrublands, and silvopastoral areas, and to a lesser extent, into agricultural areas (Table 7).

Figure 4 presents the map of LULC transitions in the study area. It is notable that silvopastoral and forest expansions were mainly located in the northern parts of the study area. On the other hand, shrubland expansions were mainly observed in the southern and central parts of the landscape. Overall, LULC transitions were observed in all parts of the Zireia landscape, but appeared to be more extensive in the northern parts.

The most significant changes in net values (absolute values) were observed in forest, grasslands, and silvopastoral areas, and to a lesser extent, in agricultural areas (

Table 8. Temporal evolution of gain, losses, net change, and swap of land use/land cover, in terms of percent, in the landscape of Mt Zireia for the period 1945 to 2020.

Land Use/Land Cover	Percentage of Change
	Gain	Loss	Total Change	Swap	Absolute Value of Net Change
Agricultural areas	1.20	4.95	6.15	2.40	3.75
Grasslands	3.89	9.31	13.20	7.78	5.42
Open shrublands	4.33	4.79	9.12	8.66	0.46
Dense shrublands	4.39	3.24	7.63	6.48	1.15
Silvopastoral areas	4.31	8.89	13.20	8.62	4.58
Forest	15.48	1.82	17.30	3.64	13.66
Barren areas	0.20	1.10	1.30	0.40	0.90
Urban areas	0.37	0.04	0.41	0.08	0.33
Lakes	0.12	0.15	0.27	0.24	0.03
Landscape	34.29	34.29	34.29	19.15	15.14

). These data also indicated that the highest losses in the area were observed in grasslands and silvopastoral areas, and the highest gains in forest areas. Net change values in forest were much higher than in comparison with their swap values, suggesting that forest expansion in new areas (quantity difference) was more significant than their simultaneous exchange of forest areas to other uses (allocation difference). On the contrary, net change and swap value changes appeared to be more balanced in grasslands and silvopastoral areas. Forest, grasslands, and silvopastoral areas were the main LULC types that underwent significant changes, and the recorded woody plant expansion in the landscape focused particularly on forest development (Table 8).

Table 9. Area percentage of the main systematic transitions of land use/land cover changes in terms of gains and losses in the landscape of Mt Zireia for the period 1945 to 2020.

Transitions	Gains (%)		Losses (%)
	Dij	Rij	Dij	R
AG to GR	−0.06	−0.05	0.20	0.22
AG to OS	−0.11	−0.09	0.49	0.88
AG to DS	−0.09	−0.08	0.34	0.45
GR to OS	1.24	1.28	1.23	1.25
GR to SI	1.58	1.59	1.79	2.28
GR to FO	−0.65	−0.17	−0.42	−0.11
OS to DS	1.34	2.87	1.17	1.81
OS to FO	−0.05	−0.03	0.00	0.00
DS to FO	−0.02	−0.01	0.87	0.70
SI to FO	4.88	2.14	3.92	1.21

AG: agricultural areas, GR: grasslands; OS: open shrublands; DS: dense shrublands; SI: silvopastoral areas; FO: forest; BA: barren areas; UR: urban areas; LA: lakes; Dij: the difference between the observed and expected transitions; Rij: the difference between the observed and expected transitions, relative to the expected transitions.

presents the percentage of the area of the main systematic transitions of LULC changes in terms of gains and losses. The largest positive or negative variation from zero (systematic transitions) appeared to be in the transitions of silvopastoral to forest (SI to FO), grasslands to silvopastoral (GR to SI), grasslands to open shrublands (GR to OS), and open to dense shrublands (OS to DS). These results suggested that forest was systematically gaining area from silvopastoral areas, while dense shrublands gained from open ones (focus on gains). On the other hand, the same results showed that grasslands were systematically losing areas to open shrublands and silvopastoral areas (focus on losses).

Overall, spatiotemporal changes in the landscape of Mt. Zireia indicated that the most important element of landscape change was the woody plant expansion into open areas. Furthermore, among the different LULC types, the ones that were mainly affected by landscape changes were forest, grasslands, and silvopastoral areas in terms of area, percent, and the annual rate of changes (Table 6), the absolute value of net change (Table 8) and the main systematic transitions (Table 9).

4.2. Logistic Regression, Probability Maps and Future Projection of Forest, Grasslands and Silvopastoral Areas

The logistic regression analyses (forward conditional–stepwise) revealed the influence of each of the 20 included independent variables on the LULC types of forest, grasslands, and silvopastoral areas (dependent variables). The area under the ROC curve (AUC) is presented in Figure 5.

AUC values for the three main LULC types were above 0.8 for forest and grasslands, which is accepted as an excellent discrimination, and 0.697 (equal to almost 0.7) for silvopastoral areas, which is marginally acceptable discrimination [77], indicating that the logistic regression models possessed significant goodness of fit.

The b-values of the independent variables are presented in

Table 10. The logistic regression coefficients (b-values and constant).

Independent Variables	Forest	Grasslands	Silvopastoral
Elevation (m)	0.0010	0.0024	−0.0007
Slopes (%)	0.0306	−0.0289	0.0122
Alluvial deposits/very deep soils	−2.8321	−2.6520	−4.1481
Hard limestone/shallow to bare soils	−0.6249		−0.4752
Limestone colluvium/deep to moderately deep soils	0.6755
Doline-deposition cones/deep soils	−19.1054	−1.3975	−18.2201
Tertiary deposits/deep to moderately deep soils		−0.3152	0.6172
Tertiary deposits/shallow soils	−0.6550	−0.6041
Schist/shallow soils	1.6790		−0.4972
Schist/deep soils	2.7273	−1.7426	−2.2187
Erosion (t/ha/year)	−0.2248	0.0324	−0.0059
Distance from unpaved roads (m)	0.0001	−0.0003
Distance from paved roads (m)	0.0005
Distance from water courses (m)	−0.0002	−0.0003
Distance to settlements (m)	−0.0001	−0.0002	−0.0003
Population density (number of people per area (ha)	−7.4041	1.5776
Sheep/cow density (number of heads per area (ha)	−0.3747	0.1150	0.2747
Goat density (number of heads per area (ha))	−3.0946
Annual mean temperature (Co)	−0.3339		−0.2338
Annual precipitation (mm)	−0.0113	0.0097	0.0052
Constant	12.8058	−11.7581	−3.7588

. The cells without data indicate the independent variables which did not show a statistically significant correlation with the LULC types. The regression coefficients of Table 10 (b-values and constant) were entered into the CLUE software environment to build a set of three probability maps and finally complete the landscape change prediction procedure. The produced maps of the probability (%) of future occurrence for forest, grasslands, and silvopastoral areas in the study area are presented in Figure 6.

According to the probability maps (Figure 6), forest possessed a higher possibility of occurrence mainly in the north-northwest parts of the study area, and to a lesser extent, in a restricted area in the south. Grasslands, on the other hand, were found to be highly possible to occupy areas of high altitude in the center part of the study region, covering the grounds of the sub-alpine zones near Mt Zireia’s summit. Silvopastoral areas received a lower chance of occurrence compared with forest and grasslands, and these chances of occurrence were scattered around the landscape.

Land use demands of the three LULC types under the BAU scenario for a 20-year prediction period (2020–2040), are presented in

Table 11. Area distribution (ha) and rate of change (%) of forest, grasslands, and silvopastoral areas in the study area for the projection period 2020–2040.

Land Use	Area (ha)		Rate of Change %
	2020	2040
Forest	13,430.98	15,474.87	15.22
Grasslands	5893.51	5228.86	−11.28
Silvopastoral areas	2850.80	1979.33	−30.57
Total	22,175.29	22,683.06	2.29

. According to Table 11, if the main factors of change (landscape changing factors–independent variables) continue to be the same between 2020 and 2040, then forest is expected to increase in some areas by 15%, and grasslands and silvopastoral areas to decrease by 11% and 31%, respectively.

Spatial allocation of the projected land use demands into the probability maps (Figure 6) produced the predicted potential map (Figure 7) of the spatiotemporal distribution for forest, grasslands and silvopastoral areas for the 2020–2040 period. Forest is expected to continue expanding in the north-northwest parts, and probably will occupy scattered new areas in the central parts of the landscape. The projected grasslands reduction, on the other hand, will most probably force the remaining grassland patches to be limited to the central parts of the landscape. Silvopastoral areas will probably continue to occupy small, scattered areas around the landscape, but with a spatial distribution uneven in size.

The above results can be confirmed by evaluating the landscape structure of the projected maps with the help of landscape metrics (

Table 12. Landscape metric evaluation for forest (FO), grasslands (GR), and silvopastoral areas (SI) in the landscape of Mt Zireia for the projection period 2020–2040.

	NumP 1			MPS 2			ED 3			MNN 4
	FO	GR	SI	FO	GR	SI	FO	GR	SI	FO	GR	SI
2020	30	155	70	448.70	38.11	40.21	14.80	13.48	8.22	348.96	229.83	419.34
2040	200	155	123	77.39	33.90	15.95	20.49	7.74	5.78	192.68	289.07	268.55

¹ Number of patches. ² Mean patch size (ha). ³ Edge density (m/ha). ⁴ Mean nearest neighbor (m).

). According to the metrics during the projected period (2020–2040), the expansion of forests will probably increase their overall fragmentation in the sense that numerous new and smaller sized patches (indicated by the NumP and MPS values) will be created in new areas across the landscape (Table 12). That increase will improve the ED value, creating new forest edges and promoting forest connectivity, as was indicated by the decrease in their MNN value.

Grassland patches, on the other hand, will probably become smaller in size (decrease in MPS value) and distant from each other, as was indicated by the increase in their MNN value. The latter observation will probably promote patch isolation of the smaller grassland units which occupy the marginal areas around their main distribution in the center of the landscape (Figure 7). Furthermore, the decrease in ED value indicates that a significant reduction in grasslands ecotone is to be expected. Silvopastoral patches are expected to become more fragmented in the future, meaning greater in numbers but smaller in size (indicated by the NumP and MPS differences). Moreover, even though silvopastoral areas would increase their overall connectivity (MNN value), they are expected to greatly reduce their edges (ED value), similar to grasslands.

5. Discussion

Spatiotemporal transition analysis of the landscape of Mt Zireia suggested that the strongest trend of landscape evolution was woody plant expansion in open areas, and grasslands reduction. Among the different LULC transitions, the most systematic ones (Table 9) were forest expansion over silvopastoral areas, of open shrublands over dense shrublands, and of grassland reduction in favor of open shrublands and silvopastoral areas. These results, combined with the data of the total LULC changes in area, percent, and the annual rate and net changes (Table 6 and Table 8), suggested that the LULC types that are mainly affected by landscape changes are forests, grasslands and silvopastoral areas. This finding, especially as far as the forest expansion/grasslands reduction trend is concerned, is in line with similar studies conducted in Greece [15,24,25,28] and other Mediterranean countries [2,47,55,81], indicating that special focus should be provided to these specific LULC interactions, especially in the rapidly-changing Mediterranean landscapes [9].

The above trend of LULC interactions can mainly be attributed to land abandonment issues related to socioeconomic conditions. Relevant inventory data from the village communities in the study area (Table 1 and Table 2) suggested that socioeconomic changes over the previous decades in the study area had the form of a decrease in local population, population aging and a significant temporal reduction in the percentage of employees in the primary economic sector. These specific types of socioeconomic changes are reported to especially occur in Mediterranean landscapes, as directly related to the abandonment of traditional management practices, such as extensive or semi-extensive pastoral activities (including transhumance pastoralism), wood product collection (e.g., coal and fuel woods) and agricultural practices in less favorable areas (e.g., crop fields in terraces) [8,47,54,55,81,82,83].

Additional inventory data related to the numbers and farms of grazing animals from the study area supports the notion that land abandonment has affected pastoral activities (Table 3 and Table 4). More specifically, the number of sedentary and transhumant grazing animals and farms for the study area were significantly reduced over the previous decades, following a similar trend of change as the number of people that were employed in the primary sector of the economy. The reduction in grazing animals was also reported to follow a similar, more general, trend of reduction for the whole country and for other south European Mediterranean countries [8,11,81]. The collection of forest products seems to be affected by the abandonment of traditional practices. Unpublished data from the PACTORES Project (www.pactores.eu (accessed on 10 December 2021)) indicated that fuel wood collection by local people within the study area has significantly reduced over the years, and in some areas has practically stopped. On the other hand, some of the forest expansion over open areas can be attributed to the afforestation policies of the local forest service to increase the area covered by high forests. Finally, agricultural activities were also affected by land abandonment, but this effect was less severe on the extent of agricultural lands. According to the spatiotemporal analysis of this research, agricultural areas scored as the fourth most important LULC change in terms of total area, percent of change, annual rate, and net changes (Table 6 and Table 8) and these changes did not appear to have a strong geographic orientation (Figure 3). Agricultural activities were mainly oriented in plains in favorable and more accessible parts of the Mt Zireia landscape, which, as similar studies have pointed out, are probably less affected by land abandonment [11,84]. All these aspects of socioeconomic effects in the current management of Mediterranean landscapes have been noted throughout the Mediterranean region of Europe [2,13,15,22,24,25,54,55,81,85] and have been identified as the main reasons for landscape change.

Future development for forest, grasslands, and silvopastoral areas based on the BAU scenario of linear extension of land use demands for 2040 and probability maps, suggested that forest will most probably continue to expand in the north-northwest parts, adding new areas scattered mainly in the central parts of the landscape. At the same time, grassland and silvopastoral areas will probably continue to reduce in area, occupying territories mainly at the central part for grasslands, or small scattered territories around the landscape for silvopastoral areas (Figure 7). Similar results of future development of forests and grasslands were very recently reported from a rural landscape study of central Greece under a similar trend of land abandonment [27,65]. Evaluation of the structural developments (Table 12) of the Mt Zireia landscape from the projected maps revealed that forest expansion into new areas, and in many cases, as small patches, will increase their overall dispersal and will create new forest edges (higher ED value). Furthermore, the decrease in MNN value will promote forest connectivity. Grasslands, on the other hand, apart from occupying one large and three smaller core areas in the center and the north-northwest part of the study area (Figure 7), will probably keep only smaller, fragmented, and isolated patches around the landscape, as indicated by the reduction in their MPS and the increase in their MNN values. Moreover, the decrease in ED value will cause a significant reduction in grassland edges. Silvopastoral patches, similar to grasslands, will became more fragmented with reduced edges. These findings correlate with the response of many other landscapes around the world, showing that forest expansion usually leads to increased forest patch connectivity promoting forest edges, while open habitat reduction usually creates the opposite trend of a reduction in connectivity and edges [86,87]. These results could be alarming for sustaining the environmental integrity of the Mt Zireia landscape, as many researchers have linked grasslands fragmentation and the loss of connectivity and boundary lengths of open habitats, to the decline of species richness and mountainous biodiversity [47,86,88]. The results of landscape metrics evaluation on the future development of LULC types can serve as evidence of the great dynamic of expansion that forest patches possess over grassland and silvopastoral patches, independent of the environment and biodiversity implications.

The findings of this study are consistent with the common pattern of woody cover expansion over open regions in many Mediterranean landscapes that suffer from land abandonment [89]. Environmental integrity, biodiversity, and cultural heritage may be positively or negatively impacted by land abandonment, which can additionally benefit forest ecosystems by fostering it at minimal cost and on a larger scale [65]. Forest recovery promotes carbon sequestration, erosion reduction, and several other ecosystem services such as climate and water regulation, wood production, and recreation [12]. On the other hand, land abandonment, especially in the Mediterranean region, results in declining biodiversity and loss of traditional cultural landscapes [7,81,84], and is often linked to an increased risk of wildfires and decreased river flows [83,90,91]. The land abandonment effect can also be associated with the loss of important cultural elements and services related to traditional pastoral activities, such as cultural heterogenic pastoral landscapes, gastronomical heritage, and folklore elements [92,93].

Developmental planning must take into consideration the spatiotemporal trends and the future projection of LULC types recorded in this study. Forest expansion over grassland and silvopastoral areas, apart from the environmental and cultural implications, would have a strong negative effect on the future of sustainable development of livestock husbandry in the study area. Grassland and silvopastoral areas are considered important natural resources necessary for applying extensive pastoral practices, especially the transhuman livestock system, and the reported threat status could have a damaging effect on keeping the ecological integrity and the social benefits that people expect from pastoral landscapes [62,94,95].