Assessing Land Use-Cover Changes and Modelling Change Scenarios in Two Mountain Spanish National Parks

Land Use-Cover Changes (LUCCs) are one of the main problems for the preservation of biodiversity. Protected Areas (PAs) do not escape this threat. Some processes, such as intensive recreational use, forest fires or the expansion of artificial areas taking place inside and around them in response to their appeal, question their environmental sustainability and their efficiency. In this paper, we analyze the LUCCs that took place between 1990 and 2006 in two National Parks (NPs) belonging to the Spanish network and in their surroundings: Ordesa and Monte Perdido (Ordesa NP) and Sierra de Guadarrama (Guadarrama NP). We also simulate land use changes between 2006 and 2030 by means of Artificial Neural Networks (ANNs), taking into account two scenarios: trend and green. Finally, we perform a multi-temporal analysis of natural habitat fragmentation in each NP. The results show that the NPs analyzed are well-preserved and have seen hardly any significant LUCCs inside them. However, Socioeconomic Influence Zones (SIZs) and buffers are subject to different dynamics. In the SIZ and buffer of the Ordesa NP, there has been an expansion of built-up areas (annual rate of change = +1.19) around small urban hubs and ski resorts. There has also been a gradual recovery of natural areas, which had been interrupted by forest fires. The invasion of sub-alpine grasslands by shrubs is clear (+2735 ha). The SIZ and buffer of the Guadarrama NP are subject to urban sprawl in forest areas and to the construction of road infrastructures (+5549 ha and an annual rate of change = +1.20). Industrial area has multiplied by 3.3 in 20 years. The consequences are an increase in the Wildland-Urban Interface (WUI), greater risk of forest fires and greater fragmentation of natural habitats (+0.04 in SIZ). In the change scenarios, if conditions change as expected, the specific threats facing each NP can be expected to increase. There are substantial differences between the scenarios depending on whether or not incentives are accepted and legal restrictions are respected.


Introduction
Protected Areas (PAs) are a key for mitigating climate change, preserving biodiversity, providing ecosystem services and fostering human wellbeing.The declaration of PAs has increased globally because of increased environmental sensitivity [1]; to counter the threats of climate change [2]; land use changes [3]; deforestation [4]; the risk of flooding [5]; the risk of forest fires [6]; habitat fragmentation [7]; the propagation of invasive species [8]; urban pressure [9]; and recreational use [10].
In 1990, PAs covered 8.6% of the Earth's surface.According to the World Data Base of PAs [11], they have expanded from 84,577 in 2003 to 217,155 in 2016.93% of them occupy 19.8 million km 2 , equivalent to 14.7% of the worldwide surface area of terrestrial ecosystems and continental and inland waters, excluding the Antarctic.The remainder are Protected Marine Areas, which cover 14.9 million km 2 , 4.12% of the world's oceans and 10.2% of marine and coastal waters under national jurisdiction [11][12][13].To reach Aichi Goal 11, the Convention on Biological Diversity recommends that by 2020, at least 17% of terrestrial areas and continental waters be recommended, as well as 10% of coastal and marine areas [14].
In Europe, PAs occupy 13.6% of the land surface and of continental waters [15].In Spain, from 1990-2013, the number of PAs multiplied by seven, and their surface area tripled [16].In the worldwide and continental context, Spain plays a relevant role in the preservation of biological diversity.Today, more than 27% of the surface area occupied by terrestrial ecosystems and their inland waters is protected under national, European or worldwide networks.It is the EU country that contributes most territory to the Natura 2000 network.
Amongst the categories defined by the International Union for Conservation of Nature (IUCN), National Parks (NPs) are the figure that is most widely used for PAs of high natural value [17].In Spain, NPs are covered by a law enacted 100 years ago [18] and recast recently [19].Since 1918, when the first two national parks were declared, the network has grown at an average annual rate of 23%.Such parks are a key in various Spanish strategies for nature protection (for example, the Spanish Strategic Plan for the Natural Heritage and Biodiversity 2011-2017).
For decades, the scientific community has been showing interest in the spatial and temporal analysis of Land Use-Cover Changes (LUCCs) [20][21][22].Cartography of LUCCs is crucial for monitoring ecosystems at different scales [23][24][25] and for assessing the impact of changes on PAs and biodiversity.Their study allows one to know the size, extent, type and trends of LUCCs and to identify the main factors of change.This is a fundamental previous step for the design of conservation policies in PAs [3], especially in NPs [26].Such an analysis is also required for assessing the efficiency of PAs [27].
Moreover, increasing importance is being placed by scientists and managers on the LUCCs taking place inside PAs and in their surroundings [28].External pressure from such transformations has different impacts on the biodiversity of PAs and their Socioeconomic Influence Zones (SIZs), reduces their efficiency and amounts to a threat for environmental sustainability [29,30].
Some aggressive LUCCs-expansion of intensive farming lands and of built areas-lead, amongst others, to the fragmentation of natural habitats [31], the loss of connectivity for habitats of ecological interest [32], a reduction in the floodable area of protected wetlands and poorer water quality in such wetlands [33].
At a local level, in the Pyrenees of the province of Huesca, rising temperatures are leading to invasion by woody vegetation and displacement of forests to higher altitudes, with a decline in grasslands [34].The reduction in traditional farming activities has been an explanatory factor for recolonization by woody species [35].The invasion of abandoned farmland and of grassland by more xerophilous woody species [36] increases the risk of forest fires.
In the Sierra de Guadarrama, close to the metropolitan area of Madrid, there is an increased recreational pressure on the PAs [37], which are attractive for their landscape and their environments.There is also great pressure from urban sprawl on ecosystems of high ecological value [38][39][40].
In order to help in the decision-making of the stakeholders, the scientific community is working on the simulation of various change scenarios using various methods and tools: logistic regression [41], neural networks [42,43], cellular automata [44][45][46] and, sometimes, combined methods [47,48] or compared [49,50].
The main objective of this research is to analyze the changes in land use that occurred, between 1990 and 2006, in two mountain Spanish National Parks and in their surroundings.In addition, considering these trends of change in recent decades, we intend to simulate different scenarios of change-trend and green-that are expected between 2006 and 2030.For this, we use a simulator based on Artificial Neural Networks (ANNs).In addition, we have analyzed the changes in the fragmentation of natural and semi-natural habitats recorded in the last few decades, as well as those that are expected if the simulated changes are met.We compare the processes that have occurred inside the two NPs with those that have occurred outside of them, in their surroundings, in order to find substantial differences relating to environmental sustainability.The ultimate aim of the research is to provide information that will be of use to the managers of affected NPs and to local administrators in their preventive environmental and territorial decisions.We also aim to increase awareness among those responsible for the Spanish network so that similar methods can be applied in other NPs.

Study Areas
Mountain ecosystems are the focus of interest in this work.For this reason, we have chosen two mountain NPs that, besides belonging to two different biogeographic regions, had contrasting geographical and socioeconomic characteristics.For this study, we selected a sample of two mountainous NPs belonging to the Spanish network: one of the oldest (Ordesa NP) and the newest (Guadarrama NP); see Figure 1.
Ordesa NP is located in a rural area with poor accessibility, along the frontier with France and far from the main urban centers in Spain.It was declared in 1918 to be representative of high mountain ecosystems and of systems linked to erosion formations and sedimentary rocks in the Pyrenee Province of the EuroSiberian Region [51].Guadarrama NP is the latest to be declared, in 2013, although several groups have shown interest in its declaration as a national park since the beginning of the 20th century.It is representative of mountain ecosystems of the Mediterranean biogeographical region.It is in a peri-urban zone, less than 50 km from Madrid and its metropolitan area, which has a population of about 6.5 million inhabitants [52].It is very close to motorways and high-speed railway lines.It is very accessible and is visited by 3,000,000 visitors a year [53].We intend to know if the processes of change of land uses are significantly different between these two NPs.
In addition, we also selected various areas surrounding them: their Peripheral Protection Zones (PPZs), Socioeconomic Influence Zones (SIZs) and external 5 km-wide buffers plus the extensions of the buffers up to the administrative limits of the municipalities that are totally or partially affected by them.In Table 1, we summarize the main characteristics of the study areas.

Materials
Our research follows the workflow shown in Figure 2. We used ArcGIS v10.3 (ESRI Inc., Redlands, CA, USA) for vector processing of the downloaded data and to draw up the buffer.For LUCC analyses, we used IDRISI-TerrSet (v.18.31,Clark Labs, Clark University, Worcester, MA, USA) [56] and Land Change Modeller (LCM) (Clark Labs, Clark University, Worcester, MA, USA) for designing the scenarios.For the calibration and validation of land use change scenarios, we used the Map Comparison Kit v3.2.3 software (RIKS BV, Maastricht, The Netherlands) [57].Finally, we used GUIDOS-MSPA (v.2.4, rev.2, Joint Research Centre, European Commission, Ispra, Italy) [58] to analyze the spatial landscape pattern.We used three datasets developed by the CORINE Land Cover (CLC) Project [59]: years 1990, 2000 and 2006.

Materials
Our research follows the workflow shown in Figure 2. We used ArcGIS v10.3 (ESRI Inc., Redlands, CA, USA) for vector processing of the downloaded data and to draw up the buffer.For LUCC analyses, we used IDRISI-TerrSet (v.18.31,Clark Labs, Clark University, Worcester, MA, USA) [56] and Land Change Modeller (LCM) (Clark Labs, Clark University, Worcester, MA, USA) for designing the scenarios.For the calibration and validation of land use change scenarios, we used the Map Comparison Kit v3.2.3 software (RIKS BV, Maastricht, The Netherlands) [57].Finally, we used GUIDOS-MSPA (v.2.4, rev.2, Joint Research Centre, European Commission, Ispra, Italy) [58] to analyze the spatial landscape pattern.We

Land Use Changes between 1990 and 2006
Firstly, although the size of buffers around PAs has been established in various ways in the literature [3,26,60], in this study, we chose to design a 5-km buffer, which suits the characteristics of the two study areas and is also in line with other previous studies [28,31].A buffer with this size increases the probability of the territories included in these "control areas" having similar geographical characteristics to the "cases" (NPs and their SIZs).On the contrary, if the buffer is too small, the probability of finding differences in land uses-land cover, inside and outside NP, would be small.
Second, we drew up an initial cross-tabulation between the CLC vector maps to identify unusual and unexpected land use changes (Figure A1).Some authors find interpretation and location errors [61,62].After checking any doubts with the aerial orthophotos provided under the National Plan for Aerial Orthophotography [63], any errors found were corrected.Subsequently, we transformed the CLC vector maps to a raster format with a 50 × 50 m pixel size.From the CLC legend, we made two different groupings.The first is a simplification of Level 3 into six categories: Urban areas (URB), Industrial areas (IND), arable land and permanent crops (Agricultural (AGR)), Heterogeneous agricultural areas (HET), Forests (FOR) and Shrubs and Grasslands (SHR-GRAS).We considered the OTH category (Other: open spaces with little vegetation, wetlands and water bodies) to be stable.Because of the singularity of the Ordesa NP, the CLC Level 3 was simplified by dividing it into six categories: Urban areas (URB), Agricultural areas (AGR), Grasslands (GRAS), Shrubs (SHR), Forests (FOR) and Others (OTH: open spaces with little vegetation, wetlands and water bodies) (Table A1).We used this grouping to simulate future scenarios.
Third, we drew up cross-tabulation matrices [42] to obtain values and maps of changes between 1990 and 2006.We calculated the annual rate of change of each use [43].We then compared the results with the PAs and with the 5-km buffer.The aim was to find some of the main processes of land use change that had already taken place: built-up land, naturalization and disturbances and exchanges in natural areas [3,64].

Simulated Scenarios of Change in 2030
Firstly, using LCM, we simulated land use in 2030 under different scenarios.In the Ordesa NP, we considered two different scenarios: (a) trend scenario and (b) green scenario.The trend or "business as usual" scenario shows what would happen if the past trend in 1990-2000-2006 were to continue until 2030.The green scenario shows what would happen if there were more active reforestation policies and if greater importance were placed on the natural environment.In both scenarios, we considered two alternatives.Under the first, we considered that land planners will respect the restrictions imposed by law in certain zones (NPs and other PAs, public utility forests or the public domain zones of rivers, roads and railway lines) (Table A2).We also considered the incentives that planners are promoting in certain zones to encourage the naturalization of land for various purposes (for example, to expand habitats of community importance).Under the second alternative, we did not take these restrictions or incentives into account to avoid distorting the simulation, for example excessive preparation for use changes in the location or the elimination of other conditioning factors [65].We also considered that the stakeholders acting in the territory might not respect legal restrictions on certain land uses or might not take incentives into account.Although such scenarios would be illegal and would not be a sustainable alternative, examples are often found.In the Guadarrama NP, we only took into account a trend scenario, both with and without restrictions and incentives.
In the model simulated with LCM, we related land use and driving factors (Table 2) by means of a Multi-Layer Perceptron neural network (MLP).The MLP classifier uses a Backward Propagation algorithm (BP), one of the most used neural networks.We have taken into account a variable number of input and output nodes depending on the designed scenarios (between 7 (in Ordesa NP) and 28 (in Guadarrama NP)).For each category, the number of pixels per class is randomly divided between training and testing: half of the sample size for training and half for testing.The first are used in the analysis and the second to validate the results.Samples used for the training process are taken from pixels that have and have not undergone the transition being modeled.We selected an automatic and dynamic training to get the proper weights both for the connection between the input and hidden layer and between the hidden and the output layer for the classification of the unknown pixels.We ran the model 13 times, changing the learning rate parameters after obtaining an accuracy rate above 50% for each transition.We selected the driving factors by consulting experts and prior studies.We used Cramer's V statistic [66] to test the explanatory power of each variable and selected the most relevant.
We did carry out calibration in order to improve the results in the Ordesa NP.Taking the sequence of maps for 1990-2000 as a base, we simulated a land use model in 2006 and compared it with the real map for 2006.The calibration compared the number of pixels and the spatial location of each land use-land cover with that simulated.Models were tuned by minimizing the disagreement between actual and simulated maps.This was done changing the matrix of transition probabilities, selecting the driving factors and adding or changing the size and/or weight of the neighborhoods.
We used the Kappa simulation statistic (Kappa Sim), which tests only the changed areas of the map [67][68][69].In addition, we have used TransLoc and Transition, which evaluate the accuracy in the location and quantity, respectively, of the pixels that experience a land use change.
Second, we drew up cross-tabulation matrices [42] to obtain values and maps of changes between 2006 and 2030.We calculated the annual rate of change of each use [43].We then compared the results with the PAs and with the 5-km buffer.The aim was to find some of the main processes of land use change that had already taken place and that could be expected in different scenarios: built-up land, naturalization and disturbances and exchanges in natural areas.

Fragmentation of Habitats of Interest
In order to find the dynamics of landscape structure, we took into account Level 1 of the CLC legend.We reclassified the maps in binary format.In the case of Guadarrama, we considered Class 1 (artificial areas) as the background and combined Classes 2, 3, 4 and 5 in a single target category (agricultural and natural areas) linked to the habitats represented in the National Park and the surroundings.In the case of Ordesa, we considered Class 1 (shrubs) as the background and only Class 2 (herbaceous vegetation) as the target category.The managers of this National Park and of other PAs around it are worried about the invasion of alpine and sub-alpine grasslands (focal habitat) by shrubs resulting from the reduction in extensive cattle-breeding and from climate change.
We calculated a Habitat Fragmentation Index (HFI) [6], in our case the sum of agricultural and natural habitats in NP and in PPZ, SIZ and the corresponding buffer.This goes, with continuous values, from 1 (greatest fragmentation) to 2 (least fragmentation).It assigns a different weight to each of the entities mapped in terms of relations between resilience and spatial coherence [70,71].There is a constant gradation from the core (greatest weight) to the islets (least weight).The index relates the number of pixels in each category or fragmentation entity to their weights.In the Ordesa NP, we used the same method to calculate a fragmentation rate for the alpine and sub-alpine grasslands, which, in this case, are habitats of interest.We studied changes expected between 2006 and 2030, considering the green scenario (GS), without incentives or restrictions (GS30), in order to show changes in the landscape if the PA regulatory instruments are not respected.

Land Use Changes between 1990 and 2006
In general and as expected, there were few changes in the Ordesa NP.This is a rural district in which land uses are strictly regulated because they belong to an NP, two RPs and other sites within the Natura 2000 Network (1.1937).Global persistence is over 98% in the study area (Table 3).The urban areas underwent the greatest annual rate of change (1.1937), and the classes with the greatest net changes were grasslands (about −2800 ha) and shrubs (+2463 ha). Figure 3 shows the spatial distribution of LUCCs taking place between 1990 and 2006.In the NP and in its PPZ, no changes were recorded.Most of them were located inside the buffer and some in its SIZ.The use changes correspond to three different processes: natural regeneration, degradation of the vegetation and increase in artificial surfaces.The most extensive changes correspond to a progression of the vegetation to the climax, from grasslands to shrubs to forest (in green colors).They can be seen in detail in Window C of Figure 4.There are two possible causes for such progressions.The first is the abandonment of agricultural land and, above all, of grasslands because of the abandonment of traditional cattle-farming.Population density is low (<5 inhabitants/km 2 ) and the rate of ageing high (>26%).From an ecological point of view, this could be considered positive.However, it is worrying for land managers because sub-alpine grasslands are focal habitats in the NP.In addition, this transition involves an increase in available potential fuel and, consequently, an increase in the risk of forest fires in an area of high ecological value.According to the CLC90-06 series and at the landscape scale, there was no invasion of grasslands by shrubs within the NP.Secondly, the progression of vegetation is the consequence of recovery from previous forest fires.
Changes associated with degradation of vegetation (yellow colors) stem from more recent forest fires and can be seen in greater detail in Window A of Figure 4 (transitions from forest or shrubs to grassland).Finally, the expansion of urban zones and infrastructure (red colors) mostly took place in the surroundings of Sabiñánigo, the center of the district with the largest population and industries.Window B of Figure 4 shows a golf course under construction.In the center of Figure 3, certain long patches can be seen, which have changed to artificial use as a result of the widening of the N-260 motorway.In the NW quadrant of the same figure, there are also small grasslands that have been converted into new ski slopes and other infrastructure associated with skiing in Formigal.
In the Guadarrama NP and its surroundings, most of the changes are concentrated in the south, in the region of Madrid and, especially, in the buffer and its SIZ (Figure 5).This is the most dynamic region from the demographic and socio-economic points of view, so it is the region that exerts the greatest pressure on the environment.
According to the CLC90-06 series and at the landscape scale, there was no invasion of grasslands by shrubs within the NP.Secondly, the progression of vegetation is the consequence of recovery from previous forest fires.
Changes associated with degradation of vegetation (yellow colors) stem from more recent forest fires and can be seen in greater detail in Window A of Figure 4 (transitions from forest or shrubs to grassland).Finally, the expansion of urban zones and infrastructure (red colors) mostly took place in the surroundings of Sabiñánigo, the center of the district with the largest population and industries.Window B of Figure 4 shows a golf course under construction.In the center of Figure 3, certain long patches can be seen, which have changed to artificial use as a result of the widening of the N-260 motorway.In the NW quadrant of the same figure, there are also small grasslands that have been converted into new ski slopes and other infrastructure associated with skiing in Formigal.
In the Guadarrama NP and its surroundings, most of the changes are concentrated in the south, in the region of Madrid and, especially, in the buffer and its SIZ (Figure 5).This is the most dynamic region from the demographic and socio-economic points of view, so it is the region that exerts the greatest pressure on the environment.In this study area, there were also similar processes, although of different intensity: increase in artificial areas, natural regeneration and degradation of vegetation (Table 4).Improved accessibility to Madrid has facilitated urban development of the towns to the south of the mountains (region of Madrid) and, also, although to a lesser extent, to the north (region of Castilla y León).Artificial areas increased by more than 5500 ha, for both residential (+1.19%) and industrial and commercial uses (+0.18%).Low-density, single-family housing predominates, mostly for use as holiday homes.These changes occurred in former forest areas (shrubs, grassland and forests) and on land occupied by heterogeneous agriculture.This phenomenon did not affect either the NP or the PPZ, but did affect the SIZ, and this amounts to one of the most worrying threats because of its proximity to the NP.In this study area, there were also similar processes, although of different intensity: increase in artificial areas, natural regeneration and degradation of vegetation (Table 4).Improved accessibility to Madrid has facilitated urban development of the towns to the south of the mountains (region of Madrid) and, also, although to a lesser extent, to the north (region of Castilla y León).Artificial areas increased by more than 5500 ha, for both residential (+1.19%) and industrial and commercial uses (+0.18%).Low-density, single-family housing predominates, mostly for use as holiday homes.These changes occurred in former forest areas (shrubs, grassland and forests) and on land occupied by heterogeneous agriculture.This phenomenon did not affect either the NP or the PPZ, but did affect the SIZ, and this amounts to one of the most worrying threats because of its proximity to the NP.There were many interchanges between the other land uses, but they all resulted in negative net change.On the one hand, there was an invasion of former agricultural zones by shrubs (0.52% of the study area).In most cases, this was marginal agricultural land with limited capacity for agricultural use because of important biophysical limitations (slopes, shallow useful soil or stony soil, risk of erosion, bioclimatic limitations, etc.).From the environmental point of view, there are two sides to this phenomenon.On the one hand, it has led to an increase in available fuel, and its extension has led to increased connectivity between forest masses, thus increasing the risk and propagation of forest fires.Conversely, it allows for greater carbon capture.In addition, the regeneration of vegetation is the natural evolution of ecosystems towards the climax.Indeed, it stimulates the structural connectivity between PAs, for example.
Forests have gained by more than 0.68% of the study area, as a result of incentives created by management plans for the PAs located within the study area, including the actual NP.However, the opposite phenomenon was stronger.During the same period, more than 0.87% of the study area occupied by forests was lost.Most of them are now occupied by shrubs and grassland.This transition mainly occurred in the NP, its PPZ and its SIZ.Many of these patches were deforested by fires during this same period [72].Although there are ample provisions of firefighters and fire-fighting equipment, fires continue to take place and to be propagated throughout this study area.An example was the fire in Abantos which affected 500 ha of the south slope of this mountain in 1999.The dense pine forests have been replaced by grassland and, more recently, by shrubs.

Simulated Scenarios of Change in 2030
Cramer's V test (Table A3) shows that, in the Ordesa NP, the variable that is most closely associated with all classes of land use is altitude (0.4247).Next come other variables that are useful for the simulation such as those relating to relief (slopes), climate (temperatures and precipitation), accessibility and, finally, changes in population density according to the censuses that are closest to the study period.However, we eliminated the slope variable because it is highly correlated with altitude.In the Guadarrama NP, access to and from Madrid is the factor that is most closely associated with all uses, especially the most intensive: urban, industrial and agricultural.On the other hand, slopes explain more of the changes in marginal agricultural and forestry uses.The distance to reservoirs and accessibility to roads are moderately associated with some land uses.The remaining variables are practically unrelated to land uses, so, since they are not likely to explain any significant change, they were not selected.
In the Ordesa NP, we made various adjustments and calibrations in order to generate the simulated map for 2006 and to compare it with the real map for 2006.Considering CLC90 as the base map of the series, the comparison results in the following statistics: Kappa Sim = 0.865, TransLoc = 0.957 and Transition = 0.904.The URB classes (Kappa Sim = 0.428) and AGR (Kappa Sim = 0.539) score the worst values, while FOR reaches the best values (Kappa Sim = 0.912).
Finally, we used two matrices for change probabilities in line with the changes occurring over recent decades and with the two simulated scenarios (Table A4).If the expected conditions prevail, the trend scenario will see few land use changes.However, the green scenario predicts that more agricultural land will be abandoned, with the development of natural ecosystems.In Figure 6, we show land use changes expected between 2006 and 2030, in a Green Scenario with Restrictions and Incentives (GS30-WRI).Again, we can see two opposing trends.On the one hand, there is an increase in the artificial area (reds), with an annual rate of change of +1.2870 (Table 5).If the expected conditions are met, urban expansion will take place in the surroundings of Sabañánigo, mainly on former agricultural land and also on scrubland.Other transitions to artificial areas can be expected in the areas that are closest to the main roads.On the other hand, natural regeneration of vegetation (green colors) is the natural evolution of ecosystems towards the climax.We can expect transitions of agricultural land (about 1000 ha) and of grassland (also about 1000 ha) towards shrubs.In ecological terms, the latter are those that most worry PA managers because focal habitats might be lost.Finally, patches of grassland and shrubs can be expected to evolve towards forest (about 2000 ha).
Environments 2017, 4, 79 13 of 30 agricultural land will be abandoned, with the development of natural ecosystems.In Figure 6, we show land use changes expected between 2006 and 2030, in a Green Scenario with Restrictions and Incentives (GS30-WRI).Again, we can see two opposing trends.On the one hand, there is an increase in the artificial area (reds), with an annual rate of change of +1.2870 (Table 5).If the expected conditions are met, urban expansion will take place in the surroundings of Sabañánigo, mainly on former agricultural land and also on scrubland.Other transitions to artificial areas can be expected in the areas that are closest to the main roads.On the other hand, natural regeneration of vegetation (green colors) is the natural evolution of ecosystems towards the climax.We can expect transitions of agricultural land (about 1000 ha) and of grassland (also about 1000 ha) towards shrubs.In ecological terms, the latter are those that most worry PA managers because focal habitats might be lost.Finally, patches of grassland and shrubs can be expected to evolve towards forest (about 2000 ha).In Figures 7 and 8, we show the differences between the two simulated change scenarios: GS30-WRI, in which restrictions and incentives are taken into account, and GS30, in which neither restrictions, nor incentives are considered.An example is shown in Windows A. If conditions are as expected, the green patches will be grassland in the GS30-WRI scenario, in line with the policy to promote the maintenance of pastures in PAs, while in GS30, they will probably be shrubs.The transition to shrubs goes against the management plan for the Sierra y Cañones de Guara Regional Park that appears in this Window.In Windows B, G30-WRI predicts a slight increase in the urban area within the town of Aínsa, because of incentives from the Government of Aragon [73], while GS30 simulates these patches as agricultural areas.Windows C represents the difference between the two models on the banks of the river Cinca, the Pineta reservoir and the A13B road.Any use change in these areas of public domain is prohibited.For this reason, GS30-WRI respects its current agricultural use, while GS30 predicts a transition to shrubs.
Environments 2017, 4, 79 14 of 30 In Figures 7 and 8, we show the differences between the two simulated change scenarios: GS30-WRI, in which restrictions and incentives are taken into account, and GS30, in which neither restrictions, nor incentives are considered.An example is shown in Windows A. If conditions are as expected, the green patches will be grassland in the GS30-WRI scenario, in line with the policy to promote the maintenance of pastures in PAs, while in GS30, they will probably be shrubs.The transition to shrubs goes against the management plan for the Sierra y Cañones de Guara Regional Park that appears in this Window.In Windows B, G30-WRI predicts a slight increase in the urban area within the town of Aínsa, because of incentives from the Government of Aragon [73], while GS30 simulates these patches as agricultural areas.Windows C represents the difference between the two models on the banks of the river Cinca, the Pineta reservoir and the A13B road.Any use change in these areas of public domain is prohibited.For this reason, GS30-WRI respects its current agricultural use, while GS30 predicts a transition to shrubs.Logically, in the Guadarrama NP, the changes expected between 2006 and the Trend Scenario (TS30) are similar to those that took place between 1990 and 2006, although the change patches are larger.In the Trend Scenario with Restrictions and Incentives (TS30-WRI), artificial areas will probably continue to increase, and new urban zones will mainly be located on the Madrid side of the buffer (Figure 9).The process of the regeneration of vegetation will affect all the zones.In the PPZ and on the northern side of the SIZ, transitions are predicted from agricultural areas to grassland and shrubs because of the abandonment of farming.In the SIZ, the herbaceous and woody vegetation is expected to progress towards forest.In the western buffer, forest patches are expected to become consolidated in response to the incentives offered by several PAs (Campo Azálvaro y Pinares de Peguerinos Special Protection Area (SPA) and the Voltoya River Basin).
In Figure 10, we show the differences between simulated land use maps for 2030 (TS30 vs. TS30-WRI).In the NP and PPZ, small patches of heterogeneous agriculture are expected to be transformed into shrubs and grassland as a consequence of the incentives offered under the NP management plan.In the buffer and in the SIZ, there are differences between the two simulations: in the first model, there are areas that could be urban, while in the second, they will be shrubs and grassland (around the A6 motorway and the Santillana reservoir and in other areas noted in Table A2b, where urbanization is restricted, so transitions from any use to urban are not possible).In the western buffer, forest patches are expected to become consolidated as a result of incentives offered by several PAs (Campo Azálvaro y Pinares de Peguerinos SPA and the Voltoya River Basin).Logically, in the Guadarrama NP, the changes expected between 2006 and the Trend Scenario (TS30) are similar to those that took place between 1990 and 2006, although the change patches are larger.In the Trend Scenario with Restrictions and Incentives (TS30-WRI), artificial areas will probably continue to increase, and new urban zones will mainly be located on the Madrid side of the buffer (Figure 9).The process of the regeneration of vegetation will affect all the zones.In the PPZ and on the northern side of the SIZ, transitions are predicted from agricultural areas to grassland and shrubs because of the abandonment of farming.In the SIZ, the herbaceous and woody vegetation is expected to progress towards forest.In the western buffer, forest patches are expected to become consolidated in response to the incentives offered by several PAs (Campo Azálvaro y Pinares de Peguerinos Special Protection Area (SPA) and the Voltoya River Basin).
In Figure 10, we show the differences between simulated land use maps for 2030 (TS30 vs. TS30-WRI).In the NP and PPZ, small patches of heterogeneous agriculture are expected to be transformed into shrubs and grassland as a consequence of the incentives offered under the NP management plan.In the buffer and in the SIZ, there are differences between the two simulations: in the first model, there are areas that could be urban, while in the second, they will be shrubs and grassland (around the A6 motorway and the Santillana reservoir and in other areas noted in Table A2b, where urbanization is restricted, so transitions from any use to urban are not possible).In the western buffer, forest patches are expected to become consolidated as a result of incentives offered by several PAs (Campo Azálvaro y Pinares de Peguerinos SPA and the Voltoya River Basin).

Fragmentation of Habitats of Interest
At the landscape scale, alpine and sub-alpine grasslands are hardly fragmented by shrubs (HFI = 1.84,Table 6) in the Ordesa NP.From the NP towards its periphery, there is a gradient of increasing fragmentation (1.21 in the buffer).If the conditions predicted in the GS30 scenario prevail, the fragmentation can be expected to increase slightly in all areas, except in the NP.There is also a clear difference between PA networks.NPs are more effective for maintaining grassland ecosystems than Regional Parks.Of special interest are Sites of Community Importance (SCIs) and SPAs, which see values very close to those of unprotected areas.This may be the result of a lack of planning and of efficient management.In Figure 11, we show an example of the retreat of grasslands seen between 1990 and 2006 and predicted between 2006 and 2030.It is clear that the core, edges, branches and bridges of grasslands are disappearing and being invaded by shrubs (in grey colors).This trend is creating increasingly serious ecological and environmental problems, against the wishes of naturalists and PA managers.

Fragmentation of Habitats of Interest
At the landscape scale, alpine and sub-alpine grasslands are hardly fragmented by shrubs (HFI = 1.84,Table 6) in the Ordesa NP.From the NP towards its periphery, there is a gradient of increasing fragmentation (1.21 in the buffer).If the conditions predicted in the GS30 scenario prevail, the fragmentation can be expected to increase slightly in all areas, except in the NP.Table 6.Expected evolution of the fragmentation index for habitats of interest (alpine and sub-alpine grasslands) in the Ordesa NP and its surroundings between 2006 and 2030, taking into account the Green Scenario (GS30).There is also a clear difference between PA networks.NPs are more effective for maintaining grassland ecosystems than Regional Parks.Of special interest are Sites of Community Importance (SCIs) and SPAs, which see values very close to those of unprotected areas.This may be the result of a lack of planning and of efficient management.In Figure 11, we show an example of the retreat of grasslands seen between 1990 and 2006 and predicted between 2006 and 2030.It is clear that the core, edges, branches and bridges of grasslands are disappearing and being invaded by shrubs (in grey colors).This trend is creating increasingly serious ecological and environmental problems, against the wishes of naturalists and PA managers.In the case of the Guadarrama NP, urban zones and artificial infrastructure have not fragmented the agricultural and natural ecosystems of the NP and of its PPZ (Table 7).They have done so slightly in the surroundings, but this process has had a greater impact on the SIZ, especially its southern slopes, than on the buffer.Urban sprawl and newly-built infrastructure (high-speed In the case of the Guadarrama NP, urban zones and artificial infrastructure have not fragmented the agricultural and natural ecosystems of the NP and of its PPZ (Table 7).They have done so slightly in the surroundings, but this process has had a greater impact on the SIZ, especially its southern slopes, than on the buffer.Urban sprawl and newly-built infrastructure (high-speed train and AP61 motorway) have fragmented ecosystems of interest in the SIZ (Figure 12), causing negative visual impacts and obstructing movement for terrestrial mammals.Although this land does not fall under their authority, the managers of the NPs should therefore carry out special monitoring of these threats and draw up corrective and protective measures together with the regional and local land planning authorities.train and AP61 motorway) have fragmented ecosystems of interest in the SIZ (Figure 12), causing negative visual impacts and obstructing movement for terrestrial mammals.Although this land does not fall under their authority, the managers of the NPs should therefore carry out special monitoring of these threats and draw up corrective and protective measures together with the regional and local land planning authorities.

Discussion
There have been few attempts at guiding future planning of PAs [2,60,73,74].PANORAMA [75] is an online platform that shows examples, in PAs of the world, of planning based on participation between local communities and managers.Our research is in line with this objective.The simulated scenarios and other variants may be a good starting-point for discussion and agreements between local communities and managers.

Discussion
There have been few attempts at guiding future planning of PAs [2,60,73,74].PANORAMA [75] is an online platform that shows examples, in PAs of the world, of planning based on participation between local communities and managers.Our research is in line with this objective.The simulated scenarios and other variants may be a good starting-point for discussion and agreements between local communities and managers.
From a methodological point of view, the main discussion revolves around the spatial and temporal dimensions of the data used as inputs in the simulations.In spite of errors [61,62], we consider that CLC is a standard source of information that is available at the European level, allowing studies to be replicated in other locations.In addition, it has a historical series from 1990.The main drawbacks are the scale, which offers limited detail for local studies and, recently, the change in method in the 2012 series.We explored the Spanish Land Occupation Information System (in Spanish SIOSE) [76], which has a larger scale (1:25,000), but its historical sequence is short (2005-2012), its legend is complex, and each polygon is defined by more than a single class [77].We also considered the possibility of using the new Spanish Forest Map on a scale of 1:50,000 [78], but its time sequence is also short (2009)(2010)(2011)(2012).It might be assumed that in order to calibrate the model better, and therefore achieve a better simulation, it would be necessary to use the longest historical sequence of data available.In this case, we could use the historical aerial photographs taken in 1945 to build a land use map for that year.Instead of using just two dates (for the start and the finish), Paegelow proposes that all available dates be used, in irregular periods to cover non-linear trends during calibration [79].However, Candau [80] and Clarke [81] show that the use of an excessive number of data amounts to an unjustified effort in processing, introduces various change trends that often counter each other, makes the simulation more difficult and also introduces sources of error in the model resulting from the different methods and sources used to draw up the maps in the series.
From the point of view of simulating change scenarios, other authors have compared different analytical techniques [49].We chose the ANN technique over others for several reasons.Operationally, it is easy to use, and the results are similar to those obtained using logistic regression techniques [31].In order to avoid conditioning the predictive model, some Cellular Automata (CA) models do not pay much attention to zoning [65], which contains legal restrictions and incentives, among other aspects.However, other CA models take into account the zoning criteria [82].
From the point of view of territorial management, this is a key aspect that we think should not be ignored.Our simulation models simultaneously consider the geographic [83], multivariable and multitemporal models.Our simulations do not take into account only variables related to the neighborhood factor, but also other biophysical and socioeconomic variables.We worked with static variables.We were only able to take into account current legislation, without considering possible changes in it.Additionally, we were unable to consider the construction of future infrastructure that might change accessibility.Our models do not predict changes in the metropolis or in its behavior.This would have required knowing infrastructure and urban plans and including them in the simulated models as dynamic variables.Such information, which may be unknown or politically sensitive, is not always available.Finally, in spite of the importance of forest fires for modelling the landscape and for many land use changes, they should not be included in simulations because we cannot know where fires are going to take place.It is only possible to know where there is a greater probability or risk of forest fires [84].
From the point of view of calibrating the models and validating the resulting maps, some authors show that the probability of being correct is fairly limited [85].In calibration, LCM does not contain as many changeable stochastic parameters as other techniques such as CA [68] or Conversion of Land Use and its Effects (CLUE) [41], which allow for a better fit.Another problem is to simulate a scenario with restrictions considering past trends since 2006, prior to the declaration of the Guadarrama NP and approval of its management plan.Finally, it is not possible to validate the results.Obviously, it is not yet possible to compare any of the simulated models with the real land use map for 2030.Some authors [86] do not validate their results for similar reasons, although they defend their usefulness as a starting point for discussion among stakeholders.Simulations do not provide a single solution for end users, but they do facilitate communication and debate among different stakeholder groups [81].Moreover, for simulations to be useful to end users, Barredo argues that models have to differentiate between a larger number of land uses, even distinguishing between land use and land cover [87].It also seems essential to include sectoral policies and plans in models, however complex they may be.Finally, Gómez-Delgado points to the difficulty of simulating the distribution of land uses of different types, based on different driving factors that compete simultaneously in a single territory [88].
The results obtained in the Guadarrama NP are in line with previous studies in nearby or similar study areas [62,65,89].The trend scenarios built are in line with the results obtained by other authors [60].These researchers modelled trend scenarios in 1260 PAs in the USA and concluded that the greatest threat for them and for their buffers is urban expansion.The main land use changes taking place in Ordesa NP are similar to the findings of other studies in the same area [34,36,90] and in other mountain systems on the planet [91,92].From the point of view of landscape structure, the results are similar to those obtained by other authors [93] from two different approaches: fragmentation versus connectivity.We confirmed differences in the habitat fragmentation rates among, and within, PA networks [31].
We also aim to increase awareness among those responsible for the Spanish network so that similar methods can be applied in other NPs.To achieve this objective, we have invited their managers to a Workshop on Monitoring and Evaluation of Sustainability in National Parks.The expert panel included managers of national parks and of other Pas, as well as regional managers of territorial planning, mayors and representatives of local action groups.We consulted their opinions to find out their assessment of the method used in this work and the feasibility of its implementation as a tool to help make spatial decisions.Globally, they appreciated the models as a starting point for the debate.

Conclusions
We can confirm that the cores of both NPs are not subject to significant LUCCs on an intermediate scale of analysis.This fact is a partial indicator of efficiency.However, their areas of influence and buffers are subject to pressure and threats that might affect the sustainability of the NPs in the future, especially those that are located close to large cities.Managers should carry out constant monitoring in order to minimize the impact of LUCCs and of visitors, in the form of soil sealing, fragmentation of natural habitats, lower water quality, increased risk of forest fires or the concentration of visitors in a few hotspots, such as visitor centers and the most frequently-used paths.Some of these are irreversible processes.

Figure 3 .
Figure 3. Land use-cover changes that took place between 1990 and 2006 in the Ordesa NP and its surroundings.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests; OTH = Others.

Figure 4 .
Figure 4. Details of land use changes that took place between 1990 and 2006 in the area of socioeconomic influence and the buffer of the Ordesa NP.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests; OTH = Others.

Figure 3 .
Figure 3. Land use-cover changes that took place between 1990 and 2006 in the Ordesa NP and its surroundings.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests; OTH = Others.

Figure 3 .
Figure 3. Land use-cover changes that took place between 1990 and 2006 in the Ordesa NP and its surroundings.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests; OTH = Others.

Figure 4 .
Figure 4. Details of land use changes that took place between 1990 and 2006 in the area of socioeconomic influence and the buffer of the Ordesa NP.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests; OTH = Others.

Figure 4 .
Figure 4. Details of land use changes that took place between 1990 and 2006 in the area of socioeconomic influence and the buffer of the Ordesa NP.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests; OTH = Others.

Figure 5 .
Figure 5. Land use-cover changes that took place between 1990 and 2006 in the Guadarrama NP and its surroundings.URB = Urban areas; IND = Industrial areas; AGR = Agricultural areas; HET = Heterogeneous agricultural areas; FOR = Forests; SHR-GRAS = Shrubs and Grasslands.

Figure 5 .
Figure 5. Land use-cover changes that took place between 1990 and 2006 in the Guadarrama NP and its surroundings.URB = Urban areas; IND = Industrial areas; AGR = Agricultural areas; HET = Heterogeneous agricultural areas; FOR = Forests; SHR-GRAS = Shrubs and Grasslands.

Figure 6 .
Figure 6.Simulated model of LUCCs between 2006 and 2030 in a Green Scenario with Restrictions and Incentives (GS30-WRI) in the Ordesa NP and its surroundings.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests.

Figure 6 .
Figure 6.Simulated model of LUCCs between 2006 and 2030 in a Green Scenario with Restrictions and Incentives (GS30-WRI) in the Ordesa NP and its surroundings.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests.

Figure 7 .
Figure 7. Expected differences in land use in the Ordesa NP and its surroundings according to the simulated scenarios for 2030: Green Scenario with Restrictions and Incentives (GS30-WRI), on the left of the legend, and Green Scenario without restrictions or incentives (GS30), on the right of the legend.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests.

Figure 7 .
Figure 7. Expected differences in land use in the Ordesa NP and its surroundings according to the simulated scenarios for 2030: Green Scenario with Restrictions and Incentives (GS30-WRI), on the left of the legend, and Green Scenario without restrictions or incentives (GS30), on the right of the legend.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests.

Figure 8 .
Figure 8. Large-scale Windows from Figure 7 (GS30-WRI, left of the legend, versus GS30, right of the legend): (Window A) incentive for maintaining grassland in the Regional Park of Sierra y Cañones de Guara; (Window B) incentive for urbanization in the town of Aínsa; (Window C) restriction on use change in the area of public domain beside water and along roads.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests.

Figure 8 .
Figure 8. Large-scale Windows from Figure 7 (GS30-WRI, left of the legend, versus GS30, right of the legend): (Window A) incentive for maintaining grassland in the Regional Park of Sierra y Cañones de Guara; (Window B) incentive for urbanization in the town of Aínsa; (Window C) restriction on use change in the area of public domain beside water and along roads.URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests.

Figure 9 .
Figure 9. Simulated model of LUCCs between 2006 and 2030 in a Trend Scenario with Restrictions and Incentives (TS30-WRI) in the Guadarrama NP and its surroundings.URB = Urban areas; IND = Industrial areas; AGR = Agricultural areas; HET = Heterogeneous agricultural areas; FOR = Forests; SHR-GRAS = Shrubs and Grasslands.

Figure 10 .
Figure 10.Expected differences in land use in the Guadarrama NP and its surroundings according to simulated scenarios for 2030: Trend Scenario without restrictions or incentives (TS30), on the left of the legend, and Trend Scenario with Restrictions and Incentives (TS30-WRI), on the right of the legend.URB = Urban areas; IND = Industrial areas; AGR = Agricultural areas; HET = Heterogeneous agricultural areas; FOR = Forests; SHR-GRAS = Shrubs and Grasslands.

Figure 9 . 30 Figure 9 .
Figure 9. Simulated model of LUCCs between 2006 and 2030 in a Trend Scenario with Restrictions and Incentives (TS30-WRI) in the Guadarrama NP and its surroundings.URB = Urban areas; IND = Industrial areas; AGR = Agricultural areas; HET = Heterogeneous agricultural areas; FOR = Forests; SHR-GRAS = Shrubs and Grasslands.

Figure 10 .
Figure 10.Expected differences in land use in the Guadarrama NP and its surroundings according to simulated scenarios for 2030: Trend Scenario without restrictions or incentives (TS30), on the left of the legend, and Trend Scenario with Restrictions and Incentives (TS30-WRI), on the right of the legend.URB = Urban areas; IND = Industrial areas; AGR = Agricultural areas; HET = Heterogeneous agricultural areas; FOR = Forests; SHR-GRAS = Shrubs and Grasslands.

Figure 10 .
Figure 10.Expected differences in land use in the Guadarrama NP and its surroundings according to simulated scenarios for 2030: Trend Scenario without restrictions or incentives (TS30), on the left of the legend, and Trend Scenario with Restrictions and Incentives (TS30-WRI), on the right of the legend.URB = Urban areas; IND = Industrial areas; AGR = Agricultural areas; HET = Heterogeneous agricultural areas; FOR = Forests; SHR-GRAS = Shrubs and Grasslands.
Park, PPZ = Peripheral Protection Zone, SIZ = Socioeconomic Influence Zone, RP = Regional Park; SPA = Special Protection Area, SCI = Sites of Community Importance.

Figure 11 .
Figure 11.Trend of the fragmentation of grasslands in a window of the Regional Park of Cañones y Sierra de Guara (SE of the Ordesa NP study area) in three years: 1990 (A), 2006 (B) and 2030 (C).

Figure 11 .
Figure 11.Trend of the fragmentation of grasslands in a window of the Regional Park of Cañones y Sierra de Guara (SE of the Ordesa NP study area) in three years: 1990 (A), 2006 (B) and 2030 (C).

Figure 12 .
Figure 12.Trend in fragmentation of semi-natural and natural ecosystems in a window for the SW of the Guadarrama NP and its surroundings between 1990 (left) and 2006 (right).The pink patches are urban and artificial (infrastructure) areas.

Figure 12 .
Figure 12.Trend in fragmentation of semi-natural and natural ecosystems in a window for the SW of the Guadarrama NP and its surroundings between 1990 (left) and 2006 (right).The pink patches are urban and artificial (infrastructure) areas.

Table 1 .
Summary of the main characteristics of the national parks analyzed.NP, National Park.

Table 1 .
Summary of the main characteristics of the national parks analyzed.NP, National Park.

Table 2 .
Summary of the characteristics and sources of the data used for LUCC analysis and simulation.

Table 3 .
Matrix and statistics for land use changes, in % of the total study area, between 1990 (rows) and 2006 (columns) in Ordesa NP and its surroundings.

Table 4 .
Matrix and statistics on land use change, in % of the total study area, between 1990 (rows) and 2006 (columns), in the Guadarrama NP and its surroundings.

Table 5 .
Matrix and statistics on expected land use changes, in % of the total study area, between 2006 (rows) and 2030 (columns) in the Ordesa NP and its surroundings, considering the Green Scenario with Restrictions and Incentives (GS30-WRI).URB = Urban areas; AGR = Agricultural areas; GRAS = Grasslands; SHR = Shrubs; FOR = Forests;

Table 5 .
Matrix and statistics on expected land use changes, in % of the total study area, between 2006 (rows) and 2030 (columns) in the Ordesa NP and its surroundings, considering the Green Scenario with Restrictions and Incentives (GS30-WRI).

Table 6 .
Expected evolution of the fragmentation index for habitats of interest (alpine and sub-alpine grasslands) in the Ordesa NP and its surroundings between 2006 and 2030, taking into account the Green Scenario (GS30).

Table 7 .
Trend in the fragmentation index of natural and semi-natural habitats in the Guadarrama NP and its surroundings between 1990 and 2006.= National Park, PPZ = Peripheral Protection Zone, SIZ = Socioeconomic Influence Zone. NP

Table 7 .
Trend in the fragmentation index of natural and semi-natural habitats in the Guadarrama NP and its surroundings between 1990 and 2006.= National Park, PPZ = Peripheral Protection Zone, SIZ = Socioeconomic Influence Zone. NP

Table A2 .
Land uses-land cover and transitions with restrictions or incentives in the Ordesa NP (a) and Guadarrama NP (b) and in their surroundings by management areas.

Table A3 .
Cramer's V test in the Ordesa NP (a) and Guadarrama NP (b) and their surroundings.
Note: Factors with greater explanatory power regarding land uses are in bold print (values > 0.15).