**Assessing Sustainable Rural Development Based on Ecosystem Services Vulnerability**

**Pascual Fernández Martínez <sup>1</sup> , Mónica de Castro-Pardo 2,\* , Víctor Martín Barroso <sup>1</sup> and João C. Azevedo <sup>3</sup>**


Received: 16 June 2020; Accepted: 6 July 2020; Published: 9 July 2020

**Abstract:** Sustainable Rural Development is essential to maintain active local communities and avoid depopulation and degradation of rural areas. Proper assessment of development in these territories is necessary to improve decision-making and to inform public policy, while ensuring biodiversity conservation and ecosystem services supply. Rural areas include high ecological value systems but the vulnerability of environmental components in development indicators has not been sufficiently pinpointed. The main objective of this work was to propose a new sustainable rural development composite indicator (nSRDI) while considering an environmental dimension indicator based on ecosystem services vulnerability and social and economic dimension indicators established using a sequentially Benefit of the Doubt-Data Envelopment Analysis (BoD-DEA) model. It aimed also to test effects of weighting methods on nSRDI. The composite indicator was applied to 10 regions (*comarcas*) in the Huesca province, Spain, producing a ranking of regions accordingly. The indicator was further tested through the analysis of the effect of an equal and optimum weighting method on scores and rankings of regions. Results showed substantial differences in nSRDI scores/rankings when vulnerability was added to the process, suggesting that the environmental dimension and the perspective from which it is conceived and applied matters when addressing sustainable rural development.

**Keywords:** sustainable rural development; regional composite indicators; vulnerability; ecosystem services; goal programming; analytic hierarchy process; data envelopment analysis; Spain

#### **1. Introduction**

Depopulation of rural areas is a major development challenge in Europe. Rural regions currently account for 28% of Europe's population (cities account for 40%; towns and suburbs account for 32%) [1] but, by 2050, more than 50% of the population of Europe will live in urban areas due to an expected increase of 24.1 million persons, as the population of rural regions is predicted to decrease by 7.9 million [2]. These demographic processes are not, however, homogeneous across Europe. The strength and persistence of rural depopulation trends in recent decades in countries in the south of Europe such as Portugal, Spain and Italy, has increased the number of rural regions of low population density in these countries. Until the 1990s, the major direct cause of depopulation was immigration to other countries (the Americas, Central Europe) and within-country rural-to-urban migration but since then, aging became the main cause of rural depopulation. The Spanish provinces of Aragon, Extremadura and Castilla-León, for example, have experienced levels of depopulation so severe that several villages are now uninhabited, while rates of aging and population loss in other provinces jeopardize their future [3]. The recent return of migrants and newcomers in rural areas, in particular during the 2010–2014 financial crisis and the current COVID-19 pandemic, is insufficient to mitigate ongoing depopulation.

In addition to social, economic and cultural effects, rural depopulation has negative consequences over nature conservation. In South European countries, one of the major effects of depopulation, while simultaneously affecting biodiversity, supply of ecosystem services and environmental quality, is an increase in wildfire hazards. Land abandonment related to depopulation changes the vertical and horizontal fuel structure in ecosystems and landscapes, thereby increasing fire hazards [4]. Abandoned landscapes tend to exhibit higher fire spread and intensity due to the spatial pattern of fuel load and to the reduction of fuel breaks (managed areas) in the landscape [5,6]. The consequences of fires of higher severity in abandoned landscapes becomes an even more serious issue when these landscapes contain systems of high ecological value which have been maintained through moderate levels of management intensity over periods of centuries to thousands of years. Several reports show a recent increase in forest fires in conservation areas. During the summer of 2009, approximately 30% of the burned areas in Europe were in Natura 2000 sites that were seriously affected and now face great challenges of recovering the original conditions [7].

Rural areas globally suffer from endemic poverty. International Fund for Agricultural Development— IFAD [8] estimated that over 60% of the poor will still be in rural areas, even in 2025. Rural poverty is perpetuated by the lack of access to essential assets like basic infrastructure, education or knowledge to access technologies and markets that could improve their productivity and income [9]. In Europe poverty is not exclusive to rural areas, but most of its poverty is found in rural areas [1]. Both poverty and degradation of high ecological value sites in rural areas need to be taken into account when rethinking rural development and defining the priority criteria to be considered in the development of tools to support rural regions and the corresponding improved assessment mechanisms in order to preserve biodiversity and increase welfare in local communities.

In Europe, the reversion of degradation processes, such as those mentioned above, and the development of rural areas have become important objectives in policy-making. Although support to low productivity agriculture dates back to 1972, strong efforts oriented towards rural development have been put forward by European institutions since early 1990s, assuming the most relevance in terms of financial support from Common Agricultural Policy (CAP) measures and the corresponding European Structural and Investment funds [10]. From a conservation perspective, European Union policies, such as CAP, have been progressively reformed, increasingly stressing goals directly and indirectly related to biodiversity and the environment. In early stages, CAP encouraged intensive farming and this caused the loss of habitats and species during the 20th century [11]. The first agri-environmental payments were introduced by the McSharry reform in 1992 [12]. However, it was not until the 2013 CAP reform that substantial changes were introduced. These were oriented towards conservation of nature, such as greening payments, to support sustainable farming, sustainable management of natural resources, climate action and balanced territorial development focused on rural areas with expected impacts on semi-natural habitats and wild species across Europe [13]. Other European Union initiatives, such as LEADER, INTERREG or LIFE programs have channeled important investments towards vulnerable regions of low productivity and low income agriculture, weak industry and weak services sectors, all of which are highly dependent on public policy [14,15]. Although these efforts have apparently been well directed and their effects positive [10], socioeconomic and conservation goals in EU policy have not always been easy to integrate. Conflicts between stakeholders with different interests in rural areas make it difficult to establish locally consensual and coordinated planning and management instruments [16] and this difficulty also affects European policies and their application. Despite efforts, integration of goals has not been sufficiently achieved [17], which requires methodological developments in particular to pinpoint and consider environmental components in decision-making processes.

The complexity and multi-functionality of rural socioecological systems, including latent and active conflicts among actors and stakeholders, make decision-making in rural areas difficult. In this context, composite indicators (CIs) can provide powerful tools to aid decision-making processes [18]. However, the way CIs are structured and built is relevant to these processes and choices concerning the right models to use depend chiefly on the desires and preferences of the analyst [19]. For this reason, it is important that the decision-maker has the possibility to choose from several indicator alternatives according to the most suitable method to apply in different decisional contexts. In this work, we aim to overcome challenges discussed early regarding the assessment of sustainable development in complex rural systems.

The objective of this paper is twofold: (i) to propose a new sustainable rural development index (nSRDI) considering vulnerability of ecosystem services as part of its environmental dimension and using a sequential BoD-DEA model for its development; and (ii) to test the effect of different weighting methods on sustainable rural development index scores and rankings resulting from its application.

The added value of this study lies on providing a composite indicator to improve assessment and decision-making processes in sustainable rural development policies, allowing ranking regions according to social and economic aspects while considering environmental value according to the vulnerability of ecosystem services.

#### **2. Background and Hypotheses**

Rural development is "the set of activities and actions of diverse actors that taken together leads to progress in rural areas" [20]. The term progress, in its diverse meanings, is in this definition key to understand the evolution of the concept and practice of rural development over time and can be used as a reference to classify the different composite indicators proposed so far. Rural development has changed considerably from its origins in early 1900s up to the 1990s, when it was guided mostly by profitability according to a technocratic and technical, exclusive, big corporation oriented model, to the late 20th century concept(s) directed to sustainability of agriculture and other activities and based on a holistic, inclusive and participatory local and territorial model [21,22]. Rural policy has changed accordingly and the conceptual framework supporting rural development has also been reviewed in the last two decades, in close connection with the development and implementation of the concept of sustainability in political and sectorial areas. Over time, rural policies have been changing from a unifunctional agricultural model focused on food production to a multifunctional agricultural model producing a range of private and public goods, (positive) environmental and cultural impacts, agreeable landscapes and quality and safe products [22]. This new model integrates agriculture and the environment.

Today, making agriculture sustainable is a global challenge and European institutions have provided scope for enhanced sustainability in instruments like the CAP. In the draft of the European Green Deal (EGD) one of the fields proposed for improvement was related to the insufficient set of indicators available and the development of indicators addressing nature conservation [23]. To align CAP with the principles of sustainability, multifunctionality and payments for public goods, 10 urgent action points were proposed in the EGD. Action 7, 8, 9 and 10 are related to the review of indicators, strengthening environmental monitoring, identification and addressing global impacts of the CAP as well as improving the governance of the CAP and its reforms [24]. Regarding the improvement of indicators, the EGD makes a particular mention to the relevance to maintaining conceptual references such as the Sustainable High Nature Value (HNV) Farming framework. One of the specific measures in Action 7 ("Revise the set of indicators") is to "reintroduce the HNV indicator". Sustainable High Nature Value (HNV) Farming is based on the importance that some rural/agriculture landscapes have for biodiversity conservation and identify, as key elements of sustainable HNV farming: socio-economic typology of HNV farmers, sustainable agricultural systems, aspects related with communities of NHV farmers as identity, motivation or social recognition and finally, the economic concept of profitability [25].

Numerous indicators have been proposed to assess "development" and "rural development". The World Bank and the OECD are among the international institutions that have dedicated efforts to measuring progress of rural areas [26] through the establishment and use of comprehensive systems of indicators. The OECD grouped a set of basic indicators in four development dimensions: Population and migration, Social well-being and equity, Economic structure and performance and Environment and sustainability [22]. The most important issues regarding rural development by the OECD tended to discard the Environmental and sustainability group and selected all the other three dimensions as key groups to define rural development indicators [22]. This has, however, changed with the adoption of new concepts regarding rural development, such as the OECD Rural Policy 3.0 framework, and the use of indicators such as the Green Development index [27].

The World Bank [28] uses a core set of indicators that captures a myriad of components of rural development and poverty. These indicators are classified in five dimensions: Basic data, Enabling environment for rural development, Broad based economic growth for rural poverty reduction, Natural resource management and biodiversity and Social well-being. This system includes a specific dimension related to natural resources.

Few works address specifically the design of indicators related to the environmental domain in the construction of rural development CIs that use vulnerability as part of environmental indicators. Environmental indicators have been defined in multiple ways and used from diverse approaches in the construction of sustainability CIs. The broad definition of sustainable development of the Brundtland Report in 1987, based on an intergenerational equity principle, and sustainability sciences based on the social dependence on natural resources but not exactly specifying ways of operationalization [29,30] have opened the field for a high variety of composite indicators developed and tested with the goal of measuring sustainability from different approaches [18,19,31–36]

The Ecosystem Service (ES) approach sets the foundations for a new way of assessing sustainability emphasizing the rational exploitation of the environment and resources vs the most strict non-use conservation idea [37]. This citation of [29] effectively describes the role of the ES science in dealing with sustainability: "In much of the world, conserving nature out of moral obligation is a luxury most simply cannot afford. Nevertheless, human well-being is intimately linked to the immediate environment and natural capital is a vital part of the economic base. In the face of a sea of poverty, demonstrating the ignored links between nature and elements of well-being safe drinking water, food, fuel, flood control, and aesthetic and cultural benefits that contribute to dignity and satisfaction, is the key to making conservation relevant and, if we are lucky, possible". In recent years ES-based indicators have been developed and applied to address sustainability. Mononen et al. [37] described the process of establishing a national ES indicator framework for Finland directed to social-ecological sustainability. Chen et al. [38] proposed an ecosystem service-based sustainability CI for the urban agglomeration of the Lake Biwa region in Japan based on 22 indicators. Díaz-Balteiro et al. [39] proposed a methodology to the dynamic aggregation of indicators of sustainable forest management based on six climate change scenarios considering five ESs: timber production, carbon sequestration, habitat and biodiversity conservation, recreation and game habitat quality in Central Spain. Chen et al. (2020) proposed a method to develop a CI to assess sustainability in eight ESs (soil retention, biodiversity maintenance, food provision, raw-material production, climate regulation, hydrological regulation, recreation services and landscape aesthetics) through value, vulnerability and spatial scale analyzing of the effect of different calculation methods on the index behavior.

Some studies remarked specifically on the importance of vulnerability of ecological components to assess sustainability [31,39–41] but only a few recent works consider specific measures of vulnerability to develop sustainability composite indicators. Li et al. [42] proposed a sustainability CI based on livelihood considering indicators of natural capital, such as land capital or drinking water quality, to develop an evaluation index system of livelihood sustainability in rural destinations in the Wuhan area in China. To identify weights of sub-indicators they used a degree of coupling between livelihood

and ecosystems. After that, they used a weighted summation to aggregate each livelihood dimension and the entropy method used to weight them.

From a theoretical perspective, all previous research has added relevant value to the background of CI in the context of sustainability and it can be assumed, from the literature, that multifunctionality as part of the rural development concept and its assessment requires consideration of environmental aspects when constructing rural development indicators. The vulnerability of these environmental aspects has not been included in this type of CIs.

From an operational/methodological perspective, it is also necessary to emphasize the great diversity of approaches and methods used in the development of sustainability CIs. Singh et al. [43] found 18 different schemes of weighting indicators and 16 different aggregation methods in 41 studies to construct sustainability CIs. Díaz-Balteiro et al. [39] and Giménez et al. [32] used Goal Programming methods to rank sustainability of forest plantations in Spain. They proposed four Goal Programming models to aggregate sustainability indicators considering preferences of decision makers using a pairwise survey and applied them to rank 30 industrial forest plantations. Castellani et al. [44] applied a Sustainable Performance Index to measure welfare and development at local scales under the framework of the European Charter for Sustainable Tourism in Protected Areas. This CI was defined as the sum of the values of 20 indicators in six groups: population, housing, services, economy and labor and finally environment and tourism. Caschili et al. [45] proposed the Composite Index of Rurality (CIR) to assess rurality in the region of Sardinia (Italy) firstly studying accessibility using an indicator constructed by a doubly constrained spatial interaction model and then proposing CIR to evaluate rurality in a regional setting employing multivariate analysis. CIR was obtained as an unweighted sum of descriptors of three pillars: demography, economics and settlement.

One of the most significant challenges in sustainable development planning is to obtain information economically, ensuring that it is thematically, spatially and temporally relevant to support policy analysis and decision-making [46]. The usefulness of CIs to improve the management of complex problems and support decision-making processes has been demonstrated in different contexts, particularly in cases related to environmental problems [18,47]. CIs based on sustainability and ecosystem services have been commonly used in the assessment of complex problems related to rural development [47]. However, the development of CIs is not straightforward. CIs development faces methodological challenges related to the process of aggregating heterogeneous information. The use of different aggregation methods can lead to very different outputs. Thereby, an erroneous decision regarding aggregation methods can fundamentally alter how CIs perform. To aggregate simple indicators, weights are often used as measures of perceived importance of each analyzed subgroup. However, due to the lack of available information about subgroup importance and unwillingness to prioritize one subgroup above another, it is common to use equal weight methods. Although this procedure is seen as being neutral, it is still a weighting decision [47].

In connection to the objectives of this research, and based on the previous literature review and the identified research gaps related to the addition of vulnerability of environmental issues in rural development assessments and the effects of weighting methods on the behavior of composite indices, two hypotheses were tested:


Testing these hypotheses can indicate whether vulnerability should be included in the assessment of rural development and whether weighting methods can substantially impact regional rural development rankings. The first hypothesis deals with a conceptual challenge regarding the role of environmental components in the assessment of sustainable rural development and the second hypothesis with a weighting methodological issue.

#### **3. Materials and Methods**

#### *3.1. Study Area*

This research was applied in Huesca, a province of Spain located in the autonomous community of Aragon in the northeast of the country. The province comprises 202 municipalities distributed over 10 *comarcas*: Alto Gallego, Bajo Cinca, Cinca Medio, Hoya of Huesca, La Jacetania, La Litera, La Ribagorza, Los Monegros, Sobrarbe and Somontano de Barbastro (Figure 1). Overall, population density (220,000 inhabitants, 15,626 km<sup>2</sup> ) is 14 inhabitants per km<sup>2</sup> but many municipalities show densities much lower than average. Absolute population per municipality ranges from 100 to 500 inhabitants in 48% of the area of the province. Only 13 municipalities out of 202 present population levels higher than 2000 inhabitants. In this study, only the rural municipalities with a population density lower than 150 inhabitants per km<sup>2</sup> following OECD criteria were addressed [48]. Moreover, only municipalities with less than 2000 inhabitants were considered to deal with the problem of definition of rural areas in some Spanish municipalities. Indicators were collected or calculated at the rural municipality level, later grouped into *comarcas*. The rural municipalities cover 1,236,976 ha, which represents 79.16% of the area of the province [49].

Although the majority of the area comprised by *comarcas* corresponds to rural municipalities, rural municipalities represent just a small proportion of the population of *comarcas*. Depopulation is a major socio-economic problem of these rural areas [50]. High and constant rates of depopulation have led numerous rural municipalities to extremely low demographic densities [3], currently less than 4 inhabitants/km<sup>2</sup> in many of these. A rural exodus in the 20th century and consequent land abandonment have contributed to the current landscape structure and functioning in these rural areas. One of the consequences of depopulation/abandonment was the decline of traditional land-use/land-cover systems and the expansion of shrubland and forest regeneration systems and the impoverishment of cultural landscapes, among others [50].

In general, Huesca can be characterized by a high degree of rurality and a high ecological value. The economy of the rural areas in the region is based on agriculture, cattle raising and nature tourism. Farmland in the rural area is around 661,600 ha, which represents 42.3% of the area of Huesca. Irrigated crops, such as alfalfa, corn and fruit tree orchards, are very important in La Litera, Cinca Medio and Bajo Cinca *comarcas* while in Somontano and Monegros vineyards, fruit tree orchards and forage crops are the dominant systems. Rural municipalities contain 63,460 ha of grasslands supporting extensive cattle and sheep grazing systems, i.e., 5.13% of their area [49]. In addition to grazing, grasslands provide a multitude of environmental services, such as carbon sequestration and storage or regulation of the water cycle [51]. Although sheep are more important in terms of the number of animals, in recent years they have declined in favor of cattle. In mountainous areas, cattle still persist extensively based in meadows. Today, local communities in rural areas are strongly dependent of tourism, as the main activity of families or as the secondary activity of farmers, breeders and small entrepreneurs whose activity is focused on local products.

Mountains dominate the landscape in Huesca to a large extent. Most of the area of the province is located in the Aragonese Pyrenees, to the north, and the Iberian System to the south. In the south of the province forests of conifers are abundant and small woodlands of *Juniperus sabina* and *Quercus ilex* emerge in semi-arid areas with crops and extensive meadows. Overall, forest and shrubland cover in the Huesca province is around 714,500 ha including 567,750 ha (45.90%) in the rural municipalities addressed in this study [49].

Regarding natural value, 25.24% of the area of the Huesca province has been classified as Sites of Community Importance (SCI) and 23.64% as Special Protection Areas (SPA) [49] under the Natura 2000 network, and nearly 9% of the province has also been classified as protected areas under national and regional conservation systems [48]. These figures are still more relevant when presented in terms of rural area of *comarcas*. For instance, 67.67% of the rural area in the *comarca* of La Jacetania is SCI, of which 44.66% is protected area, and 50.69% of the rural area of Sobrarbe is SCI, of which 30.06% is protected area [49]. One of the most important protected areas in the study area is the Ordesa and Monte Perdido National Park (OMPNP). Located in the Sobrarbe *comarca* (municipalities of Bielsa, Broto, Fanlo, Puértolas, Tella-Sin and Torla), the OMPNP was established in 1918 and expanded in 1982 to its present limits. The National Park comprises a core area of 15,608 ha and a buffer zone of 19,679 ha. This National Park includes the most singular and representative ecological values in the study area and for this reason it was selected as the reference area for the assessment of vulnerability of ecosystem services in this study.

**Figure 1.** *Comarcas* in Huesca and location of Huesca province in Spain and in the autonomous community of Aragon (**left**), Source: [52]; Protected areas in the Huesca province (**right**), Source: [49].

#### *3.2. Indicators*

− Initially, a set of 21 indicators distributed by the three sustainability dimensions (social, economic and environmental) were selected for this study (Table 1) based on the published data available for the Huesca province, taking into consideration the dimensions of sustainability and the use of sustainability indicators of sustainable rural development in the literature, following the extensive review in Section 2. These indicators were classified according to impacts on sustainable development in two forms: indicators with a positive impact (+) and indicators with negative impact (−). All indicators were positively related to sustainable development, with the exception of the Aging Index and Burned area that were considered to be negatively related to social and economic development indicators and some environmental indicators (I12, I13, I15, I16, I20 and I21) were collected from the Government of the Aragón database [49]. Indicator I14 was obtained from a review of terrestrial vertebrate species registered on the Terrestrial Vertebrate Database of the Ministry of Environment of the Government of Spain [53], I15 from the CSIC herbarium and Government of Aragon spatial databases [54], I18 from the CORINE Land Cover database in Aragon [49], and I19 was calculated using Equation (1):

$$ECr = 1 - \frac{\left(\frac{\sum\_{i=1}^{n} LS\_i \star 100}{\sum\_{i=1}^{n} ES\_i}\right)}{100} \tag{1}$$

సభ

ଵ

 

1 −

where *EC<sup>r</sup>* is the erosion control index in *comarca r*, *LS<sup>i</sup>* is the area (ha) with more than 25 t/ha/y of soil loss in municipality *I* and *ES<sup>i</sup>* is the erodible soil in municipality *i* when *comarca r* is formed by *n* municipalities. *LS* and *ES* were calculated from data of the National Inventory of Soil Erosion in Spain [55].

The selection of the environmental indicators was supported by the identification of key ecosystem services in the Ordesa y Monte Perdido National Park (OMPNP) by [49] based on a coarse correspondence established as follows: Cultivated terrestrial plants represented by I12; Reared animals by I13; Wild animals and Genetic material for animals, grouped in a single class (ESs provided by wild animals), by I14; Wild plants and Genetic material for plants, grouped in a single class (ESs provided by wild plants), by I15; Surface water and Ground water grouped in the class Water, by I16; Regulation of extreme events by I17; Bioremediation, Mediation of nuisances, Lifecycle maintenance, Water conditions and Atmospheric conditions, grouped in a single class (Regulation ESs provided by natural environments), by I18; Soil quality by I19; Physical, Intellectual and Spiritual interactions, grouped in the class Leisure of nature, by I20; and Non-use value (landscape) by I21. The correspondence between ecosystem services and the selected indicators is supported by ecosystem functions and services classification systems and indicators established for their assessment [56,57].


**Table 1.** Description of indicators used for *comarcas* in the Huesca province by dimension and type.

#### *3.3. Aggregation*

The new Sustainable Rural Development Index (nSRDI) proposed in this work was obtained in a sequential process. The first step consisted of aggregating individual indicators in independent dimension indicators (DI) for the social and economic dimensions of sustainability based on an efficiency approach using a Benefit of Doubt-Data Envelopment Analysis (BoD-DEA) model. For the environmental dimension, aggregation was based on a vulnerability approach. The second step consisted of aggregating social, economic and environmental dimension indicators into single composite index using a modified BoD-DEA model. Aggregation was conducted after a correlation analysis with the set of indicators in Table 1 was performed to look for possible linear relationships between pairs of

indicators [19]. When correlations were statistically significant, one of the indicators was removed from the analysis, reducing the dimensionality of the data. Calculations were made in SPSS v15.0 (SPSS Inc., Chicago, IL, USA) (correlation analysis) and Lingo v18.0 (Lindo Systems Inc., Chicago, IL, USA) (BoD-DEA models).

#### 3.3.1. Social and Economic Dimension Indicators

The construction of composite indicators requires a multidimensional approach to aggregate heterogeneous information in a structured manner. Data Envelopment Analysis (DEA) is used to optimize variables in complex scenarios based on efficiency. It is a linear programing technique used to assess a set of productive units using input and output variables in an uncertainty context, where the weights of these variables and the production function thatrelates them are unknown. There is conceptual similarity between that problem and the construction of CIs, in which quantitative sub-indicators are available but the actual knowledge of weights is not [58]. Therefore, DEA can be applied in the construction of CIs from the Benefit of the Doubt approach (BoD) [58]. The main objective of BoD models is to obtain the most suitable weights to each compared assessment of a decision unit from the available data for each unit [59]. Thus, BoD models assess the performance of each decision unit regarding other decision units, therefore allowing the definition of the most suitable weights for each unit. The only difference between BoD models and traditional DEA models is that the only fixed variables are output variables, considering only one dummy input variable that assumes the value 1 for each decision unit. This approach has been used frequently in the construction of CIs (e.g., [58,60]). In particular, Cherchye et al. [61] used a BoD model to construct sustainability CIs.

In the case of this research, the underlying idea is that a good relative performance of a particular *comarca* in one specific dimension indicator indicates that for that *comarca*, this dimension is relatively important. The model proposed here to aggregate a set of selected individual indicators into each of the sustainability dimension indicators is formulated as described in Equations (2)–(5):

$$\text{Dlc} = \max \sum\_{i=1}^{m} w\_{c,i} I\_{c,i} \tag{2}$$

s.t.,

$$\sum\_{i=1}^{m} w\_{c,i} I\_{j,i} \le 1 \tag{3}$$

(*n* constraints, one for each *comarca j*),

$$\frac{w\_{i}I\_{j,i}}{\sum\_{i=1}^{m}w\_{i}I\_{j,i}} \ge \alpha\_{i} \tag{4}$$

$$w\_{c,i} \ge 0\tag{5}$$

#### (*m* constraints, one for each indicator *i*),

where *j* = 1, 2, ......, *n* and *I* = 1, 2, ......, *m*, *DI<sup>c</sup>* is the dimension indicator of decision unit *c*, *wc*,*<sup>i</sup>* is the weight of the decision unit *c* regarding indicator *i*, *Ic,i* is the indicator *i* for each decision unit *c*, *Ij,i* is the indicator *i* for each *comarca j* and ∝*<sup>i</sup>* is a bound parameter that represents the contribution of indicator *i* to each *comarca j*, regarding all indicators. This constraint was added to improve the discriminatory power of the model and to avoid extreme results. The use of proportions has been successfully applied in DEA-models to avoid extreme scenarios, i.e., that all the relative weight can be assigned to a single CI, which would contribute exclusively to the overall performance value, while the other indicators would assume zero as their relative weight [58]. This process was applied to the social and economic dimensions, resulting in two independent dimension indicators.

#### 3.3.2. Vulnerability Assessment

Vulnerability can be defined as "the degree to which a system, subsystem, or system component is likely to experience harm due to exposure to a hazard, either a perturbation or a stress/stressor" [62]. In essence, vulnerability refers to the potential for loss deriving from natural or other hazards and

changes [63]. There is not a consensus regarding the best method to assess vulnerability but there is agreement on the need for using vulnerability to assess environmental components of socio-ecological systems [64]. The quantification of vulnerability is structured in 3 stages:


Let M = (*mij*)*ij* be a general matrix of positive numbers *mij* representing the importance of an item *i* over another item *j* given by a participant. There is a set of positive numbers *w<sup>l</sup>* . . . *wn*, such that *mij* = *wi wj* for every *i*, *j* = 1, . . . , *n*.

$$\text{Min } \sum\_{l} \left( n\_l^{(1)} + p\_l^{(1)} \right)^p + \sum\_{s} \left( n\_s^{(2)} + p\_s^{(2)} \right)^p + \sum\_{t} \left( n\_t^{(3)} + p\_t^{(3)} \right)^p \tag{6}$$

s.t.,

$$w\_{l\bar{l}} - m\_{l\bar{l}} + n\_{l}^{(1)} - p\_{l}^{(1)} = 0, \ l = 1, \ 2, \dots, \ n(n-1), \tag{7}$$

$$n w\_{i\bar{j}} w\_{j\bar{i}} + n\_s^{(2)} - p\_s^{(2)} = 1, \; s = 1, 2, \ldots, \frac{n(n-1)}{2}, \tag{8}$$

$$w\_{i\uparrow}w\_{j\bar{k}} - w\_{i\bar{k}} + n\_{\bar{t}}^{(3)} - p\_{\bar{t}}^{(3)} = 0, \quad \text{if } i = 1, 2, \dots, n(n-1)(n-2), \tag{9}$$

$$0.11 \le w\_{\rm ij} \le 9 \; \forall \; i, j. \tag{10}$$

where *n* (1) *l* and *p* (1) *l* are the negative and positive deviations of the goal, respectively, for the constraints that ensure the condition of similarity in the position l; *n* (2) *<sup>s</sup>* and *p* (2) *<sup>s</sup>* are the negative and positive deviations of the goal, respectively, for constraints that ensure the condition of reciprocity in the position *s*; and *n* (3) *t* and *p* (3) *t* are the negative and positive deviations of the goal, respectively, for constraints that ensure the condition of consistency at position *t*; *mij* are the components of the matrix M for each pair of criteria; *wij* are the components of matrix W formed by the weights that represent the most similar weights to the components of the original M matrix for each pair of criteria *ij*.

This model has already been successfully applied to correct inconsistencies in paired matrices for planning in protected areas [16]. After application of the GP model, consistent matrices were obtained that are as similar as possible to the original ones, while ensuring the conditions of similarity, consistency and reciprocity required by matrices built using pairwise comparisons.

(iii) Aggregation. After inconsistency correction, each matrix was normalized and aggregated into a single matrix using a geometric mean. Final weights were obtained using the eigenvalue method. These weights represent the relative vulnerability of each ecosystem service. Once vulnerability weights were quantified, environmental dimension indicators were generated as a weighted sum.

#### 3.3.3. The nSRDI Composite Indicator

The indicator nSRDI is the result of the aggregation of social, economic and environmental dimension indicators obtained in the previous steps of the process (Sections 3.3.1 and 3.3.2), using a modified BoD-DEA model (Equations (11)–(14)). Aggregation follows therefore a mixed efficiency-vulnerability approach, optimizing social and economic weights but fixing the environmental vulnerability dimension indicator (*EVDIC*). As a result, nSDRI was calculated for each *comarca*, making it possible to rank these territorial units in the entire province of Huesca, as follows:

$$m\text{SRDIc} = \max \sum\_{i=1}^{m} w\_{c,i} \text{DI}\_{c,i} + \text{EVDI}\_c \tag{11}$$

s.t.,

$$\sum\_{i=1}^{m} w\_{c,i} D I\_{j,i} \le 1 \tag{12}$$

$$\frac{w\_{\bar{l}}DI\_{j,\bar{i}}}{\sum\_{i=1}^{m}w\_{\bar{i}}DI\_{j,\bar{i}}} \ge \alpha\_{\bar{l}} \tag{13}$$

$$w\_{c,i} \ge 0\tag{14}$$

where *j* = 1, 2, ......, *n* and *i* = 1, 2, ......, *m*, *DIc,i* is the socioeconomic dimension indicator *i* for decision unit *c*, *EVDI<sup>C</sup>* is the environmental vulnerability dimension indicator for decision unit *c*, *DIj,i* is the dimension indicator *i* for *comarca j*, *wc*,*<sup>i</sup>* is the weight of the decision unit c regarding dimension indicator *i*, and ∝*<sup>i</sup>* is a bound parameter that represents the contribution of indicator *i* for *comarca j*, regarding all the indicators.

*DIj,i* are obtained through Equations (2)–(5). *EVDI<sup>C</sup>* is calculated as:

$$EVDI\_c = \sum\_{i=1}^{m} w\_i I\_{c,i} \tag{15}$$

where *w<sup>i</sup>* is the vulnerability weight related to each ecosystem service and *Ic,i* is the environmental indicator *i* for decision unit *c*.

In order to test the effect of adding vulnerability based indicators in nSRDI, the index was also calculated removing vulnerability from the assessment and calculated weights for indicators in Table 1, after removal of redundancy, using two approaches: optimal and equal. In the optimal approach, since vulnerability was not considered, the process described in Section 3.3.2 was not included in the calculation of the index and social, economic and environmental dimension indicators were all obtained from Equation (2). Therefore, the BoD-DEA model applied did not consider a priori fixed weights. In this case, the weights in all the dimensions (social, economic and environmental) are the most efficient in order to achieve the best score of the DI. In the equal approach, an average was

applied to aggregate dimension indicators. This is one of the most common aggregation methods used to construct CIs when weights are unknown.

Finally, an inter-*comarca* divergence analysis was conducted to analyze the sensitivity of nSRDI to the weighting method applied. This analysis was based on regional pairwise matrices of Euclidean distances in percentage (Equation (16)) between nSRDI values for the three methods (vulnerability, equal, optimal), resulting in three *n* × *n* diagonal and symmetric matrices when *n comarcas* were compared. Analyses were conducted in Excel 2010.

$$dst = \left(\sqrt{(\text{Cl}\_s - \text{Cl}\_l)^2}\right) \times 100\tag{16}$$

#### **4. Results and Discussion**

#### *4.1. Indicator Selection*

As a result of the correlation analysis, I1, I3, I5, I13, I14 and I20 were removed from the set of indicators due to high correlation with other indicators. I1 was correlated with I2 (−0.707; *p* < 0.05), I3 was correlated with I5 (0.935; *p* < 0.01) and I7 (0.746; *p* < 0.05) and I5 was correlated with I7 (0.863; *p* < 0.01). Regarding the environmental dimension, I12 was correlated with I13 (0.832; *p* < 0.01), I14 (−0.778; *p* < 0.01) and I18 (−0.720; *p* < 0.05), I13 was correlated with I14 (−0.947; *p* < 0.01), I18 (−0.766; *p* < 0.01), I20 (0.680; *p* < 0.05) and I21 (−0.641; *p* < 0.05), I14 was correlated with I15 (0.656; *p* < 0.05), I18 (0.637; *p* < 0.05) and I21 (0.684; *p* < 0.05), I18 was correlated with I20 (0.668; *p* < 0.05) and I20 was correlated with I21 (0.765; *p* < 0.05). Selected indicators were then normalized (Min-Max scaling) for each of the *comarcas* (Table 2).

#### *4.2. Vulnerability Weights*

As a result of the individual assessment of vulnerability of ecosystem services based on indicators, 64 pairwise comparison matrices were obtained. Of these, 32 were inconsistent and corrected with the GP model described in Equation (3) by recovering 50% of the information, after which weights were calculated for selected indicators. The most vulnerable ecosystem services were Wild plants and Genetic material for plants (ESs provided by wild plants) which correspond to the indicator Plant richness (27.62%) and Bioremediation, Mediation of nuisances, Lifecycle maintenance, Water conditions and Atmospheric conditions, grouped in class Regulation ESs provided by natural environments, represented by indicator Forest Cover (24.76%) while the less vulnerability were Cultivated terrestrial plants represented by indicator Agriculture cover (3.81%) (Table 3).


**Table 2.** Final indicators by *comarca* and dimension. Sign within parentheses indicates type of indicator: (+) when having a positive impact and (−) when having a negative impact on development. Maximum values are in bold.


**Table 3.** Environmental indicators and vulnerability weights calculated for Ordesa and Monte Perdido National Park based on expert opinion.

#### *4.3. E*ff*ect of Vulnerability on the Environmental Dimension*

Assessment of the environmental dimension indicator quantified according to the three different weighting methods (vulnerability, equal, optimal) indicates that the process of weighting affects the score of the dimension indicator (Figure 2) as well as the ranking of *comarcas* in Huesca. The application of the process involving vulnerability of ecosystem services ranks Sobrarbe (C9) first, followed by La Jacetania (C5) and La Ribagorza (C7) (Figure 2). Sobrarbe (C9) ranked first also for the other two methods but the order of the remaining regions differs among methods. In terms of scores, the weights assigned through the optimal approach resulted in indicator scores higher than when the other two methods were used. This happens because the optimal approach provides the most favorable values of the weights, while the equal and vulnerability methods share the weights for all the considered indicators. Equal weights assign the same relative importance to each indicator causing divergences between *comarcas* to depend only on the differences of indicator values. When weights are not equal, the relative importance of each indicator can change the score and ranking of *comarcas* above the original indicators. For this reason, the vulnerability approach provides higher divergences than equal weights (see Section 4.4). This effect is more relevant when vulnerability weights rely on few indicators, as is the case in this study. Optimal weights caused the highest differences among *comarcas* due to the benchmarking nature of the optimization process that assign the best punctuation to the most efficient indicators regarding all the others. As such, when the optimal approach was followed, the highest distance observed was between La Jacetania (C5) or Sobrarbe (C9), and Los Monegros (C8), that received the best score. The differences between best and worst scores were very small when equal weights were used. The vulnerability approach provided a balanced solution relative to the other two approaches and produced the largest distance between Sobrarbe and Somontano de Barbastro *comarcas*.

Although, apparently, equal weights produced scores closer to vulnerability than optimum weights, so the regional rankings based on this dimension indicator differ. Alto Gállego (C1) and Hoya de Huesca (C4) swapped positions 4 and 5. Sobrarbe (C9), La Jacetania (C5) and La Ribagorza (C7) *comarcas* were better ranked when vulnerability was added, while Bajo Cinca (C2), Somontano de Barbastro (C10) and La Litera (C6) have lower positions. This is due to the higher weights of indicators that represent high vulnerability regarding forest area and plant richness. When the optimum approach was used, the ranking changed dramatically as compared to the vulnerability approach, although the two best positions were maintained. Considering vulnerability, *comarcas* such as Somontano de Barbastro (C10) dropped from 3th to 9th in the rankings. Contrarily, La Ribagorza (C7) and Los Monegros (C8) climbed from 7th to 3th, and from the last to 6th, respectively. In the case of Sobrarbe, the score is so high that it is not affected by changes in the environmental dimension indicator.

**Figure 2.** Normalized environmental dimension indicator scores by *comarca* and corresponding ranking position according to the process of calculating weights: vulnerability based weights and equal and optimal approaches.

#### *4.4. Global Ranking Using Optimum, Equal and Vulnerability Weights*

ߙ ߙ Results of the composite indicator (nSRDI) for all sustainability dimensions considering vulnerability weights as part of the environmental dimension indicator showed La Jacetania (C5), Sobrarbe (C9) and La Ribagorza (C7) *comarcas* to be the regions of higher sustainable regional development in the Huesca province (Table 4). These results were obtained through the application of Equations (2)–(5) for the social and economic dimensions with the bound parameter α*i*= 0.05 (constraint 4) in a first round followed by a tie break (second round) of the *comarcas* with maximum scores that ranked first (score of 1) with the bound parameter set to α*i*= 0.2 to increase the discriminatory power of the model and to solve ties (indicated with asterisk in Table 4).


**Table 4.** Global new Sustainable Rural Development Index—nSRDI rankings of *comarcas* for weights calculated with the vulnerability, equal and optimal approaches.

\* Tied in the first round.

The effect of different weighting methods on nSRDI was assessed comparing *comarca* rankings obtained with the indicator under equal, optimal and vulnerability approaches (Table 4) and by the inter-*comarca* divergence matrices of nSRDI pairwise distances in the first round of the calculation, i.e., before tie-break (Figure 3) calculated with Equation (15), which helped us identifying the most relevant divergences among regions based on weighting method.

Divergences between *comarcas* were identified when vulnerability was added in the calculation of nSRDI. For instance, divergence between C8 and C5 and between C8 and C9 was 64.12%. This means that *comarca* Los Monegros (C8) presented the highest distance in score of nSRDI for these 2 regions when vulnerability was considered as part of the environmental dimension. However, as observed in the previous subsection, the strongest disparities among highest and lowest scores were observed when the optimum approach was applied (Figure 3). The highest distances observed (74.99%) were between C8 and C1, C5 and C7, *comarcas* in the highest positions in the ranking. With this method, the first round of calculations resulted in four ties for the top ranked *comarcas*. Moreover, all *comarcas*, with the exception of C2 and C8, had a score above 0.8, keeping a great distance from other *comarcas*. The lowest divergences occurred when equal weights were used. Here the highest divergence was 25.04% (C2 and C7). As already seen for the environmental indicator, the vulnerability approach provides a balanced solution between the equal and optimum approaches in the calculation of the composite index nSRDI. This means that, after the ties solved, this approach offers nSRDI results with enough discriminatory power to support a useful ranking of sustainable rural development for regions.

**Figure 3.** Pairwise *comarca* new Sustainable Rural Development Index—nSRDI divergence matrices by weighting approach (equal, optimal and vulnerability) based on distances between scores of the indicator. Dark red indicates high divergence and dark blue indicates low divergence. Values are in percentages.

#### *4.5. Hypotheses Assessment*

Results of this research support the two hypotheses presented initially (Section 2). The first, that the inclusion of vulnerability of ecosystem services in nSRDI affects *comarcas* sustainable rural development score and ranking position, is supported by results in Sections 4.3 and 4.4. This contributes to clarifying the importance of methods of describing and addressing quantitatively the environmental component in constructing composite sustainability indices. When vulnerability was added to the composite index calculation process, the scores of the environmental dimension indicator and nSRDI changed in comparison to the cases where the model excluded vulnerability. Moreover, rural sustainable development rankings also changed which suggest that the changes caused by the introduction of vulnerability were relevant. In particular, the most affected regions in the global ranking were *comarcas* with the highest ecological value in the Huesca province, as is the case of Sobrarbe (C9). This *comarca* was better ranked when vulnerability was considered. Also, La Jacetania (C5) scores were in the highest position when vulnerability was added, in comparison with equal weights, although it was also top with the optimum approach. The two *comarcas* with the highest percentage of SCI area were 67.17% (83,951 ha) in La Jacetania and 50.69% (97,318 ha) in Sobrarbe.

Economic development of rural areas with conservation areas lays mainly on tourism, often related to natural areas or sites with high nature value. Despite opportunity costs derived from land use and land cover limitations, it is likely that rural areas presenting a higher presence of conservation areas show higher levels of development which is related to the tourism industry. In fact, tourism is one of the main economic drivers of rural economies in Europe. Directly and indirectly it accounts for around 10% of European GDP and 20 million jobs [69]. In Spain, rural tourism drew 1.57 million tourists in 2016, growing 86.36% in 11 years [70,71]. In particular, nature-based tourism provides very important sources of income all over the world. Balmford et al. [72] estimated that the world's terrestrial protected areas together receive nearly 8 billion visits a year and these visits generate around US \$600 billion per year in direct in-country expenditure and US \$250 billion/year in consumer surplus. Vulnerability of ecosystem services in Huesca *comarcas* such as Sobrarbe (C9) or La Jacetania (C5), is particularly important because their main driver of development is tourism with a strongly dependency on nature conservation [73]. Increasing vulnerability jeopardizes the supply of ecosystem services over time and their contribution to local development in particular affects the cases of high ecological value sites.

Changes in the regional rankings due to the inclusion of vulnerability in the evaluation process can affect decisions over rural policies in particular when regions with high ecological value and more vulnerable ecosystem services do not receive enough institutional relevance as other rural areas. Adding vulnerability to the assessment of rural areas can improve the quality of decisions and guide resources management recognizing and boosting development of the rural areas that better provide social and economic development while conserving high ecological value and vulnerable sites.

The second hypothesis addressed, according to which the method of aggregation of sustainability components of rural development affects the ranking of regions according to their level of rural development and can emphasize divergences between regions, received support from results in Sections 4.3 and 4.4, contributing to guide the use of the most suitable decision-making tools regarding rural policies, in particular weighting methods. This hypothesis, more focused on a methodological issue, was put forth in order to analyze the effect of three weighting methods on indicators.

Despite their usefulness and increasing use in assessing complex and multidimensional problems, CIs still generate controversy. The dependency to a preliminary normalization stage and the disagreement among decision-makers on the specific weighting scheme used to aggregate indicators or dimension indicators, are the main problems that analysts face [61]. It is therefore relevant the information that can be extracted from the use of three weighting methods applied in this particular study.

Results showed that the optimum weighting approach produced the lowest discriminatory power among nSRDI scores resulting in three *comarcas* (La Jacetania, La Ribagorza and Alto Gállego) at the top of the corresponding ranking with the highest score. The vulnerability process also showed

limited discriminatory power with Sobrarbe and La Jacetania ranking first with the highest index score. The main weakness of the optimization models used for constructing CIs lay in the low discriminatory power and this could be a serious problem with basic BoD-DEA models because their capacity to rank units is weak. Fortunately, [74] found that adding constraints limiting the relative contribution of criteria can solve this problem. By imposing intuitive, widely accepted normative weight bounds, the discriminatory power increases. When these constraints are added, the value of the objective function decrease and only the decision units that have a minimum contribution in all the criteria can reach the maximum score. Thus, the bound parameter α*<sup>i</sup>* can be changed to add or remove flexibility to the model. [58,74] and [75] used this type of extended DEA models in the construction of CIs. In the present case study, the optimum and vulnerability weights were obtained setting the bound parameter α*i* to 0.05 for all the criteria (constraints in Equations (4) and (13)). When ties were observed, the model was reapplied over tied regions, increasing the parameter up to 0.2. The equal approach presented the highest discriminatory power. However, this method does not allow considering efficiency to build nSRDI, since it assigns the same weights to all criteria.

In the case of this study, the model proposed fills the information gap in the "right" set of weights by generating flexible BoD-weights for social and economic assessment and expert opinion AHP-GPweights for the environmental assessment for each *comarca*. The vulnerability approach initially provided ties for the two best positions in the ranking, but, as Figure 3 shows, provided a higher frequency of large distances (red cells) than the equal and optimum approaches. This shows wider regional distances among nSRDI values that represent robust positions in the ranking. Once ties were solved increasing the bound parameter, this method provided a complete ranking which presented the most efficient regions in terms of economic and social elements and vulnerability in terms of ecosystem services, simultaneously.

The assessment of both hypotheses emphasizes the relevance to make efforts to accurately define the relative importance of each assessed component in SRDIs.

The prominence of the resources directed to improve rural areas through the development of national and international policies requires rigorous and reliable tools to support decision-making processes. The challenge is to develop an integrated systems approach to sustainable rural development evaluation systems to support decision-processes that are effective in making sure public policy choices are well informed. The proposed model contributes to the achievement of this objective by providing a novel approach to the accurate assessment of sustainable rural development.

#### **5. Conclusions**

In this paper, a novel approach has been provided to accurately assess sustainable rural development with a potential strong contribution in the improvement of integrated system analysis and in the support of effective and informed decision-making processes and public policies.

The weighting method used in the construction of composite indicators was found to be not neutral. The proposed model presented a solution for problems of unknown information regarding weights by generating flexible BoD-weights for social and economic assessment and expert-opinion based AHP-GP-weights based on ecosystem services vulnerability for the environmental assessment.

When vulnerability of ecosystem services was added into the environmental DI, positions in the ranking changed dramatically and only the highest position was maintained. In this study, some *comarcas*, such as La Ribagorza and Los Monegros, obtained a much better score for the environmental DI based on vulnerability. Divergences between *comarcas* were identified when vulnerability was added in the calculation of nSRDI. The highest distances observed were between Los Monegros and Alto Gállego, and between La Jacetania and La Ribagorza. The vulnerability approach developed here provided a balanced solution in comparison to equal and optimum approaches in the calculation of the composite index nSRDI.

The inclusion of vulnerability of ecosystem services in nSRDI affected the *comarcas'* sustainable rural development score and ranking positions. Increasing vulnerability jeopardizes local development, especially in the areas with high ecological value sites, whose development depends greatly on nature-based tourism. Adding vulnerability in the assessment of rural areas can improve the quality of decisions and guide resources management promoting the regions that better address social and economic development while conserving high ecological value and vulnerable sites.

Based on the results obtained in this study, useful recommendations can be presented to institutions in order to improve the assessment of the sustainable rural development. One of the most relevant practical implications of this study is the finding that to obtain information regarding the analyzed indicators it is not only relevant, but also necessary to define accurately the relative importance of the assessed topics. This means that the concept of sustainability should be very well defined and structured in order to assign the proper importance to each of its dimensions: social, economic and environmental. Assigning the same relative importance to all the criteria is not a neutral decision and it may not correspond to the particular objectives of the assessment process or the political strategy, or even the interests of the analyst. Similarly, absolute flexibility regarding the relative importance of components may not be adequate when some criteria need to be limited.

This study showed some limitations related to the availability of quantitative information of vulnerability, which was solved by including expert opinion-based information in the model. This entails the need to carefully select the experts to be involved in the assessment, considering their diversity, know-how and the full knowledge of the analyzed problem. Moreover, future researcher lines should be oriented to the development of models that consider the feedback of stakeholders in different steps of the assessment process, using collaborative/participative methods, such as Delphi, in order to identify conflicts and achieve consensus in decision-making processes.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-445X/9/7/222/s1, Supplementary material S1: Questionnaire used in the Assessment of the Vulnerability of Ecosystem Services in The Ordesa and Monte Perdido National Park (originally in Spanish).

**Author Contributions:** Conceptualization, P.F.M. and M.d.C.-P.; methodology, P.F.M., M.d.C.-P. and V.M.B.; validation, J.C.A.; formal analysis, M.d.C.-P. and J.C.A.; investigation, M.d.C.-P. and V.M.B.; resources, M.d.C.-P. and V.M.B.; data curation, P.F.M., V.M.B., M.d.C.-P. and J.C.A.; writing—original draft preparation, M.d.C.-P.; writing—review and editing, M.d.C.-P. and J.C.A.; supervision, J.C.A.; project administration, P.F.M.; funding acquisition, P.F.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is partially funded by the Cátedra de Parques Nacionales UPM-URJC-UAH.

**Acknowledgments:** The authors would like to acknowledge "Organismo Autónomo de Parques Naturales" (OAPN), its director, the entire staff, and the stakeholders of the Ordesa and Monte Perdido National Park for their support in this endeavor. The authors would like also to acknowledge the valuable comments and recommendations made by four reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
