Are Agrarian Areas in Mediterranean Mountain Regions Becoming Extinct? A Methodological Approach to Their Conservation

Mediterranean mountain regions have undergone several landscape changes since the end of the 19th century due to progressive depopulation and the abandonment of cattle rearing, forestry, charcoal production and agricultural activity. Such activity favored landscape dynamics by creating grassy habitats, which in turn resulted in greater landscape diversity. This is now being lost as the forest reclaims abandoned pastures. Thus, the purpose of this work was to identify those open habitats most in need of management action to maximize biodiversity and cultural heritage conservation and minimize fire risk and management costs. These analyses show a sharp decrease of open agriculture areas, which are the habitat of many endemic species (from 46.4% to 12.3%), currently overgrown with secondary forests. Multivariate analysis and the PGP (Patch Growing Process) heuristic model indicate the areas in which the restoration of open areas (by about 8%; about 500 ha) will be the most advisable and the most beneficial, taking into account environmental, social and economic factors. The use of PGP provides for a 21% improvement in total agriculture areas. Still, the natural state of the protected Mediterranean mountain area “Alta Garrotxa” (Catalonia, Spain) is almost continuous forest. However, the management models proposed in this study offer flexible precepts to achieve the desired landscape patterns and maintain biodiversity, while conserving cultural heritage and decreasing the risk of fire.


Introduction
The appearance of a territory today is the result of management by different generations over time and changes in political and agrarian pressures that act directly or indirectly on the landscape. This alteration in management over the years (whether manual, animal or mechanized) has created exceptional cultural landscapes of a high aesthetic, cultural and ecological value [1].
However, in the early 20th century and particularly since World War II, the relationship between agriculture and nature was dramatically transformed due to an increase in production systems and the abandonment of rural areas [2]. This abandonment of traditional activity in agricultural spaces has triggered serious transformations in the landscape [3], characterized by a continuous and progressive forest and scrub colonization, which has spread across the majority of what were previously open spaces. These changes can be summarized in two types of dynamics: on the one hand, landscape homogenization due to forest expansion and on the other, the fragmentation of agricultural areas due to said expansion [4,5]. This sharp increase the amount of forest area in most industrialized countries may represent a positive process in terms of counterbalancing trends in forest overexploitation in other parts of the world. But the uncontrolled growth of new forest areas represents two problems: the loss of a mosaic landscape constructed by agricultural activity and a wide range of habitats of community interest, as well as an increased risk of forest fire [6][7][8][9][10].
The abandonment of livestock, forestry and agricultural activity due to a lack of economic profitability has led to significant landscape and environmental changes in the territory, such as the increase in forest surface area and accumulation of forest biomass, which has led to a clear homogenization of mountain landscapes, especially in the case of the mid-sized Mediterranean mountains [11]. The loss of agricultural areas (open areas), especially pastures, in subalpine and alpine areas has not been so drastically threatened due to more favorable bioclimatic conditions, which have allowed for the maintenance of quality pasture, guaranteeing a certain level of economic returns. This is much less likely to happen in the Mediterranean bioclimatic domain, which only allows for the existence of low-quality pasture. Landscape dynamics and their possible effects on the territory have become an important focus of study and work in the last decade [1,[12][13][14][15].
Expanding our knowledge of landscape patterns through spatial and temporal study makes it possible to identify the aforementioned homogenization/fragmentation dynamics, as well as environmental consequences linked to the appearance of erosion, infiltration and runoff processes, changes in the behavior of hydrological and geomorphological processes, increased fire risk and the loss of biodiversity and cultural heritage [16][17][18][19][20] European mountain landscapes, which boast a rich heritage in terms of natural and cultural values, have been threatened by the forces of change, putting them under pressure and causing them to lose their identity. Landscape ecology, a more holistic trend in the study of these landscapes, has facilitated their study, analysis and understanding, which will in turn improve their conservation [21][22][23]. It is important to remember that the historical study of this cultural heritage often allows us to understand the reason for the existence of current biodiversity, how it has evolved over the years and what the future may hold. This is why the sociocultural study of these landscapes rich in cultural symbols and which are the fruit of years of human activity cannot be considered separately from studies on the natural environment [19,21,24].
The common agricultural policy (CAP) supports the vibrancy and economic viability of rural communities through rural development measures. The European agricultural fund for rural development (EAFRD) is the funding instrument of the CAP that supports rural development strategies and projects. It also forms part of the European structural investment funds (ESIF). The EAFRD budget for the 2014-20 period amounts to around €100 billion. The budget will be spent over the course of this period, through the implementation of rural development programs which run until the end of 2023. It is distributed according to six priorities: fostering knowledge transfer and innovation in agriculture, forestry and rural areas; enhancing the viability and competitiveness of all types of agriculture and promoting innovative farm technologies and sustainable forest management; promoting food chain organization, animal welfare and risk management in agriculture; promoting resource efficiency and supporting the shift toward a low-carbon and climate resilient economy in the agriculture, food and forestry sectors; restoring, preserving and enhancing ecosystems related to agriculture and forestry; promoting social inclusion, poverty reduction and economic development in rural areas.
Within the framework of the new financial perspectives, the European Union has established criteria to guide rural development plans [25] and created the new EAFRD fund (European Agricultural Fund for Rural Development). The European Commission is interested in funding part of the costs of the Natura 2000 Network [26] through various measures that will have to be included in rural development plans approved at national or regional level. At the same time, areas of natural interest will also enjoy financial support through the new Life+ program. Rural development policy is closely related to these concepts: diversity (multifunctionality); sustainability and climate change mitigation and adaptation. This new trend to promote the multifunctionality of land use can improve the viability of currently marginalized rural mountain areas and can be consolidated through international and national policy decisions applied at the local level. Therefore, tackling the continuing loss of these agroforestry landscapes has become a priority for administrations, scientists and citizens alike. And this priority needs to be translated into greater effort when it comes to studying, protecting, maintaining and even, as far as possible, reclaiming them.
Mediterranean mountain areas have undergone a progressive process of land abandonment and this economic and demographic marginalization has only occasionally been interrupted in recent decades by the appearance of new activities related to tourism and second residential uses. Agricultural areas in mountain regions are subject to a continuous process of fragmentation. Fragmentation can have profound impacts on the landscape, especially with respect to the persistence of certain species [16,27,28]. The effects of habitat loss and fragmentation are considered the main threat to global biodiversity loss [16,[29][30][31]. The loss of open spaces in the study area has had a major impact on certain taxa of vascular flora of interest, including Polygala vayredae, Allium pyrenaicum, Erinacea anthyllis, Dryopteris remota, Narcissus poeticus and Iberis sempervirens, which could all disappear due to this process of landscape uniformity. As for fauna, a significant number of bird species included in the Birds Directive are also threatened by this process, as may be the case with Aquila chrysaetos, Pernis apivorus, Bubo bubo, Lanius collurio, Circaetus gallicus, Hieraaetus pennatus, Falco peregrinus and Caprimulgus europaeus.
Other negative effects of depopulation and homogenization in the Mediterranean mountain landscapes include the degradation of the region's rich cultural heritage (e.g., country houses, archaeological sites, churches and hiking paths), the loss of folk knowledge associated with traditional economic activities (making charcoal, hunting, collecting medicinal plants and picking fruit and traditional agriculture systems, among others) and a decrease in the general appreciation of the landscape in abandoned areas.
Traditional ecology indices [15] have often been used to monitor and evolve these fragmentation processes. Requiring experience for their interpretation, they are sometimes insensitive to the ecological processes that take place within the landscape. The results of applying quantitative methods in landscape ecology are grouped into so-called landscape metrics. These indices provide interesting numerical data on the composition and configuration of landscapes, the proportion of each type of land cover or use and the surface area and shape of landscape elements. In addition, landscape indices allow for a useful and interesting comparison between different landscape configurations appearing at different times in the same territory and the definition of future scenarios [19,[32][33][34].
Uncontrolled growth of new forest area also can be seen of having a potential threat to the landscape and that agrarian activity can actually even benefit the landscape structure. Applying a Multi-Criteria Analysis and a Patch-Growing Process (PGP) method, this study develops patches of different sizes and suitability values, which represent open spaces that should be preferably maintained and/or recovered.

Study Area
Alta Garrotxa stretches across the north of the Garrotxa region, the northeast of Ripollès and the northwest of Alt Empordà. From a geographical point of view, this area is the easternmost part of the southern Pyrenees. Pyrenees are a mountain chain of southwestern Europe that consists of flat-topped massifs and folded linear ranges. It stretches from the shores of the Mediterranean Sea on the east to the Bay of Biscay on the Atlantic Ocean on the west. The Pyrenees form a high wall between France and Spain. The range is some 430 km long (see Figure 1). It is a space delimited by two main basins-the Fluvià, which is the larger of the two and the Muga, which is only found at its head. The watershed between the two basins also defines the border between Spain and France. It has two clearly differentiated areas in climatic terms. On the one hand, the north-western area, with abundant rainfall and relatively low winter temperatures corresponding to an alpine biome; and on the other, the south-eastern area, with more moderate rainfall and milder temperatures throughout the winter, corresponding to a Mediterranean biome. Between the two extremes there is an entire climatic range, which is favored by the outstanding geomorphological complexity of the territory. Generally speaking, the dominant vegetation in this sector is the holm oak, in its lowland (Viburno-Quercetum ilicis "Quercetum ilicis galloprovinciale") and mountain variants (Asplenio-Quercetum ilicis "Quercetum mediterraneo-montanum"). These holm oak groves, which are very thick, often cover the rock slopes and shelves with a dense carpet, while in those areas touched by the sun, Martinique oak groves grow (Buxo-Quercetum pubescentis) up to the top of the highest peaks. With regard to non-forest habitats, in Alta Garrotxa we can differentiate between two main types of pasture-calcareous dry pasture (Aphyllantion) and calcareous mesophilic pasture (Festuco-Brometea).
The landscape of Alta Garrotxa is unique due to its relatively isolated position surrounded by mountains, the splendor of its relief and leafy vegetation. This has led to it maintaining a high diversity of habitats and most of this landscape unit is therefore included in Catalonia's Areas of Natural Interest Plan.
The protected natural area of Alta Garrotxa is considered to be an area of natural interest (Bill 12/1985 on Natural Areas), a level of protection equivalent to IUCN category VI. In 2005, it became part of the Natura 2000 Network.
From a biogeographical point of view, this territory corresponds to eminently forested mountain habitats of a Mediterranean (holm oak) and sub-Mediterranean (oak) nature, with the presence of beech trees in the shade.

Methods
This research is divided in 3 sections (see Figure 2). The first section (Analysis prior modeling) consist of calculating landscape configuration and suitability scenarios. The second section (Model parameters) analyze the landscape patterns of all the protected natural areas of the Catalan Pyrenees. Prior to running the model, questions such as how much surface area is to be recovered, what size the patches should be, what the configuration should be (in terms of density, distance, shape, compactness) and the distribution of these generated patches were answered. Finally, third section present the application of heuristic techniques, in particular, in the location of areas (site location) that meet certain requirements. Generally speaking, the dominant vegetation in this sector is the holm oak, in its lowland (Viburno-Quercetum ilicis "Quercetum ilicis galloprovinciale") and mountain variants (Asplenio-Quercetum ilicis "Quercetum mediterraneo-montanum"). These holm oak groves, which are very thick, often cover the rock slopes and shelves with a dense carpet, while in those areas touched by the sun, Martinique oak groves grow (Buxo-Quercetum pubescentis) up to the top of the highest peaks. With regard to non-forest habitats, in Alta Garrotxa we can differentiate between two main types of pasture-calcareous dry pasture (Aphyllantion) and calcareous mesophilic pasture (Festuco-Brometea).
The landscape of Alta Garrotxa is unique due to its relatively isolated position surrounded by mountains, the splendor of its relief and leafy vegetation. This has led to it maintaining a high diversity of habitats and most of this landscape unit is therefore included in Catalonia's Areas of Natural Interest Plan.
The protected natural area of Alta Garrotxa is considered to be an area of natural interest (Bill 12/1985 on Natural Areas), a level of protection equivalent to IUCN category VI. In 2005, it became part of the Natura 2000 Network.
From a biogeographical point of view, this territory corresponds to eminently forested mountain habitats of a Mediterranean (holm oak) and sub-Mediterranean (oak) nature, with the presence of beech trees in the shade.

Methods
This research is divided in 3 sections (see Figure 2). The first section (Analysis prior modeling) consist of calculating landscape configuration and suitability scenarios. The second section (Model parameters) analyze the landscape patterns of all the protected natural areas of the Catalan Pyrenees. Prior to running the model, questions such as how much surface area is to be recovered, what size the patches should be, what the configuration should be (in terms of density, distance, shape, compactness) and the distribution of these generated patches were answered. Finally, third section present the application of heuristic techniques, in particular, in the location of areas (site location) that meet certain requirements.

Analysis Prior to Modeling
In this study was mapped, quantified and described the composition and configuration of the landscape with a high accuracy (25 m pixel) through digital mapping and landscape indices [35] in the period 1957-2017. Data analysis was done using the ArcGIS 10.5 Licensed Software. To quantify and value the different land use and land cover changes produced between 1957 and 2017, detailed cartographic material was required. For the first period, I used black and white aerial photographs captured by the United States Army Map Service and these were enlarged to an approximate scale of 1:7500. These aerial photographs had to be corrected using ERDAS IMAGINE (Orthobase module) and then, through photointerpretation, the different land uses and covers were digitalized. For the second period I used a black and white ICGC (Cartographic and Geological Institute of Catalonia, Barcelona, Spain) orthophotomap (scale 1:5000). To verify the orthophoto map interpretation, the information was complemented with fieldwork.
Also, the results of previous research published by Varga et al. (2018) [35] were used and made part of the project "Land use and cover change in Alta Garrotxa-evolution, consequences and future scenarios" supported by Consorci of Alta Garrotxa. A set of different scenarios have been designed using a multi-criteria evaluation and geospatial information available for the study area to identify the key areas for management action and to predict the potential effects on agricultural lands by prioritizing one or another management objective.
Multicriteria assessment methods are often integrated into GIS and used in optimal location analysis for a given activity, allowing continuous suitability maps to be obtained quickly. One of the key aims of multicriteria analysis is to assist decision-making centers in describing, evaluating, ordering, hierarchizing, selecting or rejecting objects, based on an evaluation expressed by scores, values or intensities of preference and according to various attributes or criteria [36,37].
In our study, five types of criteria were established related to biodiversity, landscape structure and perception, cultural heritage, fire hazard and management cost. Different scenario maps were projected using multicriteria evaluation and geospatial information available for the study area. Each scenario map was classified into seven suitability categories (Category 1 being of low suitability and Category 7 high suitability for agrarian lands). A neutral scenario (see Figure 3) was used, whereby each of the five aforementioned criteria is awarded a weight of 20%.

Analysis Prior to Modeling
In this study was mapped, quantified and described the composition and configuration of the landscape with a high accuracy (25 m pixel) through digital mapping and landscape indices [35] in the period 1957-2017. Data analysis was done using the ArcGIS 10.5 Licensed Software. To quantify and value the different land use and land cover changes produced between 1957 and 2017, detailed cartographic material was required. For the first period, I used black and white aerial photographs captured by the United States Army Map Service and these were enlarged to an approximate scale of 1:7500. These aerial photographs had to be corrected using ERDAS IMAGINE (Orthobase module) and then, through photointerpretation, the different land uses and covers were digitalized. For the second period I used a black and white ICGC (Cartographic and Geological Institute of Catalonia, Barcelona, Spain) orthophotomap (scale 1:5000). To verify the orthophoto map interpretation, the information was complemented with fieldwork.
Also, the results of previous research published by Varga et al. (2018) [35] were used and made part of the project "Land use and cover change in Alta Garrotxa-evolution, consequences and future scenarios" supported by Consorci of Alta Garrotxa. A set of different scenarios have been designed using a multi-criteria evaluation and geospatial information available for the study area to identify the key areas for management action and to predict the potential effects on agricultural lands by prioritizing one or another management objective.
Multicriteria assessment methods are often integrated into GIS and used in optimal location analysis for a given activity, allowing continuous suitability maps to be obtained quickly. One of the key aims of multicriteria analysis is to assist decision-making centers in describing, evaluating, ordering, hierarchizing, selecting or rejecting objects, based on an evaluation expressed by scores, values or intensities of preference and according to various attributes or criteria [36,37].
In our study, five types of criteria were established related to biodiversity, landscape structure and perception, cultural heritage, fire hazard and management cost. Different scenario maps were projected using multicriteria evaluation and geospatial information available for the study area. Each scenario map was classified into seven suitability categories (Category 1 being of low suitability and Category 7 high suitability for agrarian lands). A neutral scenario (see Figure 3) was used, whereby each of the five aforementioned criteria is awarded a weight of 20%.

Model Parameters
In order to establish these parameters, landscape patterns between forest and non-forest habitats were studied for the period 1987-2017, using landscape ecology indices and the Fragstats public domain program [15]. To carry out this study, raster data (30-m cell size) were obtained for land use and cover (22 categories) from the Department of the Environment and Housing. These data were captured using the Thematic Mapper (TM) sensor from the Landsat satellite. The characterization of landscape patterns was focused on the Catalan Pyrenees.
The aim was to group those protected natural areas present in the same geographical unit (in this case the Catalan Pyrenees) to more effectively and realistically compare the landscape patterns experienced in the last fifteen years with those of the study area.
A total of thirteen variables (landscape indices) were used to assess landscape patterns and changes in landscape composition and configuration between forested and non-forested areas in the seven protected natural areas every five years from 1987 to 2002. Those variables were [15] see Table  A1: Another aspect that must be taken into account in this analysis is the conjunction of two welldifferentiated biogeographical regions throughout the Pyrenees: The Alpine and the Mediterranean (see Figure 4). To assess whether or not these patterns are of a purely biogeographical nature, we first conducted a multivariate analysis to simplify or reduce the dimension-a principal component analysis (PCA) on a total of thirteen variables (landscape indices) used to assess landscape patterns and changes in landscape composition and configuration between forested, transitional and nonforested areas in the seven protected natural spaces for the following years 1987, 1992, 2002 and 2017.

Model Parameters
In order to establish these parameters, landscape patterns between forest and non-forest habitats were studied for the period 1987-2017, using landscape ecology indices and the Fragstats public domain program [15]. To carry out this study, raster data (30-m cell size) were obtained for land use and cover (22 categories) from the Department of the Environment and Housing. These data were captured using the Thematic Mapper (TM) sensor from the Landsat satellite. The characterization of landscape patterns was focused on the Catalan Pyrenees.
The aim was to group those protected natural areas present in the same geographical unit (in this case the Catalan Pyrenees) to more effectively and realistically compare the landscape patterns experienced in the last fifteen years with those of the study area.
A total of thirteen variables (landscape indices) were used to assess landscape patterns and changes in landscape composition and configuration between forested and non-forested areas in the seven protected natural areas every five years from 1987 to 2002. Those variables were [15] see Table A1: Another aspect that must be taken into account in this analysis is the conjunction of two well-differentiated biogeographical regions throughout the Pyrenees: The Alpine and the Mediterranean (see Figure 4). To assess whether or not these patterns are of a purely biogeographical nature, we first conducted a multivariate analysis to simplify or reduce the dimension-a principal component analysis (PCA) on a total of thirteen variables (landscape indices) used to assess landscape patterns and changes in landscape composition and configuration between forested, transitional and non-forested areas in the seven protected natural spaces for the following years 1987, 1992, 2002 and 2017. To assess whether the landscape patterns displayed in the protected natural spaces included in the study are purely biogeographical in nature, I first applied -a multivariate method of simplification or size reduction -a principal component analysis (PCA). PCA algorithm have been used to compress the dataset (13 landscape metrics) and a total of four components were captured for each of the four periods considered, comprising 95%-98% of the total variance presented by the thirteen indices used to characterize. The KMO (Kaiser-Meyer-Olkin) descriptors and the Barlett sphericity test were also calculated, both presenting reasonable values for accepting the analysis. Then, a hierarchical cluster analysis (CA) was applied creating two groups (Alpine and Mediterranean) between protected natural areas based on landscape metrics (see Figure 4). Clustering allows us to identify which landscape patterns are alike and potentially categorize them therein. K-means clustering is the simplest and the most used clustering method for splitting a dataset into a set of k groups.
Finally, linear combinations of variables, known as discriminant functions, of the dependent variables that maximize the separation between the groups are used to identify the relative contribution of the p variables that best predict group membership. In our case, to know which landscape metrics separated the two groups in the best way, a discriminant analysis (DA) was carried out.
As can be seen in the Table 1, the discriminant analysis allows us to distinguish between those landscape metrics that are significant over the four time intervals (1987, 1992, 2002 and 2017) and those that are only significant for certain time intervals when forming groups. To assess whether the landscape patterns displayed in the protected natural spaces included in the study are purely biogeographical in nature, I first applied -a multivariate method of simplification or size reduction -a principal component analysis (PCA). PCA algorithm have been used to compress the dataset (13 landscape metrics) and a total of four components were captured for each of the four periods considered, comprising 95%-98% of the total variance presented by the thirteen indices used to characterize. The KMO (Kaiser-Meyer-Olkin) descriptors and the Barlett sphericity test were also calculated, both presenting reasonable values for accepting the analysis. Then, a hierarchical cluster analysis (CA) was applied creating two groups (Alpine and Mediterranean) between protected natural areas based on landscape metrics (see Figure 4). Clustering allows us to identify which landscape patterns are alike and potentially categorize them therein. K-means clustering is the simplest and the most used clustering method for splitting a dataset into a set of k groups.
Finally, linear combinations of variables, known as discriminant functions, of the dependent variables that maximize the separation between the groups are used to identify the relative contribution of the p variables that best predict group membership. In our case, to know which landscape metrics separated the two groups in the best way, a discriminant analysis (DA) was carried out.
As can be seen in the Table 1, the discriminant analysis allows us to distinguish between those landscape metrics that are significant over the four time intervals (1987, 1992, 2002 and 2017) and those that are only significant for certain time intervals when forming groups.

Heuristic Models
Managers and researchers have focused on locating spaces that contain a particular endangered species, a particular habitat of interest, the whole of an ecosystem, areas with high scenic and lucrative values or any other type of natural resource. Heuristics are concerned with providing exploratory methods or algorithms during problem-solving, solutions being discovered by evaluating the progress made in the search for a final result. In mathematical programming, a heuristic method is a procedure that determines good or near-optimal solution to an optimization problem. Heuristic methods are used to speed up the process of finding a good enough solution.
In the study of nature reserve design there is a long-standing debate on whether it is better to create a single large reserve or many smaller ones [38,39]. This problem is referred to in the literature as SLOSS (single large or several small). Two heuristic models have been designed to provide a solution to this problem [40]. One is the so-called Mentor (stepping-stone strategy), which aims to allocate passerine habitats among existing reserves. The basis of this model is to establish a threshold in reference to the cost of the distance between existing reserves; when this cost of distance exceeds the threshold, then the model will create new passerine habitats to connect the reserves. The other model is called the Enlarge strategy and focuses on analyzing the suitability of the zones surrounding the existing reserve and looking for the most suitable zone to increase its surface area.
Land use optimization, a multifaceted process, requires complex decisions, including selection of land uses, forecasting land use allocation percentage and assigning locations to land uses [41,42]. One of the more recently developed heuristic methods [43,44] was chosen for the potential creation of open spaces, as it can be used to determine which locations are most feasible for selection based on a suitability matrix. This process can be used to define hundreds of possible places that meet certain purposes, requirements or criteria. The method is known as PGP (Patch-Growing Process) and is characterized by being able to generate representative fragments from a seed cell. It tends to select neighboring cells with higher aptitude values and takes into account parameters related to their compactness throughout the process of building a fragment. Despite great progress in these techniques, heuristic methods cannot guarantee the optimal solution to the problem. In addition to this, they are also unable to calculate the distance of the proposed result relative to what remains at the final optimum [43,45,46].
The main aim of this research is to develop a model that is capable of searching for all feasible sites and then select from the resulting set the one that optimizes a set of criteria (PLAND, PD, MPS, GYRATE, ENN), see Table A1. Each patch begins from an initial cell and continues to accrete neighboring cells until it reaches a desired total habitat value, see Figure 5. Patch-Growing Process is simply defined as the adhesion of neighboring cells to a preestablished or random seed cell, until the desired total habitat value and determined compactness are obtained. This method of identifying/generating feasible patches, in addition to focusing on the values of the input raster, also includes the parameters of size, shape and compactness, which can be manipulated by the user. The Patch-Growing Process [44] allows users to generate patches using a set of rules or a finite set of instructions to solve the problem. The result of applying the model (PGP) is a set of patches that meet the requirements of suitability, shape, size and density. In the case of this research, the aim of using PGP was to determine the maximum number of patches that should be preserved as open space.
The process (PGP) consists of six phases: • Defining the parameters: in this first section, users enter parameters such as the suitability value for the seed cell, the shape of the initial neighborhood matrix (W), the importance of compactness (N), acceptance percentage, size and number of patches. These parameters are described in detail in the next section.

•
Selecting a seed cell that meets the previously selected suitability value. In the event that the cell is close to the edge of the study area, the model will choose another one farther from the edge so that the generated patch does not extend beyond the study area.

•
Depending on the selected radius of activity (neighborhood matrix), the next step will be to adhere neighboring cells to the seed cell.

•
Calculating the final suitability value (CSi) of the adjoining neighboring cells: where Si = the suitability value for cell i, ei = the number of edges that cell i shares with the existing patch, and N = the weight attached to sharing edges with the existing patch.
• Arrangement of neighboring cells and selection: values are ordered from highest to lowest and the cells that will take part in the second iteration are selected from the percentage of acceptance X. In Figure 5, the final suitability values (CSi) for the four neighboring cells are 1, 2, 4 and 3. With X = 0.5, then 50% of the cells will be selected in order of more to less suitability. In the example, this would be cells 4 and 3. These cells adhere to the seed cell and the previous steps are repeated. Neighboring cells are then sought and the final suitability values (CSi) for each cell are recalculated. As Figure 5 shows, the new neighboring cells have the values 2, 3, 2, 3, 4 (the margin is displayed with dashed lines because, unlike the others, this cell shares two edges), 1 and 4. The cells to be incorporated in this case will be 3, 3, 4 and 4. • This process will be repeated until the area or minimum suitability value of the patch to be generated are complete. Patch-Growing Process is simply defined as the adhesion of neighboring cells to a pre-established or random seed cell, until the desired total habitat value and determined compactness are obtained. This method of identifying/generating feasible patches, in addition to focusing on the values of the input raster, also includes the parameters of size, shape and compactness, which can be manipulated by the user. The Patch-Growing Process [44] allows users to generate patches using a set of rules or a finite set of instructions to solve the problem. The result of applying the model (PGP) is a set of patches that meet the requirements of suitability, shape, size and density. In the case of this research, the aim of using PGP was to determine the maximum number of patches that should be preserved as open space.
The process (PGP) consists of six phases: • Defining the parameters: in this first section, users enter parameters such as the suitability value for the seed cell, the shape of the initial neighborhood matrix (W), the importance of compactness (N), acceptance percentage, size and number of patches. These parameters are described in detail in the next section.

•
Selecting a seed cell that meets the previously selected suitability value. In the event that the cell is close to the edge of the study area, the model will choose another one farther from the edge so that the generated patch does not extend beyond the study area.

•
Depending on the selected radius of activity (neighborhood matrix), the next step will be to adhere neighboring cells to the seed cell.

•
Calculating the final suitability value (CSi) of the adjoining neighboring cells: where Si = the suitability value for cell i, ei = the number of edges that cell i shares with the existing patch, and N = the weight attached to sharing edges with the existing patch.
• Arrangement of neighboring cells and selection: values are ordered from highest to lowest and the cells that will take part in the second iteration are selected from the percentage of acceptance X. In Figure 5, the final suitability values (CSi) for the four neighboring cells are 1, 2, 4 and 3. With X = 0.5, then 50% of the cells will be selected in order of more to less suitability. In the example, this would be cells 4 and 3. These cells adhere to the seed cell and the previous steps are repeated. Neighboring cells are then sought and the final suitability values (CSi) for each cell are recalculated. As Figure 5 shows, the new neighboring cells have the values 2, 3, 2, 3, 4 (the margin is displayed with dashed lines because, unlike the others, this cell shares two edges), 1 and 4. The cells to be incorporated in this case will be 3, 3, 4 and 4. • This process will be repeated until the area or minimum suitability value of the patch to be generated are complete.

Evaluation of the Open Space Fragmentation Process
A set of landscape indices were applied to the two land use and land cover maps from 1957 and 2017. This provided information on the composition (see Table 2) and configuration of the landscape, while allowing us to assess changes to the landscape over the last fifty years and note the fragmentation of non-forest habitats [35]. In 1957, there was a balance between forest and non-forest habitats, even though there had been quite a lot of deforestation as a result of high levels of charcoal-related activity. Approximately 48% of the 1957 landscape corresponded to forest habitats and pasture and around 52% non-forest habitats [4,5,9]. While in 2017, non-forest habitats have been reduced to 17%.

Model Parameters
A total of four components were captured for each of the four periods considered, comprising 95-98% of the total variance presented by the eleven indices used to characterize landscape patterns. The KMO (Kaiser-Meyer-Olkin) descriptors and the Barlett sphericity test were also calculated, both presenting reasonable values for accepting the analysis.
Following this, after capturing the principal components, a second analysis was performed, in this case a cluster analysis (CA). From these four heterogeneous components, natural spaces were grouped into two hypothetical groups (Alpine region and Mediterranean region).
The cluster analysis was needed to verify whether different behaviors occur in the landscape structure between the different protected natural spaces (PNS) depending on the biogeographic region in which they are located. This being the case, as evidenced in the Figure 6, it is very important to consider only the PNS data from the biogeographic region of the study area when establishing the input parameters for the model, so that the landscape structure resulting from executing the PGP approximates the reality of these PNS as much as possible.
Once the protected natural spaces had been characterized according to their landscape pattern, we could establish the necessary input parameters to indicate the landscape structure we wished to obtain in the PGP (Patch-Growing Process). These parameters will vary depending on the biogeographical region in which the landscape pattern is to be improved through the generation of new patches. In this specific case, as the study area was within the Mediterranean region (see Table 3), I used the specific parameters for this region, as shown in the following table: Once the protected natural spaces had been characterized according to their landscape pattern, we could establish the necessary input parameters to indicate the landscape structure we wished to obtain in the PGP (Patch-Growing Process). These parameters will vary depending on the biogeographical region in which the landscape pattern is to be improved through the generation of new patches. In this specific case, as the study area was within the Mediterranean region (see Table  3), I used the specific parameters for this region, as shown in the following table: Table 3. Configuration analysis for Mediterranean mountain regions. S (patches with an average area of 1 ha), size M (patches with an average area of 5 ha) and size L (patches with an area average of 50 ha).

Heuristic Models (Patch-Growing Process)
The following map (Figure 7) shows the patches generated using the PGP application after entering the input parameters (which were calculated based on the size of the patches and the biogeographic region) and using the "Neutral" suitability weights. This map shows those new spaces generated using PGP, which correspond to the areas with the highest level of suitability in the neutral matrix. If we add the agricultural areas that remained active in 2017, we obtain the following map:

Heuristic Models (Patch-Growing Process)
The following map (Figure 7) shows the patches generated using the PGP application after entering the input parameters (which were calculated based on the size of the patches and the biogeographic region) and using the "Neutral" suitability weights. This map shows those new spaces generated using PGP, which correspond to the areas with the highest level of suitability in the neutral matrix. If we add the agricultural areas that remained active in 2017, we obtain the following map: Figure 8 represents the total surface area of agricultural lands (20%) that could be maintained, recovered and generated from the neutral suitability matrix and the landscape configuration parameters for the protected natural spaces present in the Mediterranean region.    Table 4 shows the variables (landscape indices) that best represent the changes in the agrarian lands scenario proposal generated by means of PGP with respect to the situation in 2017. The PGP proposal reveals an improvement in patch density (PD), mean patch size (MPS) and edge margin density (ED). The three variables have increased in comparison with the situation in 2017. This means a greater number of open space patches distributed throughout the study area, with larger sizes than the current ones (2017) and the recovery of edge habitats, which were becoming extinct as a result of the predominant single patch of woodland that comprises the vast majority of the study area.     Table 4 shows the variables (landscape indices) that best represent the changes in the agrarian lands scenario proposal generated by means of PGP with respect to the situation in 2017. The PGP proposal reveals an improvement in patch density (PD), mean patch size (MPS) and edge margin density (ED). The three variables have increased in comparison with the situation in 2017. This means a greater number of open space patches distributed throughout the study area, with larger sizes than the current ones (2017) and the recovery of edge habitats, which were becoming extinct as a result of the predominant single patch of woodland that comprises the vast majority of the study area.  Table 4 shows the variables (landscape indices) that best represent the changes in the agrarian lands scenario proposal generated by means of PGP with respect to the situation in 2017. The PGP proposal reveals an improvement in patch density (PD), mean patch size (MPS) and edge margin density (ED). The three variables have increased in comparison with the situation in 2017. This means a greater number of open space patches distributed throughout the study area, with larger sizes than the current ones (2017) and the recovery of edge habitats, which were becoming extinct as a result of the predominant single patch of woodland that comprises the vast majority of the study area. In the case of the GYRATE variable, has increased compared to 2003, improving the compactness of the patches. In contrast, the value of the PAFRAC variable has decreased compared to the situation in 2017, which means that there is a trend towards regular-shaped patches. The following two variables provide significant improvements in configuration of the landscapes, such as a considerable reduction in the distance between the patches (ENN) compared to 2017, which guarantees a cohesion (COHESION) between them, as reflected in the increase in this value.

Impacts and Improvements of the Model Proposed Using PGP (Patch-Growing Process) on the Current Landscape Structure
However, I also observe that the improvements experienced through the use of PGP remain very far from the situation in 1957. This finding was anticipated and is explained by the overexploitation of natural resources that existed in the study area during the period, with half of the study area being open space compared to the 13% it represents nowadays.
Finally, I have the analysis corresponding to landscape, see Table 5. The percentage of surface area that corresponded to open spaces in 1957 was 47%, compared to the 13% that these environments currently represent. The use of PGP provides for a 21% improvement in total open space. Another very important improvement on 2017 presented by the scenario created by PGP from the point of view of diversity is the increase in landscape heterogeneity, as is clearly reflected in the SHDI and SHEI variables. And a further improvement produced by the patches generated by PGP is the breaking up (disintegration) of space with respect to uniformity and continuity (due to the forest mass), as observed in the reduction of the CONTAG variable [24,47].

Discussions and Conclusions
The Patch-Growing Process (PGP) model has enabled us to use a raster map to create patches of different sizes and suitability values with certain standards regarding shape and compactness. These patches represent the open spaces that should preferably be maintained and/or recovered. The results obtained from the PGP modeling would mean an increase in the proportion of surface area currently occupied by open spaces and a more heterogeneous landscape structure. This increase in open space would be around 8% (500 ha) of total open space in the study area. Such an increase in absolute data also represents a significant improvement in the landscape structure indices, including a reduction in the isolation of fragments and a recovery of the edge habitats thanks to the increase in surface area. In short, I must emphasize that this represents a recovery of landscape heterogeneity to counter the uniformity of the current landscape within the study area.
Many inhabitants have the feeling of being enclosed within a 'wild' landscape that they no longer recognize and that they cannot control any longer. The value for maintaining traditional landscapes was highlighted together with the landscape's importance for crop and livestock production, such as aesthetic and identity values, indicating the need for a careful balance between forest a non-forest habitats.
Analyzing perceptions and preferences regarding the demand for ecosystem services in Mediterranean mountain areas can help to visualize this contrast and help to adapt agricultural and management practices, to local users' and society's needs. Local assessments of ecosystem services support an understanding of the multi-functional and multi-beneficial character of agricultural landscapes and may indicate potential losses in ecosystem services caused by changing agricultural areas.
The modeling presented can be extrapolated, with the relevant adaptations, to other territories with a similar dynamic of agricultural abandonment and progressive forestry. Interventions at the territorial level and the management of these open spaces must become a basic priority for improving landscape and biological production and diversity when it comes to natural spaces, as well as reducing fire risk and improving landscape perceptions. Landscapes affected by the abandonment of traditional activities contribute to the appearance of woody species and the accumulation of dead plant matter. Abandoned land therefore becomes more vulnerable to fire, which is also likely to spread more easily.
Policymaker may come to think of this process of abandonment as an opportunity to restore wildlife in protected natural areas. Defenders of rewilding can argue that traditional agriculture is not environmentally friendly and that rewilding is beneficial for ecosystems and humanity. Furthermore, this type of passive management allows for much lower maintenance costs than other management options. However, it has been largely disregarded as a management option because the uncontrolled growth of the forest raises problems involving an increasing risk of wildfire.
The difficulty of recovering this type of agricultural lands in protected natural areas due to the restrictive conservation regulations, makes this research a necessity to rethink the regulations to the current needs of these territories. Also, this research has highlighted a current need to reclaim agricultural areas to facilitate the little livestock and subsistence farming that takes place in these protected natural areas. Models have been designed to indicate the areas in which the restoration of open areas will be the most advisable and the most beneficial, considering environmental, social and economic factors.
Agrarian habitat patches have created higher landscape diversity, which is now being lost as the forest reclaims the abandoned pastures. The sustainable development of this territory would have to make compatible objectives of conservation of biodiversity and the preservation of their Mediterranean features with support to agricultural activities.

Contagion (CONTAG)
Equals minus the sum of the proportional abundance of each patch type multiplied by the proportion of adjacencies between cells of that patch type and another patch type, multiplied by the logarithm of the same quantity, summed over each unique adjacency type and each patch type; divided by 2 times the logarithm of the number of patch types; multiplied by 100 (to convert to a percentage).
g ik · ln(P i )

Shannon's Diversity Index (SHDI)
Equals minus the sum, across all patch types, of the proportional abundance of each patch type multiplied by that proportion. Note, Pi is based on total landscape area (A) excluding any internal background present.

Shannon's Evenness Index (SHEI)
Equals minus the sum, across all patch types, of the proportional abundance of each patch type multiplied by that proportion, divided by the logarithm of the number of patch types. In other words, the observed Shannon's Diversity Index divided by the maximum Shannon's Diversity Index for that number of patch types. Note, Pi is based on total landscape area (A) excluding any internal background present. Table A1. Cont.

Patch Richness Density (PRD)
Equals the number of different patch types present within the landscape boundary divided by total landscape area (m2), multiplied by 10,000 and 100 (to convert to 100 hectares). Note, total landscape area (A) includes any internal background present. PRD = m A (10000)(100)

Patch Density (PD)
Equals the number of patches of the corresponding patch type divided by total landscape area (m2), multiplied by 10,000 and 100 (to convert to 100 hectares). Note, total landscape area (A) includes any internal background present. PD = n i A (10000)(100)

Radius of Gyration (GYRATE)
Equals the mean distance (m) between each cell in the patch and the patch centroid. GYRATE = z r h ijr Z

Fractal Dimension (FRAC)
Equals 2 times the logarithm of patch perimeter (m) divided by the logarithm of patch area (m2); the perimeter is adjusted to correct for the raster bias in perimeter. FRACT = 2 ln p ij ln a ij

Euclidean Nearest Neighbor Distance (ENN)
Equals the distance (m) to the nearest neighboring patch of the same type, based on shortest edge-to-edge distance. Note that the edge-to-edge distances are from cell center to cell center.